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
lmpx.com only provides a reader for public news (NNTP) servers. It is not affiliated with the servers or forums shown here and is not responsible for the content of articles, which is written by their respective authors.