IDL编程学习之地图投影案例解析

15195775117

第一章

因为很多数据产品需要在地图上展示，如何把数据放到地图上(GIS化)就是个问题。
幸好IDL提供了投影转换功能，可以实现我们的美好愿望。
首先打开IDL自带的png文件
file=filepath('avhrr.png',subdirectory=['examples','data'])
file是png文件的路径：
FILE            STRING    = 'C:\Program Files\Exelis\IDL83\examples\data\avhrr.png'
读取这个png：
data=read_png(file,r,g,b)
data是个720列360行的字节型数组（图片都被读为字节型数组？不太明确）
DATA            BYTE      = Array[720, 360]
这里的r,g,b为：
R               BYTE      = Array[256]
G               BYTE      = Array[256]
B               BYTE      = Array[256]
这个结构有些匪夷所思！
表示红绿蓝通道的数组分别为：
r[data],g[data],b[data]
它们分别是：
<Expression>    BYTE      = Array[720, 360]
<Expression>    BYTE      = Array[720, 360]
<Expression>    BYTE      = Array[720, 360]
值得注意的是，fig=image(data)输出的图是黑白的；iimage,red=r[data],green=g[data],blue=b[data]输出的是彩图。
案例要求对红绿蓝三通道的数据进行“重采样”：
red0=rebin(r[data],360,180)
green0=rebin(g[data],360,180)
blue0=rebin(b[data],360,180)
与之前的像素相比，缩小到了之前的1/4。
画出图来：
iimage,red=red0,green=green0,blue=blue0,dimensions=[500,600],view_grid=[1,3]（这里view_grid相当于layout）。

15195775117

第二章.古蒂投影案例

创建(=初始化?)投影Interrupted Goode(古蒂等面积投影)
sMap=map_proj_init('Interrupted Goode')

这是，smap就是一个古蒂等面积投影了，意味着地球将按该投影展示。
对red0通道进行投影转换：
red1=map_proj_image(red0,map_structure=smap,mask=mask,uvrange=uvrange,xindex=xindex,yindex=yindex)

其中，map_structure用于获取地图投影方式，就是之前的smap；
mask(后一个)是掩膜区域文件，uvrange(后一个)是笛卡尔坐标范围，xIndex(后一个)和yIndex(后一个)是x和y的索引转换对应关系，这4个是被赋值的变量，与where()函数中的count类似，等号后边的完全可以写成其他的变量名。
MASK            BYTE      = Array[360, 180]
UVRANGE         DOUBLE    = Array[4]
XINDEX          LONG      = Array[360, 180]
YINDEX          LONG      = Array[360, 180]

利用xindex和yindex转换green0通道和blue0通道：
green1=map_proj_image(green0,xindex=xindex,yindex=yindex)
blue1=map_proj_image(blue0,xindex=xindex,yindex=yindex)

GREEN1          BYTE      = Array[360, 180]
BLUE1           BYTE      = Array[360, 180]

困惑！用红色通道生成的转换关系来转换绿蓝通道，是为了什么呢？
为什么不把绿蓝通道进行一样的转换呢？
以下将进行下对比。
以上这种情况绘图如下：
iimage,red=red1,green=green1,blue=blue1,alpha=mask*255b,/view_next(图放在下一个位置)

根据我的想法，
red1=map_proj_image(red0,map_structure=smap,mask=mask)
green1=map_proj_image(green0,map_structure=smap)
blue1=map_proj_image(blue0,map_structure=smap)
iimage,red=red1,green=green1,blue=blue1,alpha=mask*255b,/view_next


两者显然有区别，这是为什么呢？

15195775117

第三章.摩尔魏特投影案例

与第二章一个套路
mapstruct=map_proj_init('Mollweide',/gctp)（GCTP是某个投影集合包）
red2=map_proj_image(red1,uvrange,image_structure=smap,map_structure=mapstruct,$
  mask=mask,xindex=xindex2,yindex=yindex2)
注意到出现了2个投影,image_structure=smap(古蒂投影),map_structure=mapstruct(摩尔魏特投影)。
green2=map_proj_image(green1,xindex=xindex2,yindex=yindex2)
blue2=map_proj_image(blue1,xindex=xindex2,yindex=yindex2)
iimage,red=red2,green=green2,blue=blue2,alpha=mask*255b,/view_next

15195775117

使用过的代码见附件

sunzy025

咱们这里搞IDL确实是少