Forward projection geodetic to X-Y

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Forward projection geodetic to X-Y

Lester Anderson
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



Reply | Threaded
Open this post in threaded view
|

Re: Forward projection geodetic to X-Y

PhilipNienhuis
(mapping pkg maintainer here)

I'd love to have projection functions in the mapping package but I haven't
been able to find the time to investigate thoroughly; nor has any volunteer
emerged yet. So it's a slowly evolving project.

There is an OF package octproj that actually has all required building
blocks. The future mapping pkg projfwd and projinv functions would be just
wrapper functions for the functions op_fwd and op_inv in octproj.

Against that background I'm also looking for a suitable WKT string parser
with a suitable license or write it myself (I have some skeleton lying
around), but op_fwd and op_inv (and op_transform) seem to be able to process
them.

HTH

Philip



--
Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html