经纬度和矩阵行列标的转换（基于MAP工具箱）及一个副高西伸脊点程序

斥鷃

好久好久好久好久好久好久没发帖了。
总想着在毕业论文搞完以后整理整理程序发几篇帖子，不过做完发现貌似没有现成能拿来发的，而后期也忙着给老师做东西，事情也就放下了，假期玩了这么多天想想再发篇帖子吧，就简单介绍一下MATLAB里面矩阵行列标和经纬度转换，顺便发个西伸脊点的程序吧~
以上废话，太久没发帖，废话偏多，见谅哈~

有关截取特点经纬度位置的函数，论坛里面沙颖凯的帖子《[分享]p坐标地形垂直剖面的MATLAB实现》一帖有给出小沙同学自编的GetLatLonInfo.m函数。建议初学者拿来看看，代码比较简单不过思路比较清楚。我这篇帖子主要从MATLAB里MAP工具箱的几个函数入手来实现截取特定经纬度，这样做的主要好处在于：
1.你熟悉这几个函数，尤其是投影矩阵的构造以后，MAP工具箱的函数都容易上手使用，方便处理shp数据；
2.出现需要多次转换经纬度和行列标的时候会更方便，尤其出现数次截取特定经纬度数据的情况。

进入正题：
一、数据的预处理（MATLAB中矩阵的绘图方式）
利用MATLAB对空间数据不加任何参数直接绘图时（如contour(data)），需要让行标代表纬度，列标代表经度，也就是说数据横着看过去，一行全是一个纬度的，竖着看过去，一列全是一个经度的，这个跟我们一般思维也是相似的。千万不要小看这点，很多数据读入时并不按这种方式排列。常见的例子就是grads的二进制数据了，第一维是经度，第二维是纬度，直接读取并reshape过来的话行标就是经度，列标就是纬度，需要转置才能正确用MATLAB画图。
所以建议在用MATLAB处理地理数据时，将其处理为行标代表纬度，列标代表经度的形式，既方便绘图，也方便在后期计算进行统一处理，这也是以下函数进行工作的一个基础，这是必须的一步。

二、投影矩阵
一个空间数据矩阵的每个数据对应着行列标，但是其实际对应的应该是经纬度位置。从行列标到经纬度，是两套不同的体系，如何实现它们之间的转换呢？MAP工具箱里面就是用投影矩阵R来实现，R是一个3×2的转换矩阵，假设row代表行标，col代表列标，lon代表经度，lat代表纬度，则有转换公式：

[lon lat] = [row col 1] * R

据此，我们可以利用这个矩阵来转换两种坐标。

补充内容（可略过）
细里讲，R的构造方法是（为了深入理解投影矩阵我们讲一下构造，但是具体构造用makerefmat函数实现就行，不需要刻意记住）：

0	Δlat
Δlon	0
lon0-Δlon	lat0-Δlat

其中Δlon，Δlat是经纬度的空间分辨率；lon0，lat0是西南角的经纬度（即经纬度各自的最小值）；这样，对于西南角那个点，应该是第1行第1列，上面提到的矩阵公式展开就是这个样子：

lon=1(row)*0+1(col)*Δlon+1*(lon0-Δlon)=lon0

lat=1(row)*Δlat+1(col)*0+1*(lat0-Δlat)=lat0

其他经纬度位置也是类似的，于是就把行列标转成了经纬度，实质只是做了一个线性变换而已，注意：Δlon和Δlat可以取负，这样就把升序变成了降序，取的lon0和lat0也要相应改变。

三、经纬度的截取和转换

1.makerefmat——构造投影矩阵

R = MAKEREFMAT(X11, Y11, DX, DY)%先经后纬

其中DX，DY是空间分辨率；X11， Y11是起始点的经纬。

举例：

如我取的矩阵数据是2°×2°分辨率，90-180°E，0-46°N的范围，则：

>> R=makerefmat(90,0,2,2)
R =
     0     2
     2     0
    88    -2

以上代码构造了一个起始经纬度为(90°E，0°)，分辨率为2°×2°的投影矩阵，利用它可以将数据矩阵中每个数据对应于地球上的经纬度位置，大家会问那终止经纬度如何表现呢？其实不用表现，数据矩阵的行列标最大值决定了最大可以取到什么位置。

