How to get coordinate system of shapefile?

[kjeremy]

I am reading in a shape file using to test my code:

           var reader = new ShapeDataReader(path);
            var bounds = reader.ShapefileBounds;
            foreach (var thing in reader.ReadByMBRFilter(bounds))
            {
                var attrs = thing.Attributes;
                var geo = thing.Geometry;
                if (geo != null)
                {
                    var coords = geo.Coordinates;
                }
            }

path is the path to the .shp file. Is this correct? Or should I be using ShapefileDataReader? I've seen examples online of both.

One thing is that I can't figure out how to get the coordinate system via code. I need to project it to WGS84. The srid on the geometry is 0 so I'm not sure how to get the information I need. Do I need to do something special to read the .prj file to get the information I need?

[oshawa-connection]

FWIW, I'm also having this issue. I cannot find the implementation in this repo that looks at the prj stream. I'm guessing that if the shapefile has a custom projection, you won't be able to specify the SRID to NetTopologySuite (as there won't be a standard one).

Assuming no custom projections: you can read in the .prj file, its just stored as a plain old PROJ string which can be used to lookup the SRID. E.g.

• 4326 : "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.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]" "+proj=longlat +datum=WGS84 +no_defs "

• 27700 : "PROJCS["OSGB 1936 / British National Grid",GEOGCS["OSGB 1936",DATUM["OSGB_1936",SPHEROID["Airy 1830",6377563.396,299.3249646,AUTHORITY["EPSG","7001"]],TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489],AUTHORITY["EPSG","6277"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4277"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","27700"]]" "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs "

This is an ultra janky solution but it could work.

If there are custom projections, one possible idea that I haven't looked into much is https://github.com/NetTopologySuite/ProjNet4GeoAPI which can interpret PROJ strings and re-project. https://gis.stackexchange.com/questions/165022/transforming-point-using-nettopologysuite

@airbreather How accurate is this?