本文共 7501 字，大约阅读时间需要 25 分钟。

Himawari-8的影像是个大圆盘，看起来就是这样：

要把这个大圆盘展开，需要对它进行几何校正。 IDL几何校正的思路如下： 1.利用行列号生成像元对应坐标（cols,lines->lon,lat） 2.利用坐标文件生成GLT 3.应用GLT文件对影像进行校正 几何校正的代码见下方 一些解释： 1.生成像元坐标时，cols,rows,LOFF,COFF,LAFC,CAFC参数根据分辨率不同有所区别。 2.生成GLT直接在ENVI里面做，注意像元大小（0.01或0.02）和旋转角度（0.0）设置，GLT生成好后对所有全云盘影像都可以用，一次生成，一直能用。

PRO Himawari8_LOC  COMPILE_OPT idl2  e=envi(/headless)  ;输入分辨率，10——1km，20——2km  resolution='10'    if resolution eq '20' then begin    rows=5500 & cols=5500    pixel_size=0.02    COFF=2750.5 & CFAC=20466275    LOFF=2750.5 & LFAC=20466275  endif else if resolution eq '10' then begin    rows=11000 & cols=11000    pixel_size=0.01    COFF=5500.5 & CFAC=40932549    LOFF=5500.5 & LFAC=40932549  endif  root_Dir=e.root_dir  demfile=filepath('GMTED2010.jp2',root_Dir=root_Dir,subdirectory=['data'])  demraster=e.OpenRaster(demfile)  demdata=demraster.getdata()  demspatial=[-180.00014,0.0083333,83.99986,-0.008333]  locdata=make_array(cols,rows,3,/float)  ;i,j从0开始校正后的影像偏差大概2个像元，这里从1开始找补一点偏差  for i=1,cols do begin    for j=1,rows do begin      ll=cl2ll(i,j,COFF,CFAC,LOFF,LFAC)      locdata[i-1,j-1,*]=ll      lon=ll[0]      if ll[0] le 360.0+demspatial[0] then begin        lon=ll[0]-360.0         endif      h=demdata[long((lon-demspatial[0])/demspatial[1]),long((ll[1]-demspatial[2])/demspatial[3])]      locdata[i,j,2]=h/1000.0    endfor  endfor  ;填补数据  fill_data=filldata(locdata=locdata)  ;导出保存在ENVI中生成GLT  fillraster=ENVIRaster(fill_data,URI='**loc_R'+resolution+'.dat')  fillraster.save  locraster.Close  fillraster.Close      e.CloseEND

将行列号转为经纬度的代码参照Himawari官网的代码，看截图：

翻译成IDL的代码是这样：

;行列号转经纬度function cl2ll,c,l,COFF,CFAC,LOFF,LFAC  COMPILE_OPT idl2  d2r=0.01745329252d  r2d=57.295779513d  scal=65536;2^(-16)=1/(2^16)  x=d2r*(c-COFF)*65536/CFAC  y=d2r*(l-LOFF)*65536/LFAC  ;'Distance from Earth's center to virtual satellite' =42164km  ;R(eq)^2 / R(pol)^2=1.006739501  ;'Coefficient for sd'=1737122264  rs= 42164.0  sd=sqrt((rs*rs*cos(x)*cos(x)*cos(y)*cos(y))-(cos(y)*cos(y)+ 1.006739501*sin(y)*sin(y))*1737122264)  sn=(rs*cos(x)*cos(y)-sd)/(cos(y)*cos(y)+ 1.006739501*sin(y)*sin(y))  s1=rs-sn*cos(x)*cos(y)  s2=sn*sin(x)*cos(y)  s3=-sn*sin(y)  sxy=sqrt(s1*s1+s2*s2)  lon=atan2(s1,s2)*r2d+140.7  lat=atan2(sxy,1.006739501*s3)*r2d  return,[lon,lat]   END

里面有用一个atan2函数，这个是另一个博文里写的：

