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