Netopology wkt for coordinates

1433
4
12-09-2020 08:13 AM
AndréGrce
New Contributor
I am trying to read coordinates from a shapefile. But I get wrong coordinates.
I am not sure if I use the correct wkt strings. Is there any function i can use to get the wkt string for a given shapefile?
 
I use these librarys:
 
using NetTopologySuite.IO;
using ProjNet.CoordinateSystems;
using ProjNet.CoordinateSystems.Transformations;
 
And this is my code:
 
CoordinateSystemFactory coordinateSystemFactory = new CoordinateSystemFactory();
string UTMWGS84_WKT2 = "PROJCS[\"WGS 84 / UTM zone 32N\", 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[\"Transverse_Mercator\"], PARAMETER[\"latitude_of_origin\", 0], PARAMETER[\"central_meridian\", 9], PARAMETER[\"scale_factor\", 0.9996], PARAMETER[\"false_easting\", 500000], PARAMETER[\"false_northing\", 0], AUTHORITY[\"EPSG\", \"32632\"], AXIS[\"Easting\", EAST], AXIS[\"Northing\", NORTH]]";
CoordinateSystem source = coordinateSystemFactory.CreateFromWkt(@UTMWGS84_WKT2);
string UTMWGS84_WKT = "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[\"Transverse_Mercator\"]," + "PARAMETER[\"latitude_of_origin\",0]," + "PARAMETER[\"central_meridian\",9]," + "PARAMETER[\"scale_factor\",0.9996]," + "PARAMETER[\"false_easting\",500000]," + "PARAMETER[\"false_northing\",0]," + "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]"; 
CoordinateSystem target = coordinateSystemFactory.CreateFromWkt(UTMWGS84_WKT);
CoordinateSystem target2 = coordinateSystemFactory.CreateFromWkt(UTMWGS84_WKT);
CoordinateTransformationFactory trf = new CoordinateTransformationFactory();
CoordinateTransformation tr = (CoordinateTransformation)trf.CreateFromCoordinateSystems(source, target);
var filename = @"C:\\myshapefile.shp";
 
using (var reader = new ShapefileDataReader(filename, NetTopologySuite.Geometries.GeometryFactory.Default))
{
for (var i = 0; i < dbaseHeader.NumRecords; i++)
        {
        reader.Read();
                string id = reader.GetString(1);
 
                NetTopologySuite.Geometries.Geometry geometry = reader.Geometry;
                NetTopologySuite.Geometries.Coordinate coordinate = new NetTopologySuite.Geometries.Coordinate();
                NetTopologySuite.IO.WKTWriter wktWrirer = new NetTopologySuite.IO.WKTWriter();
 
                coordinate = geometry.Coordinate;
 
                Console.WriteLine("Id " + i + ": " + id);
                Console.WriteLine("X Mercator " + i + ": " + coordinate.X);
                Console.WriteLine("Y Mercator " + i + ": " + coordinate.Y);
 
                double[] coordinateMercator = new double[] { coordinate.X, coordinate.Y };
                double[] coordinateWGS84 = tr.MathTransform.Transform(coordinateMercator);
 
                Console.WriteLine("X WGS84 " + i + ": " + coordinateWGS84[0]);
                Console.WriteLine("Y WGS84 " + i + ": " + coordinateWGS84[1]);
}
}
0 Kudos
4 Replies
JoshuaBixby
MVP Esteemed Contributor
0 Kudos
JoshuaBixby
MVP Esteemed Contributor

Given you are on an Esri online community, I will assume you have access to ArcGIS Desktop software.  You can use ArcPy to extract WKT from shapes in a shape file.

import arcpy

shp_file = # path to shape file

with arcpy.da.SearchCursor(shp_file, "SHAPE@WKT") as cur:
    for row in cur:
        print(row[0])
0 Kudos
AndréGrce
New Contributor

Thank you Joshua, 

I can get wkt strings for shape geometries, but I need wkt for the coordinate system used. Otherwise I can not transform the x/y to lat/long coodinates. Preferrably a function from a C# class library.

0 Kudos
JoshuaBixby
MVP Esteemed Contributor

There are a couple different ways to get the WKID of the shapes.  One way is through the spatialReference object of the dataset:

import arcpy
>>>
shp_file = # path to shape file
arcpy.Describe(shp_file).spatialReference.factoryCode
0 Kudos