在matlab中地图边界与掩膜(去掉边界外区域)的实现（基于shape文件）一文还有提及这个函数的另一个用法，更多用法详见help makerefmat，不继续做翻译工作了哈~

2.latlon2pix——将经纬度位置转换为行列标

为了方便截取矩阵数据中的特定经纬度的数据，我们必须弄清特定的经纬度在哪个行列标范围内，利用这个函数可以实现这个功能。
[ROW, COL] = LATLON2PIX(R,LAT,LON)%先行后列，先纬后经
其中R是投影矩阵，LAT、LON、ROW、COL分别是纬度、经度、行标及列标。

举例：
例如我想知道在上例2°×2°分辨率，90-180°E，0-46°N的data数据，但是要用97-107°E，21-30°N的数据，具体应该取哪行哪列？则：
>> [row col]=latlon2pix(R,[21 30],[97 107])
row =
   11.5000   16.0000
col =
    4.5000    9.5000

大家会发现出现了小数点，这是因为我给出的经纬度范围不在网格点上的原因，根据得到的数据，我们可以扩大选取data(11:16,4:10)，这样就取了96-108°E，20-30°N的数据，涵盖了截取范围，用于后续处理。这种数据的截取对减少运算量从而提高运算效率有很重要的意义，而不单单是出图有用！

3.pix2latlon——将行列标转换为经纬度位置
这一函数在确定某些特殊位置的时候起到很重要的作用，比如计算副高北界指数、西伸脊点指数等等。
[LAT, LON] = PIX2LATLON(R,ROW,COL)%先行后列，先纬后经
各符号代表与上例相同。

举例：
例如在上述2°×2°分辨率，90-180°E，0-46°N的data数据中，我找到了一个最大值位置，这个最大值位置在15行10列，则：
>> [lat lon]=pix2latlon(R,15,10)
lat =
    28
lon =
   108
即数据的最大值出现在108°E，28°N位置。

四、一个参考脚本例子
最后，贡献一个确定西太副高西伸脊点位置的函数，里面用了上面几种转换。
函数主要思路是截取西太副高西伸脊点定义区域，并确定区域内各纬度最西段5880线所在位置，最后确定西太副高西伸脊点的位置并输出该点经纬度，其经度位置便是我们常用的西太副高西伸脊点指数了。

bailing

感谢分享，不错噢

伊峰

顶顶顶顶顶顶顶顶顶顶顶顶顶顶顶顶顶顶顶

F平方

讲的非常清晰，受教～～

沙颖凯

好帖前排混脸熟 ( ′ ▽ ` )ﾉ 奇怪原来怎么没发现这贴

不过makerefmat最好还是用RasterSize吧，我试过自己的和用RasterSize得到的会有偏差。

斥鷃

沙颖凯 发表于 2014-9-1 12:34
好帖前排混脸熟 ( ′ ▽ ` )ﾉ 奇怪原来怎么没发现这贴

不过makerefmat最好还是用RasterSize吧， ...

可能主要因为选定特定经纬度范围来设置投影时候，要是直接用RasterSize方法的话，网格点的位置会有缩进，第一个点的经纬度位置算出来也并不在第一个点的实际位置，用X1，Y1这种模式做的话网格点的位置会是正确的，但是要多考虑一下边界上两个点的处理，常常会出现最后一个点或者第一个点算进去以后，构造画图的矩阵大小（用诸如meshgrid做的投影容易比矩阵大一些）和投影不匹配的问题。

新八不是废材男

  这个好贴怎么没在精华里面呢 @二爷名声在外

二爷名声在外

新八不是废材男 发表于 2014-9-3 09:29
这个好贴怎么没在精华里面呢 @二爷名声在外

因为整理的时候这个帖子还没发，之后一直没精力整理，希望理解一下，整理真是一件很累的事情

戒子1990

很好的帖子，就是关注度低点，改个名字好点

斥鷃

二爷名声在外 发表于 2014-9-3 09:39
因为整理的时候这个帖子还没发，之后一直没精力整理，希望理解一下，整理真是一件很累的事情

这个真的是，我以前自己也做过整理帖，实在是比较麻烦，你那个整理的特别细致了，辛苦了~