坐标转换和插值之matlab程序

 

Matlab版本：8.1.0.604 (R2013a)

这里需要用到的命令有：maps, defaultm, projfwd, projinv

Matlab版本：8.1.0.604 (R2013a)

1，  查看有哪些地图投影方式：

运行maps查看有哪些投影方式

>> maps

 

MapTools Projections

 

CLASS                 NAME                                ID STRING       

Cylindrical           Balthasart Cylindrical              balthsrt        

Cylindrical           Behrmann Cylindrical                behrmann        

Cylindrical           Bolshoi Sovietskii Atlas Mira*      bsam            

Cylindrical           Braun Perspective Cylindrical*      braun           

Cylindrical           Cassini Cylindrical-Standard        cassinistd      

Cylindrical           Cassini Cylindrical                 cassini         

Cylindrical           Central Cylindrical*                ccylin          

Cylindrical           Equal Area Cylindrical              eqacylin        

Cylindrical           Equidistant Cylindrical             eqdcylin        

Cylindrical           Gall Isographic                     giso            

Cylindrical           Gall Orthographic                   gortho          

Cylindrical           Gall Stereographic*                 gstereo         

Cylindrical           Lambert Cylindrical                 lambcyln        

Cylindrical           Mercator Cylindrical                mercator        

Cylindrical           Miller Cylindrical*                 miller          

Cylindrical           Plate Carree                        pcarree         

Cylindrical           Transverse Mercator                 tranmerc        

Cylindrical           Trystan Edwards Cylindrical         trystan         

Cylindrical           Universal Transverse Mercator (UTM)    utm             

Cylindrical           Wetch Cylindrical*                  wetch           

Pseudocylindrical     Apianus II*                         apianus         

Pseudocylindrical     Collignon                           collig          

Pseudocylindrical     Craster Parabolic                   craster         

Pseudocylindrical     Eckert I*                           eckert1         

Pseudocylindrical     Eckert II                           eckert2         

Pseudocylindrical     Eckert III*                         eckert3         

Pseudocylindrical     Eckert IV                           eckert4         

Pseudocylindrical     Eckert V*                           eckert5         

Pseudocylindrical     Eckert VI                           eckert6          

Pseudocylindrical     Flat-Polar Parabolic                flatplrp        

Pseudocylindrical     Flat-Polar Quartic                  flatplrq        

Pseudocylindrical     Flat-Polar Sinusoidal               flatplrs        

Pseudocylindrical     Fournier                            fournier        

Pseudocylindrical     Goode Homolosine                    goode           

Pseudocylindrical     Hatano Assymmetrical Equal Area     hatano          

Pseudocylindrical     Kavraisky V                         kavrsky5        

Pseudocylindrical     Kavraisky VI                        kavrsky6        

Pseudocylindrical     Loximuthal*                         loximuth        

Pseudocylindrical     Modified Sinusoidal (Tissot)*       modsine          

Pseudocylindrical     Mollweide                           mollweid        

Pseudocylindrical     Putnins P5*                         putnins5        

Pseudocylindrical     Quartic Authalic                    quartic         

Pseudocylindrical     Robinson*                           robinson        

Pseudocylindrical     Sinusoidal                          sinusoid        

Pseudocylindrical     Wagner IV                           wagner4         

Pseudocylindrical     Winkel I*                           winkel          

Conic                 Equal Area Conic (Albers)-Standard    eqaconicstd     

Conic                 Equal Area Conic (Albers)           eqaconic        

Conic                 Equidistant Conic-Standard          eqdconicstd     

Conic                 Equidistant Conic                   eqdconic        

Conic                 Lambert Conformal Conic-Standard    lambertstd      

Conic                 Lambert Conformal Conic             lambert         

Conic                 Murdoch I Conic*                    murdoch1        

Conic                 Murdoch III Minimum Error Conic*    murdoch3        

PolyConic             Polyconic-Standard                  polyconstd      

PolyConic             Polyconic                           polycon         

PolyConic             Van Der Grinten I*                  vgrint1         

PseudoConic           Bonne                               bonne           

PseudoConic           Werner                              werner          

Azimuthal             Breusing Harmonic Mean*             breusing        

Azimuthal             Equal Area Azimuthal (Lambert)      eqaazim         

Azimuthal             Equidistant Azimuthal*              eqdazim         

Azimuthal             Globe                               globe           

Azimuthal             Gnomonic*                           gnomonic        

Azimuthal             Orthographic*                       ortho           

Azimuthal             Stereographic                       stereo          

Azimuthal             Universal Polar Stereographic       ups             

Azimuthal             Vertical Perspective*               vperspec        

Pseudoazimuthal       Wiechel Equal Area*                 wiechel         

Modified Azimuthal    Aitoff*                             aitoff          

Modified Azimuthal    Briesemeister*                      bries           

Modified Azimuthal    Hammer*                             hammer          

 

* Denotes availability for sphere only

 

2，  定义投影结构

这里对mercator投影结构进行定义,然后对mstruct进行补全

>>mstruct=defaultm('mercator');

%定义地球坐标系WGS84

>>mstruct.geoid=almanac('earth','wgs84','meters');

>>mstruct=defaultm(mstruct);

 

3，  进行投影和数据转换

地理坐标投影到直角坐标:

>>[x,y] =projfwd(mstruct,lat,lon);

直角坐标投影到地理坐标:

>>[lat,lon]=projinv(mstruct,x,y);

  

 

