Example of application
We need to tell convertCoordinates what the input coordinate system is, and what the output coordinate system is (CS1 and CS2). We can do so by entering the name, type and /or EPSG reference code using keyword value pairs. An example below:
x = 165e3; y = 420e3; [lon,lat] = convertCoordinates(x,y,'CS1.name','Amersfoort / RD New',... 'CS2.code',4326)
lon =
5.5321
lat =
51.7686
CS1.name refers to the name of Coordinate System 1. Coordinate System 2 is identified by a code: 4326. Alternatively, also the name (WGS 84) of the second coordinate system can be entered:
try [lon,lat] = convertCoordinates(x,y,'CS1.name','Amersfoort / RD New',... 'CS2.name','WGS 84'); catch err = lasterror; disp(err.message) end
Error using ==> ConvertCoordinatesFindCoordRefSys at 81 3 coordinate reference systems can be found with options: code:'', name:'WGS 84', type:'' code: '4326', name: 'WGS 84', type: 'geographic 2D' code: '4978', name: 'WGS 84', type: 'geocentric' code: '4979', name: 'WGS 84', type: 'geographic 3D'
Because the name 'WGS 84' does not refer to a unique coordinate system, an error message is returned. From the message we can see that there are three systems known by that name. we want to convert coordinates to geographic 2D. If we also specify the type, this leads to a unique coordinate system.
[lon,lat] = convertCoordinates(x,y,'CS1.name','Amersfoort / RD New',... 'CS2.name','WGS 84','CS2.type','Geographic 2D')
lon =
5.5321
lat =
51.7686
Finding the name or code of a coordinate system
If convertCoordinates cannot find an exact match for the coordinate system name, it returns an error message with approximate matches. This is useful if you do not know the exact code or name of a coordinate system. Make sure to enter a name that does not match any entry, as for instance 'amersfoo'.
try [lon,lat] = convertCoordinates(x,y,'CS1.name','amersfoo',... 'CS2.code',4326) catch err = lasterror; disp(err.message) end
WARNING: no exact match of coordinate system name is found, aproximate match tried Error using ==> ConvertCoordinatesFindCoordRefSys at 81 5 coordinate reference systems can be found with options: code:'', name:'amersfoo', type:'' code: '4289', name: 'Amersfoort', type: 'geographic 2D' code: '7415', name: 'Amersfoort / RD New + NAP height', type: 'compound' code: '28991', name: 'Amersfoort / RD Old', type: 'projected' code: '28992', name: 'Amersfoort / RD New', type: 'projected' code: '62896405', name: 'Amersfoort (deg)', type: 'geographic 2D'
Apparently there are five coordinate systems who's name contains 'amersfoo'. Amersfoort / RD New has the code 28992.
[lon,lat] = convertCoordinates(x,y,'CS1.code',28992,'CS2.code',4326)
lon =
5.5321
lat =
51.7686
Speeding up the routine in batch mode
The slowest part of the routine is reading the EPSG database. If the routine is called several times in a batch mode, it is faste to load the data only once, and the pass it to the function. Compare the following options.
% without preloading EPSG tic for ii = 1:10 [lon,lat] = convertCoordinates(x,y,'CS1.code',28992,'CS2.code',4326); end toc
Elapsed time is 3.795889 seconds.
% with preloading EPSG tic EPSG = load('EPSG'); for ii = 1:10 [lon,lat] = convertCoordinates(x,y,EPSG,'CS1.code',28992,'CS2.code',4326); end toc
Elapsed time is 0.421173 seconds.
Evaluating the output
To evaluate the output, call the funtion with OPT:
[lon,lat,OPT] = convertCoordinates(x,y,EPSG,'CS1.code',28992,'CS2.code',4326); disp(OPT);
CS1: [1x1 struct]
proj_conv1: [1x1 struct]
datum_trans: [1x1 struct]
proj_conv2: 'no projection conversion necessary'
CS2: [1x1 struct]
From the OPT structure, all relevant parameters and conversion options can be retrieved. This can be done with the Array editor, or (easier), by using the OET function var2evalstr.
var2evalstr(OPT)
ans =
OPT.CS1.name = 'Amersfoort / RD New';
OPT.CS1.code = 28992;
OPT.CS1.type = 'projected';
OPT.CS1.geoRefSys.name = 'Amersfoort';
OPT.CS1.geoRefSys.code = 4289;
OPT.CS1.coordSys.name = 'Cartesian 2D CS. Axes: easting; northing (X;Y). Orientations: east; north. UoM: m.';
OPT.CS1.coordSys.code = 4499;
OPT.CS1.ellips.code = 7004;
OPT.CS1.ellips.name = 'Bessel 1841';
OPT.CS1.ellips.inv_flattening = 299.1528;
OPT.CS1.ellips.semi_major_axis = 6377397.155;
OPT.CS1.ellips.semi_minor_axis = 6356078.9628;
OPT.CS1.UoM.name = 'metre';
OPT.CS1.UoM.code = 9001;
OPT.CS1.datum.code = 6289;
OPT.CS1.datum.name = 'Amersfoort';
OPT.proj_conv1.code = 19914;
OPT.proj_conv1.name = {'RD New'};
OPT.proj_conv1.method.code = 9809;
OPT.proj_conv1.method.name = 'Oblique Stereographic';
OPT.proj_conv1.param.codes = [8801 8802 8805 8806 8807];
OPT.proj_conv1.param.name = {'Latitude of natural origin' 'Longitude of natural origin' 'Scale factor at natural origin' 'False easting' 'False northing'};
OPT.proj_conv1.param.value = [52.0922 5.23155 0.999908 155000 463000];
OPT.proj_conv1.param.UoM.codes = [9110 9110 9201 9001 9001];
OPT.proj_conv1.param.UoM.name = {'sexagesimal DMS' 'sexagesimal DMS' 'unity' 'metre' 'metre'};
OPT.datum_trans.code = 15934;
OPT.datum_trans.name = {'Amersfoort to WGS 84 (3)'};
OPT.datum_trans.direction = 'normal';
OPT.datum_trans.alt_code = [1112 1672 15934];
OPT.datum_trans.alt_name = {'Amersfoort to WGS 84 (1)' 'Amersfoort to WGS 84 (2)' 'Amersfoort to WGS 84 (3)'};
OPT.datum_trans.alt_direction = {'normal' 'normal' 'normal'};
OPT.datum_trans.alt_deprecated = {'FALSE' 'FALSE' 'FALSE'};
OPT.datum_trans.params.value = [565.237 50.0087 465.658 1.9725 -1.7004 9.0677 4.0812];
OPT.datum_trans.params.codes = 8605:8611;
OPT.datum_trans.params.ind = 5:11;
OPT.datum_trans.params.name = {'X-axis translation' 'Y-axis translation' 'Z-axis translation' 'X-axis rotation' 'Y-axis rotation' 'Z-axis rotation' 'Scale difference'};
OPT.datum_trans.params.UoM.codes = [9001 9001 9001 9109 9109 9109 9202];
OPT.datum_trans.params.UoM.ind = [1 1 1 51 51 51 66];
OPT.datum_trans.params.UoM.sourceN = {'metre' 'metre' 'metre' 'microradian' 'microradian' 'microradian' 'parts per million'};
OPT.datum_trans.params.UoM.sourceT = {'length' 'length' 'length' 'angle' 'angle' 'angle' 'scale'};
OPT.datum_trans.params.UoM.fact_b = ones(1,7);
OPT.datum_trans.params.UoM.fact_c = [1 1 1 1e+006 1e+006 1e+006 1e+006];
OPT.datum_trans.params.UoM.targetC = [9001 9001 9001 9101 9101 9101 9201];
OPT.datum_trans.params.UoM.targetI = [1 1 1 43 43 43 65];
OPT.datum_trans.params.UoM.targetN = {'metre' 'metre' 'metre' 'radian' 'radian' 'radian' 'unity'};
OPT.datum_trans.method_code = 9607;
OPT.datum_trans.method_name = 'Coordinate Frame rotation';
OPT.datum_trans.ellips1 = 'CS1';
OPT.datum_trans.ellips2 = 'CS2';
OPT.proj_conv2 = 'no projection conversion necessary';
OPT.CS2.name = 'WGS 84';
OPT.CS2.code = 4326;
OPT.CS2.type = 'geographic 2D';
OPT.CS2.geoRefSys.name = 'WGS 84';
OPT.CS2.geoRefSys.code = 4326;
OPT.CS2.coordSys.name = 'Ellipsoidal 2D CS. Axes: latitude; longitude. Orientations: north; east. UoM: degree';
OPT.CS2.coordSys.code = 6422;
OPT.CS2.ellips.code = 7030;
OPT.CS2.ellips.name = 'WGS 84';
OPT.CS2.ellips.inv_flattening = 298.2572;
OPT.CS2.ellips.semi_major_axis = 6378137;
OPT.CS2.ellips.semi_minor_axis = 6356752.3142;
OPT.CS2.UoM.name = 'degree (supplier to define representation)';
OPT.CS2.UoM.code = 9122;
OPT.CS2.datum.code = 6326;
OPT.CS2.datum.name = 'World Geodetic System 1984';
The contents of the OPT structure depends on the conversions and transformations performed. Consider this example (of a very unlikely coordinate conversion)
[x,y,OPT] = convertCoordinates(1,1,EPSG,... 'CS1.name','IRENET95 / Irish Transverse Mercator',... 'CS2.name','NAD27 / Alberta 3TM ref merid 111 W'); disp(OPT)
CS1: [1x1 struct]
proj_conv1: [1x1 struct]
datum_trans: 'no direct transformation available'
datum_trans_to_WGS84: [1x1 struct]
WGS84: [1x1 struct]
datum_trans_from_WGS84: [1x1 struct]
proj_conv2: [1x1 struct]
CS2: [1x1 struct]
- CS1 identifies the input coordinate system
- Proj_conv1 contains the transformation parameters from Cartesian coordinate system 1 to a geodetic (lat/lon) system
- Apparently no direct transformation is available from the GRS 1980 ellipsoid of coordinate system 1 ('OPT.CS1.ellips.name') to that of CS2 (Clarke 1866). Therefore the coordinates are first transformed from the GRS 1980 to the WGS 84 ellipsoid, and only then to Clarke 1866
- Finally the coordinates are converted to the NAD27 coordinate reference system
Multiple conversion options
For some coordinate conversions, more than one routine is available. And example is the RD to WGS 84 conversion. We can evaluate the available methods from the OPT structure.
[lon1,lat1,OPT] = convertCoordinates(1e5,5e5,EPSG,'CS1.code',28992,'CS2.code',4326); var2evalstr(OPT.datum_trans)
ans =
variable1.code = 15934;
variable1.name = {'Amersfoort to WGS 84 (3)'};
variable1.direction = 'normal';
variable1.alt_code = [1112 1672 15934];
variable1.alt_name = {'Amersfoort to WGS 84 (1)' 'Amersfoort to WGS 84 (2)' 'Amersfoort to WGS 84 (3)'};
variable1.alt_direction = {'normal' 'normal' 'normal'};
variable1.alt_deprecated = {'FALSE' 'FALSE' 'FALSE'};
variable1.params.value = [565.237 50.0087 465.658 1.9725 -1.7004 9.0677 4.0812];
variable1.params.codes = 8605:8611;
variable1.params.ind = 5:11;
variable1.params.name = {'X-axis translation' 'Y-axis translation' 'Z-axis translation' 'X-axis rotation' 'Y-axis rotation' 'Z-axis rotation' 'Scale difference'};
variable1.params.UoM.codes = [9001 9001 9001 9109 9109 9109 9202];
variable1.params.UoM.ind = [1 1 1 51 51 51 66];
variable1.params.UoM.sourceN = {'metre' 'metre' 'metre' 'microradian' 'microradian' 'microradian' 'parts per million'};
variable1.params.UoM.sourceT = {'length' 'length' 'length' 'angle' 'angle' 'angle' 'scale'};
variable1.params.UoM.fact_b = ones(1,7);
variable1.params.UoM.fact_c = [1 1 1 1e+006 1e+006 1e+006 1e+006];
variable1.params.UoM.targetC = [9001 9001 9001 9101 9101 9101 9201];
variable1.params.UoM.targetI = [1 1 1 43 43 43 65];
variable1.params.UoM.targetN = {'metre' 'metre' 'metre' 'radian' 'radian' 'radian' 'unity'};
variable1.method_code = 9607;
variable1.method_name = 'Coordinate Frame rotation';
variable1.ellips1 = 'CS1';
variable1.ellips2 = 'CS2';
By default convertCoordinates chooses the newest non-deprecated method. (It is assumed that the method with the highest code is the newest). This can be overridden by defining a datum_trans.code:
[lon2,lat2,OPT] = convertCoordinates(1e5,5e5,EPSG,'CS1.code',28992,'CS2.code',4326,... 'datum_trans.code',1672); lat1-lat2
ans = 8.4067e-008