IDL与geotiff制作相同地理范围的地图

﹏ヽ暗。殇╰゛Y 2022-03-14 09:28 124阅读 0赞

pro geotiff_map_test
      compile_opt idl2
      forward_function TS_De_Anomaly
      forward_function TS_Anomaly
    ;超快速的geotiff
    ;  file = 'G:\sea_ice_study\smos_tb_2018_12_01_geo.tif'
    ;  CO=cgGeoMap(file, /Continents, /Grid,/Display,CCOLOR='RED')
    ;但是不能出图细调了，我就放弃了
      geoTiff_filename = 'G:\sea_ice_study\AMSR\AMSRAMSR_U2_L3_SeaIce6km_B01_20121201_geo_scale_horizontal_average.tif';'G:\sea_ice_study\smos_tb_2016_12_01_geo.tif';'G:\sea_ice_study\MODIS\mosaic_modis_bt_2012_12_01_new_band32.tif';;
      geoscale='G:\sea_ice_study\MODIS\mosaic_modis_bt_2012_12_01_new_band32.tif';'G:\sea_ice_study\AMSR\AMSRAMSR_U2_L3_SeaIce6km_B01_20121201_geo_scale_horizontal_average.tif'
    ;这里的MODIS范围最小，其他的都很大，所以我就直接用MODIS的地理范围了
      ; Open a display window.
      outpsfilename='G:\sea_ice_study\'
      filename='AMSRAMSR_U2_L3_SeaIce6km_B01_20121201_geo_scale_horizontal_average_c'
      cgSet_TTFont, 'Times'
       cgPS_Open,outpsfilename+filename+'.ps', Font=1
    ;这里是设置字体
      cgDisplay, 400,520,/FREE
      
      ; Set up the colors for the Geotiff image plot. Missing values are 0.
      cgLoadCT, 25, /Brewer, /Reverse, NColors=254, Bottom=1
    
      TVLCT, cgColor('white', /Triple), 0
    
      ; Get the image and map coordinate information at the same time.
      mapCoord = cgGeoMap(geoTiff_filename, Image=image)
    ;读取要plot的数据
      ; Missing values are identified as -10000.0 in this floating image.
      missingIndices = Where(image eq 0, missingCnt)
      IF missingCnt GT 0 THEN image[missingIndices] = !Values.F_NaN
    
      ; Find the image range and scale the data accordingly.
      minRange = 140;80;200(Floor(Min(image, /NaN) / 1e3) *  1e3 ) > 0
      maxRange = 280;Ceil(Max(image, /NaN) / 1e3) *  1e3
    
      scaledImage = cgImgScl(image,MinValue=minRange, MaxValue=maxRange)
      mapnormal=cgGeoMap(geoscale)
    ;这里读取的是MODIS的地理范围
      mapnormal -> GetProperty, XRange=xrange, YRange=yrange
      mapnormal -> SetProperty, XRange=[xrange[0], xrange[1]], $
        YRange=[yrange[0], yrange[1]], Position=[0.04848528, 0.21332673 , 0.94416315 ,  0.96667329],/Draw
    ;这里的position注意是要用MODIS的
    ;    cgMap_Grid, Map= mapnormal, /Label, LatDel=10, LonDelta=60,CHARSIZE=0.8, /cgGrid
        scaledImage=cgCliptoMap(scaledImage,[xrange[0],yrange[0],xrange[1],yrange[1]],Map=mapCoord)
    ;这里是为了把大的图像裁剪了
        mapCoord -> SetProperty, XRange=[xrange[0], xrange[1]], $
          YRange=[yrange[0], yrange[1]]
    ;改了地理范围了
    ;    0.00400000     0.104011      1.00000     0.945989
      ; Display the GeoTiff image with continental outlines and a map grid overlaid.
    ;  cgImage, scaledImage, MapCoord=mapCoord, /Keep_Aspect ,$;XRange=xrange, YRange=yrange, $ 
    ;    POSITION=[0.004,0.051,1.0,0.999], OPosition=opos;
        
         cgImage, scaledImage,Position=[ 0.04848528, 0.21332673 , 0.94416315 ,  0.96667329],MapCoord=mapCoord,  OPosition=opos, $;
              /Overplot,/Keep_Aspect;XRange=xrange, YRange=yrange,
              cgText, 0.7, 0.91, 'acquisition date: 1 Dec. 2012', ALIGNMENT=0.5,CHARSIZE=1.8, /NORMAL
    ;  cgImage, scaledImage,  Position=[0.00400000,0.104011,0.99,0.945989], OPosition=imagePos, $
    ;  XRange=xrange, YRange=yrange, /Overplot,/Keep_Aspect
       mapCoord -> SetProperty, XRange=xrange, YRange=yrange, Position=opos, /Draw
      ; 0.00000000     0.077448774       1.0000000      0.92255123
      ; Add map annotations.
    ;   mapCoord -> SetProperty,/OnImage, /Draw;   
      cgMap_Continents, Map= mapnormal;mapCoord;
      cgMap_Grid, Map= mapnormal, Label=1, LonDel=60, LatDel=10,CHARSIZE=0.8;,/cgGrid
    ;  cgMap_Grid, LatDel=30, LonDel=60, LineStyle=0, Color='Black', /Label, /No_Grid
      ; Add a color bar.
      cgColorbar, Range=[minRange, maxRange], NColors=254, Bottom=1, $
        Position=[0.085, 0.055+0.05, 0.925,  0.080+0.055], Title='AMSR2 89H Brightness Temperature (K)', $;'MODIS band32 Brightness Temperature (K)'SMOS 
        tLocation='top', OOB_Low=[0,50,173], OOB_High=[173,0,50], FONT=0, oob_factor=5,DIVISIONS=5
        cgPS_Close
        outpngfilename = outpsfilename +filename+'.png'
        cgPS2Raster, outpsfilename+filename+'.ps', outpngfilename, density=300, width=500
    end

这里我主要目前还比较不会的是position的定义，由于选用了/Keep\_Aspect，总是使得最后的图像plot的位置偏离原先定义的位置，导致上面的地理边界的矢量出现偏离，因此我都会运行一次得到opos就是图像实际的位置再重新把position写了。（希望有大神出来指导）