Friday 1 May 2015

.net - How to perform transformations using ProjNet?



While writing some unit tests I thought it would be convenient/possible to use ProjNet to convert from EPSG:4326 into several other coordinate systems (3875,3395). I figured it would be as easy as grabbing the WKT values from spatialreference.org:



Add a reference to ProjNet and running the following code:


var cf = new ProjNet.CoordinateSystems.CoordinateSystemFactory ( );
var f = new ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory ( );
var sys4326 = cf.CreateFromWkt ( wkt4326 );
var sys3857 = cf.CreateFromWkt ( wkt3857 );
var sys3395 = cf.CreateFromWkt ( wkt3395 );
var transformTo3875 = f.CreateFromCoordinateSystems ( sys4326, sys3857 );
var transformTo3395 = f.CreateFromCoordinateSystems ( sys4326, sys3395 );


But 3857 and 3395 could not be parsed. I googled around and found a workable 3857 though I haven't confirmed it produces a correct result. I continue to get an error parsing 3395 (Expecting ('PROJECTION') but got a 'UNIT' at line 1 column 296.).


Perhaps the "OGC WKT" can be modified in a predicable way to work with ProjNet? Or perhaps ProjNet is not the best .NET solution for re-projecting coordinates?


const string wkt4326 = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]";
const string wkt3857 = "PROJCS[\"Popular Visualisation CRS / Mercator\", GEOGCS[\"Popular Visualisation CRS\", DATUM[\"Popular Visualisation Datum\", SPHEROID[\"Popular Visualisation Sphere\", 6378137, 0, AUTHORITY[\"EPSG\",\"7059\"]], TOWGS84[0, 0, 0, 0, 0, 0, 0], AUTHORITY[\"EPSG\",\"6055\"] ], PRIMEM[\"Greenwich\", 0, AUTHORITY[\"EPSG\", \"8901\"]], UNIT[\"degree\", 0.0174532925199433, AUTHORITY[\"EPSG\", \"9102\"]], AXIS[\"E\", EAST], AXIS[\"N\", NORTH], AUTHORITY[\"EPSG\",\"4055\"] ], PROJECTION[\"Mercator\"], PARAMETER[\"False_Easting\", 0], PARAMETER[\"False_Northing\", 0], PARAMETER[\"Central_Meridian\", 0], PARAMETER[\"Latitude_of_origin\", 0], UNIT[\"metre\", 1, AUTHORITY[\"EPSG\", \"9001\"]], AXIS[\"East\", EAST], AXIS[\"North\", NORTH], AUTHORITY[\"EPSG\",\"3785\"]]";
const string wkt3395 = "PROJCS[\"WGS 84 / World Mercator\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],PROJECTION[\"Mercator_1SP\"],PARAMETER[\"central_meridian\",0],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],AUTHORITY[\"EPSG\",\"3395\"],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]";

Answer



I found that manually conforming my WKT to the BNF described here eliminates the errors. Content repeated below:


 =  |  | 
= PROJCS["", , , {,}* ]

= PROJECTION[""]
= PARAMETER["", ]
=

= GEOGCS["", , , ]
= DATUM["", ]
= SPHEROID["", , ]
= NOTE: semi-major axis is measured in meters and must be > 0.
=
= PRIMEM["", ]

=

=
=
= UNIT["", ]
=

= GEOCCS["", , , ]

So for my specific problem (EPSG:3395) the resulting WKT would be:



PROJCS["WGS 84 / World Mercator",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],

PROJECTION["Mercator_1SP"],
PARAMETER["latitude_of_origin",0], // required, missing from spatialreference.org
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","3395"]]

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...