Cs2csFromMatlab

From Marine Science
Jump to: navigation, search

Main Page > How To Do Stuff > Coordinate Transforms > Cs2csFromMatlab

From Windows, using FWTools:

% Hamish & Cerys, 9 May 2008
%  This script is public domain code
% setup environment for FWTools
FWTOOLS_DIR='C:\PROGRA~1\FWTOOL~1.0';
setenv('FWTOOLS_DIR', FWTOOLS_DIR);
setenv('PATH', [FWTOOLS_DIR '\bin;' FWTOOLS_DIR '\python;' getenv('PATH')]);
setenv('PYTHONPATH', [FWTOOLS_DIR '\pymod']);
setenv('PROJ_LIB', [FWTOOLS_DIR '\proj_lib']);
setenv('GEOTIFF_CSV', [FWTOOLS_DIR '\data']);
setenv('GDAL_DATA', [FWTOOLS_DIR '\data']);
setenv('GDAL_DRIVER_PATH', [FWTOOLS_DIR '\gdal_plugins']);
% save two column array to a text file:
save coords_LL.txt coords -ASCII
!cs2cs -f "%.3f" +init=epsg:4326 +to +init=epsg:2131 < coords_LL.txt > coord_NTai2k.txt
% load in result  format is eastings,northings,elevation  (x,y,z)
load 'coord_NTai2k.txt' 
% remove unneeded cruft
%delete('coord_LL.txt')
%delete('coord_NTai2k.txt')


If path stuff is already taken care of (e.g. Linux):


% coordinate conversions using cs2cs     Hamish Bowman 19 Oct 2006
%  This script is public domain code
%   http://proj.maptools.org
%
%  New Zealand EPSG Codes
%
% Projection            Datum   EPSG Code
% NZ Map Grid           NZGD49  27200
% NZ Trans. Merc.       NZGD2k  2193
% N. Taieri Circuit     NZGD2k  2131
% Lat/Lon               WGS84   4326
%
% If converting NZGD_1949 <-> NZGD_2000 use the LINZ NTv2 ("nad") distortion grid
%  to perform the datum transform instead of the default 7 term scaling parameter
%   default 7 term accuracy ~ 1-4m
%   accuracy using distortion grid ~ 10-30cm
%   accuracy using 3 term scaling ~ 5m+
% get the grid here:
%   http://bambi.otago.ac.nz/public/general/mapping/nzgd2kgrid0005.gsb
%
% to use the nzgd49 distortion grid to go from lat/lon WGS84 to NZ Map Grid:
%  cs2cs -f "%.3f" +init=epsg:4326 +to +init=epsg:27200 +nadgrids=nzgd2kgrid0005.gsb
%  170d43'E 45d46'S
%  2332265.548     5490967.665 -0.000
% notice how the third (z) result is now ~ 0.0 suggesting low error in the other two axes 


% start with 170d43'E  45d46'S  (170.716666667  -45.766666667)
fix.degE = 170;
fix.minE = 43.5000;
fix.degS = 45;
fix.minS = 46.0000
% I'm not sure how to pass data via stdin and stdout to external programs in
%  Matlab, so write to a file for external communications
fd1 = fopen('coord_LL.txt', 'wt');
% format is LONGITUDE LATITUDE (x,y  eastings,northings)
%   use 170d45.1234'E or 170.716666667 format, whatever's easier
fprintf(fd1, '%dd%.9fE %dd%.9fS\n', ...
  fix.degE, fix.minE, fix.degS,fix.minS );
fclose(fd1);
% convert lat/lon WGS84 to North Taieri Circuit 2000
!cs2cs -f "%.3f" +init=epsg:4326 +to +init=epsg:2131 < coord_LL.txt > coord_NTai2k.txt
% load in result  format is eastings,northings,elevation  (x,y,z)
load 'coord_NTai2k.txt'