;atan2函数function atan2,x,y ;此函数和matlab中的atan2输入的xy是相反的x=x*1d0 & y=y*1d0 absy = 0d & absx = 0d & val = 0dIF(x eq 0 and y eq 0)THEN BEGIN return, 0 ; ERROR! ENDIFabsy = abs(y) & absx = abs(x)IF(absy - absx eq absy)THEN BEGIN IF(y lt 0)THEN BEGIN return, -!pi/2 ENDIF ELSE BEGIN return, !pi/2 ENDELSE ENDIFIF(absx - absy eq absx)THEN BEGIN val = 0d0 ENDIF ELSE BEGIN val = atan(y/x) ENDELSEIF(x gt 0)THEN BEGIN return, val ENDIFIF(y lt 0)THEN BEGIN return, val-!pi ENDIFreturn, val+!piEND

计算出来的经纬度是个大圆盘，下面看图，生成GLT需要把圆盘边角全部填满，使数组没有无效值（0/NAN）和错值（经度大于360.0，纬度超出[-90.0,90.0]等）。

这里我定义填补后的经度范围为[60,222]，纬度填补后的范围是[-90.0,90.0]，填补代码如下：

;填补数据function filldata,locdata=locdata    cols=(locdata.dim)[0]  rows=(locdata.dim)[1]  ;填补经度数据  ;经度  londata=locdata[*,*,0]  ;中间列  lon0=londata[fix(cols/2),*]  lon0idx=where(finite(lon0))  ;缺值补值  lon0top=lon0idx.min()  lon0bottom=lon0idx.max()  londata[fix(cols/2),*]=[[make_array(1,lon0top,value=lon0[lon0top])],[rotate(lon0[lon0top:lon0bottom],1)],$    [make_array(1,rows-lon0bottom-1,value=lon0[lon0bottom])]]    ;中间行  lon1=londata[*,fix(rows/2)]  lon1idx=where(finite(lon1))  ;缺值位置  lon1left=lon1idx.min()  lon1right=lon1idx.max()  ;等差数列补值  lon1left_ps=(lon1[lon1left]-60.0)/lon1left  lon1right_ps=(222.0-lon1[lon1right])/(cols-1-lon1right)  print,lon1left_ps,lon1right_ps  ;缺处补值  londata[*,fix(rows/2)]=[lon1[lon1left]-(lon1left-findgen(lon1left))*lon1left_ps,$    lon1[lon1left:lon1right],lon1[lon1right]+(findgen(cols-lon1right-1)+1)*lon1right_ps]  ;纬度补值  latdata=locdata[*,*,1]  ;中间列  lat0=latdata[fix(cols/2),*]  lat0idx=where(finite(lat0))  lat0top=lat0idx.min()  lat0bottom=lat0idx.max()  lat0top_ps=(90.0-lat0[lat0top])/lat0top  lat0bottom_ps=(lat0[lat0bottom]+90.0)/(rows-1-lat0bottom)  print,lat0top_ps,lat0bottom_ps  latdata[fix(cols/2),*]=[[rotate(lat0[lat0top]+(lat0top-findgen(lat0top))*lat0top_ps,1)],$    [rotate(lat0[lat0top:lat0bottom],1)],[rotate(lat0[lat0bottom]-(findgen(cols-lat0bottom-1)+1)*lat0bottom_ps,1)]]     ;中间行  lat1=latdata[*,fix(rows/2)]  lat1idx=where(finite(lat1))  lat1left=lat1idx.min()  lat1right=lat1idx.max()  latdata[*,fix(rows/2)]=[make_array(lat1left,1,value=lat1[lat1left]),lat1[lat1left:lat1right],make_array(cols-lon1right-1,1,value=lat1[lat1right])]  ;全部补值     for i=0,cols-1 do begin            lonidx=where(finite(londata[i,*]))      lontop=lonidx.min()    lonbottom=lonidx.max()      if (lontop gt 0) then begin      londata[i,0:lontop-1]=rebin([londata[i,lontop]],1,lontop)    endif     if (lonbottom lt rows-1) then begin             londata[i,lonbottom+1:-1]=rebin([londata[i,lonbottom]],1,rows-lonbottom-1)     endif      latidx=where(finite(latdata[*,i]))    latleft=latidx.min()    latright=latidx.max()    if (latleft gt 0) then begin             latdata[0:latleft-1,i]=rebin([latdata[latleft,i]],latleft,1)    endif    if (latright lt cols-1) then begin           latdata[latright+1:-1,i]=rebin([latdata[latright,i]],cols-latright-1,1)    endif          endfor   locdata[*,*,0]=londata  locdata[*,*,1]=latdata    return,locdata  END

