Forward projection geodetic to X-Y
Lester Anderson <[email protected]>
| Newsgroups | gmane.comp.gnu.octave.general |
|---|---|
| Message-ID | <CAE3taFDrjS_6F_Ou9=7YjnJvo1D1+1fPf-HhKRj2rufu7hRSKw@mail.gmail.com> |
Hello,
I am converting an existing Matlab code to work under Octave, but cannot
see a direct equivalent to the projfwd function of Matlab.
I am using a subset of packages to deal with inputs and the mapping:
pkg load netcdf
pkg load mapping
pkg load io
[RSTAT, RINFO,
BANDS]=gdalread('D:/Matlab_files/Bai_crustal_model/tifForProjection.tif');
proj=getfield(RINFO,"Projection");
When I look at the value of proj, it is a character string and not a proj5
format:
PROJ = PROJCS["World_Mercator",GEOGCS["WGS
84",DATUM["WGS_1984",SPHEROID["WGS
84",6378137,298.257223563,AUTHORITY["EPSG","703
0"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["
Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_nort
hing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]
Original Matlab code, gdalread (above) is nearest equivalent:
proj=geotiffinfo('Data\tifForProjection.tif');
The code requires that input Longitude and Latitude are converted to the
projected units (metres) of the defined projection (Mercator):
% create longitude and latitude grid
GLong = [minLongitude:(2/60):maxLongitude];
GLat = [minLatitude:(2/60):maxLatitude];
%[GridLong,GridLat] = meshgrid(GLong,GLat);
[GridLong,GridLat] = meshgrid(GLong,GLat);
%[GridLong,GridLat] = projfwd(proj,GridLat,GridLong); - convert X,Y to
projected units
[GridLong,GridLat] = projfwd(GridLong,GridLat,proj);
My query is what the Octave equivalent would be for projfwd (mapping
package?) and what proj in this case should be represented as, since I do
not have access to Matlab. Should proj be the proj4 structure?
The Mercator projection (proj4) would be:
+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
Further on in the code the inverse projection is required, going from
projected units to Longitude and Latitude; could not see a simple option in
Mapping that used the projection information.
Any pointers would be useful.
Thanks
Lester