4，例子：读取netcdf地形文件，并转换到Mercator投影的直接坐标系中，由于转换后坐标间隔变成非均匀的，不利于某些程序处理，这里进行了均匀插值。

首先准备好地形网格文件，该文件是利用GMT软件的grdcut从gebco_08中截取的,该文件信息如下。

>> ncdisp('data\topo.nc');

Source:

           E:\FedoraShare\work\MohoInversion\matlabCode\bouguerCorrection\data\topo.nc

Format:

           netcdf4

Global Attributes:

           Conventions = 'COARDS, CF-1.5'

           title       = 'Produced by grdcut'

           history     = 'grdcut -R112/122/25/35 /mnt/hgfs/FedoraShare/data/topoData/gebco_08/gebco_08.nc -Gtopo.nc'

           GMT_version = '5.1.1_r12476 [64-bit]'

           node_offset = 1

Dimensions:

           x = 1200

           y = 1200

Variables:

    x

           Size:       1200x1

           Dimensions: x

           Datatype:   double

           Attributes:

                       long_name    = 'user_x_unit'

                       actual_range = [1.12e+02 1.22e+02]

    y

           Size:       1200x1

           Dimensions: y

           Datatype:   double

           Attributes:

                       long_name    = 'user_y_unit'

                       actual_range = [2.50e+01 3.50e+01]

    z

           Size:       1200x1200

           Dimensions: x,y

           Datatype:   single

           Attributes:

                       long_name    = 'user_z_unit'

                       _FillValue   = NaN

                       actual_range = [-7.37e+02  2.06e+03]

 

主程序：

clear;

%% add paths for functions

addpath BasicFunctions

clear;

%定义读取经纬度范围(degree)

latMin=25;

latMax=35;

lonMin=112;

lonMax=122;

%%读取地形网格文件

%降采样率

samp=1;%不降采样

%samp=2;%每两个点采样一个点

%文件位置

topoGridFile = 'Data\topo.nc';

%读取变量

[topoGrid,lon,lat] = getcdf(topoGridFile, samp,[ lonMin lonMax latMin latMax ]);

 

%%坐标系变换并均匀插值

%定义变换方式

mstruct=defaultm('mercator');

%定义地球坐标系

mstruct.geoid=almanac('earth','wgs84','meters');

mstruct=defaultm(mstruct);

%坐标系变换

[lon,lat]=meshgrid(lon,lat);

[x,y] =projfwd(mstruct,lat,lon);

%定义直角坐标范围(m)和插值点间隔

xMin=x(1,1);

xMax=x(1,end);

yMin=y(1,1);

yMax=y(end,1);

dx=1000;%m

dy=1000;%m

%dx=5000;%m

%dy=5000;%m

 

%重新插值成均匀网格

X=xMin:dx:xMax;

Y=yMin:dy:yMax;

[X,Y]=meshgrid(X,Y);

topoGrid_merc=interp2(x,y,topoGrid,X,Y,'linear');

% topoGrid_merc=interp2(x,y,topoGrid,X,Y,'spline');

 

X=X-xMin;

Y=Y-yMin;

subplot(1,2,1);

ll=[-500,-100,0,50,100,150,500,1000,1500,2000];

contour(X,Y,topoGrid_merc,ll);

caxis([-500,2000]);

colorbar('location','SouthOutside');

subplot(1,2,2);

contour(lon,lat,topoGrid,ll);

colorbar('location','SouthOutside');

 

其中用到的BasicFunctions中的子程序getcdf：

function [t,x,y] = getcdf(file,Sampling,Dims);

% [t,x,y] = getcdf(file,Sampling,Dims);

%   Sampling is an integer where you can specify a desampling factor to get

%   a low res version of the data

%   Dims is an [Xmin Xmax Ymin Ymax] windowing option

 

  

    if nargin<2

        Sampling=1;

    end

   

    % Open the rotated grid

    ncid=netcdf.open(file,'NOWRITE');

    x = netcdf.getVar(ncid,0,'double');

    y = netcdf.getVar(ncid,1,'double');

   

    if nargin==3

        IndX = find(x>Dims(1)&x

        IndY = find(y>Dims(3)&y

        x = x(IndX);

        y = y(IndY);

        t = double(netcdf.getVar(ncid,2,[IndX(1)-1 IndY(1)-1],[length(IndX) length(IndY)],[Sampling Sampling]))';

    else

        x=x(1:Sampling:end);

        y=y(1:Sampling:end);

        t = netcdf.getVar(ncid,2,[0 0],[length(x) length(y)],[Sampling Sampling],'double')';

    end

   

    netcdf.close(ncid);

 

   

end

 

 

转换和插值结果：

http://s4/mw690/004f6djqgy6KOWPQ9Pl63&690图一：左图是坐标系转换和线性插值到1km后的结果，右边是原坐标系未插值的结果。

 

 

http://s14/mw690/004f6djqgy6KOWU29el4d&690图二：左图是坐标系转换和线性插值到5km后的结果，右边是原坐标系未插值的结果。

 

 

http://s16/mw690/004f6djqgy6KOWVjAar4f&690

 

图三：左图是坐标系转换和样条插值到5km后的结果，右边是原坐标系未插值的结果。

 

 

该坐标转换程序为自己摸索，本例使用的是正轴墨卡托投影和WGS84坐标系，并非北京54坐标系/西安80坐标系/UTM投影，其他坐标系和投影方式还须进一步探索