最近在利用GMT绘图过程中,为了丰富图件同时也为了增加美观,需要绘制Alaska区域的断层线。
在USGS下载到的断层文件是shp格式,GMT不能直接读取,因此需要我们进行格式转换,matlab就可以进行转换,貌似geoda也可以,不过没用过。
受制于matlab编程水平,虽然转换过程有些cumbersome,但最终也算是转换成功了。
1、首先用到的函数:shaperead,以qfaults.shp为例,
bbox = [-170 47;-130 65]; %框定经纬度坐标,只对该范围内的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即可。