GSHHS coastline to Surfer bln

无径之林

此贴转自：http://blog.sina.com.cn/s/blog_65f2ef3d0100s5q5.html

1.利用GSHHS自带的matlab程序获取想要的区域经纬度。如：

%  This Matlab script extracts coastline data from GSHHS database.

job='seagrid';            % Prepare coastlines for SeaGrid tool
%job='ploting';            % Prepare coastlines for NCAR ploting programs

database='full';          % Full resolution database
�tabase='high';          % High resolution database
�tabase='intermediate';  % Intermediate resolution database
�tabase='low';           % Low resolution database
�tabase='crude';         % crude resolution database

switch job,
  case 'seagrid'
    Oname='sb_coast.mat';
  case 'ploting'
    Oname='ch_coast.cst';
end,

switch database,
  case 'full'
    Cname='gshhs_f.b';
    name='gshhs_f.b';
  case 'high'
    Cname='gshhs_h.b';
    name='gshhs_h.b';
  case 'intermediate'
    Cname='gshhs_i.b';
    name='gshhs_i.b';
  case 'low'
    Cname='gshhs_l.b';
    name='gshhs_l.b';
  case 'crude'
    Cname='gshhs_c.b';
    name='gshhs_c.b';
end,

Llon=119;               % Left   corner longitude
Rlon=125;              % Right  corner longitude
Blat=27;                 % Bottom corner latitude
Tlat=37;                % Top    corner latitude

spval=9999.0;               % Special value

%----------------------------------------------------------------------------
%  Extract coastlines from GSHHS database.
%----------------------------------------------------------------------------

disp(['Reading GSHHS database: ',name]);
[C]=r_gshhs(Llon,Rlon,Blat,Tlat,Cname);

disp(['Processing read coastline data']);
switch job,
  case 'seagrid'
    [CC]=x_gshhs(Llon,Rlon,Blat,Tlat,C,'patch');
  case 'ploting'
    [C]=x_gshhs(Llon,Rlon,Blat,Tlat,C,'on');
end,

%----------------------------------------------------------------------------
%  Save extrated coastlines.
%----------------------------------------------------------------------------

lon=C.lon;
lat=C.lat;

switch job,
  case 'seagrid'
    save(Oname,'lon','lat');
  case 'ploting'
    x=lon;
    y=lat;
    ind=find(isnan(x));
    if (~isempty(ind)),
      x(ind)=C.type;
      y(ind)=spval;
    end,
    fid=fopen(Oname,'w');
    if (fid ~= -1),
      for i=1:length(x),
        fprintf(fid,'.6f  .6f\n',y(i),x(i));
      end,
      fclose(fid);
    end,
end,


2. 再将其转化为surfer的bln文件。代码：
clear all
% this is for GSHHS getted coast transformation to surfer ones!

% the yzcoast_coast.mat is getted from r_gshhhs.m function in the extract.m

load sb_coast.mat
lon(find(lon<119))=119;
lon(find(lon>125))=125;

lat(find(lat<27))=27;

lat(find(lat>37))=37;
nlon=lon;
nlat=lat;
i=length(lon);
while i>1
    temp=i-1;
    j=temp;
    while ~isnan(lon(j))
        j=j-1;
    end
    nlon(i)=temp-j;
    nlat(i)=0;
    i=j;
end


% i=length(lat);
% while i>1
%     temp=i-1;
%     j=temp;
%     while ~isnan(lat(j))
%         j=j-1;
%     end
%     nlat(i)=0;
%     i=j;
% end



a=[flipud(nlon) flipud(nlat)];



a=a(1:length(a)-1,:);

fid=fopen('sb_coast.bln','a');
for i=1:length(a)
    if a(i,2)==0
       fprintf(fid,'%d %d\n',a(i,:));
    else
        fprintf(fid,'%3.6f %3.6f\n',a(i,:));
    end
end
fclose(fid)

mofangbao

不错，surfer的底图来源方法也越来越多了

夏天爱裤衩

兄台试过这个程序吗，我试了一下，好像有不少错误

poet411

一直在找这个转化方法，不知道能行吗。多谢。。。。。

poet411

程序不能用啊。