利用MATLAB进行shp文件转换并绘制断层线

利⽤MATLAB进⾏shp⽂件转换并绘制断层线
最近在利⽤GMT绘图过程中，为了丰富图件同时也为了增加美观，需要绘制Alaska区域的断层线。

在USGS下载到的断层⽂件是shp格式，GMT不能直接读取，因此需要我们进⾏格式转换，matlab就可以进⾏转换，貌似geoda也可以，不过没⽤过。

受制于matlab编程⽔平，虽然转换过程有些cumbersome，但最终也算是转换成功了。

1、⾸先⽤到的函数：shaperead，以qfaults.shp为例，
bbox = [-17047;-13065]; %框定经纬度坐标，只对该范围内的shp断层坐标进⾏读取，要尽可能的缩⼩该区域，否则参数太多容易卡掉⽽且faultname甚⾄可能都读取不出来s = shaperead('qfaults.shp', 'BoundingBox',bbox,'Attributes',{'Lon','Lat','faultname'})
s1=rmfield(s,'Geometry')； %删除掉结构体元素：Geometry
s1=rmfield(s1,'BoundingBox')；%删除掉结构体元素：BoundingBox
得到s1如下：是⼀个结构体。

2.由于s1是⼀个1421*1的结构体，不好操作，需要把它转成元胞数组cell
s2 = struct2cell(s1);
如下，变成了更熟悉的类似于矩阵的cell数组，s2{1,1}是⼀个1*43的数组，s2{1,2}是⼀个1*11的数组。

第⼀⾏代表经度，第⼆⾏代表纬度。

接下来需要把这些数组合并。

3.合并经纬度。

%合并经度
Lon = s2{1,1};
for i = 2:1412
Lon = [Lon,s2{1,i}];
end
Lon = Lon';
%合并纬度
Lat = s2{2,1};
for i = 2:1412
Lat = [Lat,s2{2,i}];
end
Lat = Lat';
%获得经纬度坐标
coordinate= [Lon,Lat];
4.最后利⽤metrix2txt转换为txt即可。

接下来根绝txt来绘制断层线。

5.我们打开txt后会发现有很多NaN，这也是我们⽤GMT绘图需要的！把txt⽤excle打开，把第⼀列中的NaN替换为>，然后把第⼆列中的NaN替换为空。

6.直接利⽤
gmt psxy AlaskaFault.txt -J$J -R$R -w0.2p,red,solid -K -O >> PS
注意，由于GMT5废⽌了-m，所以只要txt⽂件中有‘>'GMT就会⾃动绘制多条线段⽽不必在加‘-m’选项。