MATLAB绘制地图
1使用向量绘制地图
1.1绘制全球海岸线
向量数据可以表示一个地图。这种向量存在的形式是一系列的经纬度或投影坐标对，它们代表一个点集、一个线条或者多边形。例如，描绘出行政区域边界的点、公路系统、城市的中心或者以上三个集合放在一起，都可以被用于绘制地图。在这种表示中，地理数据以向量格式存在，以它为基础绘制的地图被称作向量地图。在地图工具箱环境中，向量数据包含一系列的有序的地理（经纬度）或投影（x，y）坐标对（又被称作二元组）。相继的数对被认为按顺序连接，间断点可以这样处理：构造另一个独立的向量，或者在点集的每个断点处插入一个分隔符号（通常是NaN，意思是Not-a-Number）。下面看一个向量地图数据的例子。
绘制全球海岸线
键入命令：
load coast 载入数据，数据是海岸线
whos 显示当前工作空间中的变量列表
屏幕显示：
Name Size Bytes Class Attributes
lat 9865x1 78920 double
long 9865x1 78920 double
键入命令：
axesm mercator
framem
plotm(lat,long)
[lat(1:20) long(1:20)]
可看到数据的地图
1.2查看向量坐标的数据
键入命令：
>> [lat(20:40) long(20:40)] 屏幕显示：
ans =
-80.9200 -146.5000
-80.6700 -145.5000
-80.3300 -148.0000
-80.0000 -150.0000
-79.6700 -152.5000
-79.2500 -155.0000
-78.8300 -157.0000
-78.7478 -157.2554
-78.6654 -157.5072
-78.5828 -157.7554
-78.5000 -158.0000
-78.4806 -157.6658
-78.4608 -157.3327
-78.4406 -157.0008
-78.4200 -156.6700
-78.5000 -154.5000
-78.1700 -154.5000
-78.1700 -154.5000
-78.0800 -156.6700
-77.8300 -158.0000
-77.5000 -158.3300
1.3修改颜色
这有没有提示你这些地点代表哪个大陆的海岸线？为了看见这些向量点代表的海岸线，键入下面的命令将其显示为黄色：
键入命令
plotm(lat(20:40),long(20:40),'y')
2点、线和多边形
它们代表着地理要素。向量地理数据可以用来代表点要素（例如城市和路标）、线要素（例如河流和公路）和面要素（或者说多边形，例如湖泊）。
在下一节，我们结合地理数据结构体来说明它们之间的区别
3 地理数据结构体
在上面的例子中，地理数据是以简单的向量的形式存在的。但是通常，我们会有
一系列的地理要素（例如一系列的河流），而每个要素都有自己的属性（例如每个河流都有自己的长度），为了方便对地理数据的操作，MATLAB将数据打包成地理数据结构体（geographic data structures）。
一个地理数据结构体是一个MATLAB结构体数组，其中每个元素对应着一个地理要素。每一个要素由坐标和属性来表示。一个具有地理坐标（也就是经纬度）的地理数据结构体被称作geostruct；具有地图坐标（投影后的x和y）的结构体被称作mapstruct。地理数据结构体只能用于向量数据，不能被用于处理栅格数据。
地理数据结构体通常在从shapefile中导入向量地理数据的时候产生。简单地说，地理地理数据结构体将地理要素的空间位置信息和属性数据打包在了一个变量中。
4 .绘制世界地图的例子
wm=worldmap('world'); 其投影方式以及经纬度范围自动适应全球的区域
land=shaperead('landareas','UseGeoCoords',true); 读取相应的shapefile文件到工作空间中,landareas中是polygon类型的数据
geoshow(wm,land,'FaceColor',[0.5 0.7 0.5]); 绘制地图
lakes=shaperead('worldlakes','UseGeoCoords',true); 读取世界河流数据geoshow(lakes,'FaceColor','blue'); 显示河流
rivers=shaperead('worldrivers','UseGeoCoords',true);
geoshow(rivers,'Color','blue');
cities=shaperead('worldcities','UseGeoCoords',true);
geoshow(cities,'Marker','.','Color','red');
5使用栅格数据绘制地图
5.1参考矩阵
栅格地理数据包括空间参考的数据网格和图像，它们在matlab中以矩阵的形式存在。尽管栅格地理数据看起来与其他实数矩阵一样，但是其不同在于它与空间位置相关。每一个数据像素占据地球上一个已知的区域块。那么如何把矩阵某个元素的行列位置与其对应的经纬度联系起来呢？
栅格地理数据包括空间参考的数据网格和图像，它们在matlab中以矩阵的形式存在。尽管栅格地理数据看起来与其他实数矩阵一样，但是其不同在于它与空间位置相关。每一个数据像素占据地球上一个已知的区域块。那么如何把矩阵某个元素的行列位置与其对应的经纬度联系起来呢
在以下的叙述中，我们规定，矩阵的每一行对应相同的纬度，每一列对应着相同
的经度。
如何根据矩阵中某个元素的行数和列数来计算其对应的经度呢？很自然地，用下面的方法来计算：
long = long11+(col-1)dlong = coldlong+long11-dlong??
类似地，用下式根据某个元素的行数和列数来计算其对应的纬度：
lat=lat11+(row-1)dlat=rowdlat+lat11-dlat??
把上面的两个式子写成矩阵的形式就是
0dla[long,lat]=[row,col,1]dlong=（0 dlat;dlong o;long11-dlong lat11-dlat）
我们把最后面的那个3*2矩阵叫做参考矩阵，也就是说只要知道了参考矩阵，就可以在矩阵元素的行列号和其对应的经纬度之间转换。把参考矩阵写成R，上式就简化为：
[long,lat]=[row,col,1]*R
在Mapping工具箱中，可以用函数makerefmat来获得参考矩阵。其调用:
R = makerefmat(lon11, lat11, dlong, dlat)
其中long11是矩阵第（1，1）个元素的经度，lat11是矩阵第（1，1）个元素的纬度。dlong是矩阵相邻两列间后一列相对于前一列的经度变化，dlat是矩阵相邻两行间后一行相对于前一行的纬度的变化。例如在matlab命令窗口输入：
R = makerefmat(5, 10, 3, 4)
函数makerefmat还可以以下方式调用。在命令窗口输入：
R = makerefmat('RasterSize', [20,20], 'Latlim', [30, 50],'Lonli m', [10, 30])
其中[20,20]指定了栅格矩阵的大小，[30, 50]指定了其纬度的范围，[10, 30]指定了其经度的范围。结果如下：
注意以这种方式调用函数makerefmat时，默认的经纬度方向是
矩阵从上到下对应从南到北，从左到右对应从东到西。也就是矩阵的第（1，1）元素是区域的西南角。
2.3 一个简单的栅格数据的绘图
想象一个极端粗糙的世界地图，其中每一个元素代表60度经纬度。这样地图的
矩阵大小就是3×6。
Z = [1 2 3 4 5 6; 7 8 9 10 11 12; 13 14 15 16 17 18];%地理数据3*6
R = georasterref('RasterSize', size(Z), ...
'Latlim', [-90 90], 'Lonlim', [-180 180]);%地理栅格数据参考对象(类) figure('Color','white')
ax = axesm('MapProjection', 'eqdcylin');%设定地图等距离圆柱投影方式
axis off%关闭本地坐标轴系统
setm(ax,'GLineStyle','--', 'Grid','on','Frame','on')%指定网格线形,绘
制frame框架
setm(ax,...
'MlabelLocation', 60,...%每隔60度绘制经度刻度标签
'PlabelLocation',[-30 30],...%只在指定值处绘制纬度刻度标签
'MeridianLabel','on',...%显示经度刻度标签
'ParallelLabel','on',...%显示纬度刻度标签
'MlineLocation',60,...%每隔60度绘制经度线
'PlineLocation',[-30 30],...%在指定值处绘制纬度线
'MLabelParallel','north' ...%将经度刻度标签放在北方,即上部
);
geoshow(Z, R, 'DisplayType', 'texturemap');%显示地理数据
colormap('autumn')
colorbar
对于georasterref对象,最主要的是要告诉它:
(a)栅格的大小:'RasterSize';
(b)栅格数据表示的地理范围: 'Latlim'和'Lonlim'.
否则对于原始数据Z,无法衍射到图形上的.
5.2 参考向量
当一个栅格元素的经度范围和纬度范围相同，也就是dlong和dlat相等时，可以使用更加紧凑的表示方法，即一个三元参考向量。一个参考向量定义了像素的大小和栅格数据的西北角：
refvec = [cells-per-degree, north-lat, west-lon]
在MAT文件中，这个变量通常被称作refvec或者maplegend。第一个元素cells-per-degree，描述了每一个栅格元素的角范围（假如每一个栅格元素覆盖5个经纬度，cells-per-degree就是0.2）。注意假如一个栅格元素的经度范围和纬度范围不一致，就不能用参考向量，而必须定义参考矩阵。
第二个元素，north-lat指定了栅格的最北端纬度，而第三个元素west-lon指定了栅格最西端的经度。换句话说，north-lat和west-lon指定了栅格的西北角。
注意矩阵的(1,1)元素始终是在栅格的西南角，但是当栅格或者图像用参考矩阵来描述的时候，这就不一定了。例如假如参考向量为[1 90 0]，则这个参考向量的解释是：
z 每一个数据网格元素代表1个经度和纬度；z 地图的北边界位于北纬90度（北极）；z 地图的西边界位于0度（本初子午线）。
所有栅格数据必须有一个参考向量或者参考矩阵，即使它们覆盖了整个地球。
二、下面的图使用参考向量(Reference vector).
对于参考向量,其实就是三个值,第一个值表示每个经纬度有多少个单元(cell),第2和第3个值表示数据中第一个点的经纬度坐标，一般是地图右上角的经纬度坐标:
refvec = [cells-per-degree,north-lat,west-lon]
注意经度和维度的网格划分是一样的时候，才可以使用参考向量！例如cells-per-degree=1,表示每隔1个经度和1个维度表示一个网格点。
clear;clc
load geoid;%地理栅格数据
load coast;%返回[lat,long],海岸线数据
figure; axesm robinson %建立地图坐标系,选择地图投影方式
geoshow(geoid,...
geoidrefvec,...%使用参考向量(Reference vector.),注意有些数据库里面是geoidlegend
'DisplayType','texturemap')
colorbar('southoutside')
geoshow(lat,long,'color','k')