[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Forward projection geodetic to X-Y

From: Lester Anderson
Subject: Forward projection geodetic to X-Y
Date: Mon, 1 Mar 2021 08:23:20 +0000


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');

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

Original Matlab code, gdalread (above) is nearest equivalent:

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.


reply via email to

[Prev in Thread] Current Thread [Next in Thread]