填上边角的经纬度文件长这样：

在ENVI中生成GLT文件：

注意pixel_size和rotation改好，rotation是0，pixel_size根据分辨率确定，2KM填0.02，1KM填0.01 生成的GLT长这样： 给影像应用GLT进行校正——ENVI里进行 也可以用IDL代码：

;应用GLT校正e=envi(/headless)inputraster=e.openraster('*****H8全云盘影像')fid=ENVIRasterToFID(inputraster)out_name='***给投影完了的文件起个名'#不想起名就用e.GetTemporaryFilename()ENVI_DOIT,'ENVI_GEOREF_FROM_GLT_DOIT',FID=fid,GLT_DIMS=glt_dims,GLT_FID=glt_id,OUT_NAME=out_name,pos=lindgen(bands)

校正后的影像长这样：

小小的偏移(i,j循环从0开始) 更小的偏移（i,j循环从1开始），找补回来一点了： 到这里影像投影完成了，再来美化一点点就是把它裁剪一下了：

outraster=e.OpenRaster(out_name)#out_name就是上面GLT校正完的影像SpatialRef = outraster.SPATIALREFSpatialRef.ConvertLonLatToMap, 80.0, 60.0, MapX, MapYSpatialRef.ConvertLonLatToMap, 200.0, -60.0, MapX2, MapY2subraster = outraster.Subset(SPATIALREF=SpatialRef,SUB_RECT=[MapX, MapY2, MapX2, MapY])

裁剪完输出一下：

subraster .Export,'****给个输出路径','ENVI'

到这里影像的投影就完事了，没别的需求的就不用往下看了。

. . . 别的需求就是生成角度数据集，具体就是：卫星天顶角、方位角，太阳天顶角、方位角。 但是这里都不讲这些，计算卫星天顶角、方位角看这： https://blog.csdn.net/qq_33339770/article/details/103023725 计算太阳天顶角、方位角看这： https://blog.csdn.net/qq_33339770/article/details/104710136 这里要对生成的坐标文件做一点处理，才能应用在角度数据集生成里，具体的操作就是：重投影+裁剪。 执行的步骤是和对影像的操作是一样的，先应用GLT投影成这样： 然后用上面的裁剪代码把它裁了，使同分辨率的坐标文件和同分辨率的影像完全贴合对应，这样就能作为角度数据集生成的输入文件来计算各个位置的卫星/太阳角度了。 . . . 一点废话 Himawari8官网还给出了经纬度反算行列号的代码，看图： 我也翻译成了IDL代码：

;经纬度转行列号function ll2cl,lon,lat,COFF,CFAC,LOFF,LFAC  COMPILE_OPT idl2  d2r=0.01745329252d  r2d=57.295779513d  scal=1d0/65536;2^(-16)=1/(2^16)  if lon gt 180.0 then begin    lon -= 360.0  endif  if lon lt -180.0 then begin    lon += 360.0  endif  ;phi = arctan( (Rpol^2)/(Req^2) * tan(lat) )  ;;rpol^2/req^2 =0.993305616  phi=atan(0.993305616 * tan(lat*d2r))  ;Re = (Rpol) / sqrt( 1 - (Req^2 - Rpol^2) / Req^2 * cos^2(phi) )  ;;earth's polar radius=6356.7523km  Re=6356.7523/sqrt(1- 0.00669438444*cos(phi)*cos(phi))  ;'Distance from Earth's center to virtual satellite' =42164km  ;sub_lon=140.7  r1=42164-Re*cos(phi)*cos(lon*d2r-140.7*d2r)  r2=-Re*cos(phi)*sin(lon*d2r-140.7*d2r)  r3=Re*sin(phi)  if r1*(r1-42164)+(r2*r2)+(r3*r3) le 0 then begin    rn=sqrt(r1*r1+r2*r2+r3*r3)    x=atan2(r1,-r2)*r2d    y=asin(-r3/rn)*r2d    c=long(COFF + x * CFAC * scal)    l=long(LOFF + y * LFAC * scal)    return,[c,l]  endifEND

但是我用代码试了，反算不回去ORZ，有0.2到0.5的差。

转载地址：https://blog.csdn.net/qq_33339770/article/details/102957857 如侵犯您的版权，请留言回复原文章的地址，我们会给您删除此文章，给您带来不便请您谅解！