wfletcher / EasyGIS.NET

EasyGIS.NET source
GNU Lesser General Public License v3.0
68 stars 24 forks source link

Example of Two .PRJ Files Declaring EPSG 27700 #125

Closed dusko23 closed 12 months ago

dusko23 commented 12 months ago

There is a caveat in creating shp.CoordinateReferenceSystem. I have uploaded two prj files. Id is set to Nothing in the case of one prj file and set correctly to 27700 for the other. Both shape files are downloaded from OSG Open Data HUB. OSGB_EPEG_27700_error.zip

dusko23 commented 12 months ago

This "nice try" did not work either.

image

Id remains Nothing. Even if it would work, shp.CoordinateReferenceSystem is Read Only.

wfletcher commented 12 months ago

I have had a look at your prj files and cannot reproduce your issue. Creating a CRS from both files identify with ID 27700

Creating a CRS from ESRI WKT (which does not include the 27700 ID) is correctly identified with ID 27700.

PROJCS["British_National_Grid",GEOGCS["GCS_OSGB_1936",DATUM["D_OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",400000.0],PARAMETER["False_Northing",-100000.0],PARAMETER["Central_Meridian",-2.0],PARAMETER["Scale_Factor",0.9996012717],PARAMETER["Latitude_Of_Origin",49.0],UNIT["Meter",1.0]]


 [Test]
        public void TestEPSG27700WithDifferentWktAreEqual()
        {
            string wkt1 = @"PROJCRS[""OSGB36 / British National Grid"",BASEGEOGCRS[""OSGB36"",DATUM[""Ordnance Survey of Great Britain 1936"",ELLIPSOID[""Airy 1830"",6377563.396,299.3249646]],UNIT[""degree"",0.0174532925199433]],CONVERSION[""British National Grid"",METHOD[""Transverse Mercator""],PARAMETER[""Latitude of natural origin"",49],PARAMETER[""Longitude of natural origin"",-2],PARAMETER[""Scale factor at natural origin"",0.9996012717],PARAMETER[""False easting"",400000],PARAMETER[""False northing"",-100000]],CS[Cartesian,2],AXIS[""(E)"",east],AXIS[""(N)"",north],UNIT[""metre"",1],USAGE[SCOPE[""Engineering survey, topographic mapping.""],AREA[""United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore.""],BBOX[49.75,-9,61.01,2.01]],ID[""EPSG"",27700]]";
            string wkt2 = @"PROJCS[""British_National_Grid"",GEOGCS[""GCS_OSGB_1936"",DATUM[""D_OSGB_1936"",SPHEROID[""Airy_1830"",6377563.396,299.3249646]],PRIMEM[""Greenwich"",0.0],UNIT[""Degree"",0.0174532925199433]],PROJECTION[""Transverse_Mercator""],PARAMETER[""False_Easting"",400000.0],PARAMETER[""False_Northing"",-100000.0],PARAMETER[""Central_Meridian"",-2.0],PARAMETER[""Scale_Factor"",0.9996012717],PARAMETER[""Latitude_Of_Origin"",49.0],UNIT[""Meter"",1.0],AUTHORITY[""EPSG"",27700]]";

            var crs1 = EGIS.Projections.CoordinateReferenceSystemFactory.Default.CreateCRSFromWKT(wkt1);
            var crs2 = EGIS.Projections.CoordinateReferenceSystemFactory.Default.CreateCRSFromWKT(wkt2);

            Assert.Multiple(() =>
            {
                Assert.AreEqual(crs1.Id, "27700", "crs1 Id should be 27700");
                Assert.AreEqual(crs2.Id, "27700", "crs2 Id should be 27700");
                Assert.IsTrue(crs1.IsEquivalent(crs2), "crs1 and crs2 should be equivalent");
            });

            //test recreating from the CRS wkt
            string crs1Wkt = crs1.WKT;
            string crs2Wkt = crs2.WKT;

            crs1 = EGIS.Projections.CoordinateReferenceSystemFactory.Default.CreateCRSFromWKT(crs1Wkt);
            crs2 = EGIS.Projections.CoordinateReferenceSystemFactory.Default.CreateCRSFromWKT(crs2Wkt);

            Assert.Multiple(() =>
            {
                Assert.AreEqual(crs1.Id, "27700", "crs1 Id should be 27700");
                Assert.AreEqual(crs2.Id, "27700", "crs2 Id should be 27700");
                Assert.IsTrue(crs1.IsEquivalent(crs2), "crs1 and crs2 should be equivalent");
            });

            //test using ESRI WKT format
            string esriWkt = crs1.GetWKT(PJ_WKT_TYPE.PJ_WKT1_ESRI, true);
            Console.Out.WriteLine(esriWkt);
            crs1 = EGIS.Projections.CoordinateReferenceSystemFactory.Default.CreateCRSFromWKT(esriWkt);

            esriWkt = crs2.GetWKT(PJ_WKT_TYPE.PJ_WKT1_ESRI, true);
            crs2 = EGIS.Projections.CoordinateReferenceSystemFactory.Default.CreateCRSFromWKT(esriWkt);

            Assert.Multiple(() =>
            {
                Assert.AreEqual(crs1.Id, "27700", "crs1 Id should be 27700");
                Assert.AreEqual(crs2.Id, "27700", "crs2 Id should be 27700");
                Assert.IsTrue(crs1.IsEquivalent(crs2), "crs1 and crs2 should be equivalent");
            });

            //test using GDAL WKT format

            string gdalWkt = crs1.GetWKT(PJ_WKT_TYPE.PJ_WKT1_GDAL, false);
            crs1 = EGIS.Projections.CoordinateReferenceSystemFactory.Default.CreateCRSFromWKT(gdalWkt);

            gdalWkt = crs2.GetWKT(PJ_WKT_TYPE.PJ_WKT1_GDAL, true);
            crs2 = EGIS.Projections.CoordinateReferenceSystemFactory.Default.CreateCRSFromWKT(gdalWkt);

            Assert.Multiple(() =>
            {
                Assert.AreEqual(crs1.Id, "27700", "crs1 Id should be 27700");
                Assert.AreEqual(crs2.Id, "27700", "crs2 Id should be 27700");
                Assert.IsTrue(crs1.IsEquivalent(crs2), "crs1 and crs2 should be equivalent");
            });

        }
dusko23 commented 12 months ago

I see you have hardcoded WKT. Please try and read from the files. The files are obtained from OSGB authority.

wfletcher commented 12 months ago

There is no difference. Are you sure you attached the files you tested with?

string prjFile1 = System.IO.Path.Combine(System.IO.Path.GetDirectoryName(Assembly.GetExecutingAssembly().Location), "data", "27700", "TV69_line_ERROR.prj");
            string prjFile2 = System.IO.Path.Combine(System.IO.Path.GetDirectoryName(Assembly.GetExecutingAssembly().Location), "data", "27700", "WatercourseLink_OK.prj");
            crs1 = EGIS.Projections.CoordinateReferenceSystemFactory.Default.CreateCRSFromPrjFile(prjFile1);
            crs2 = EGIS.Projections.CoordinateReferenceSystemFactory.Default.CreateCRSFromPrjFile(prjFile2);
            Assert.Multiple(() =>
            {
                Assert.AreEqual(crs1.Id, "27700", "crs1 Id should be 27700");
                Assert.AreEqual(crs2.Id, "27700", "crs2 Id should be 27700");
                Assert.IsTrue(crs1.IsEquivalent(crs2), "crs1 and crs2 should be equivalent");
            });
dusko23 commented 12 months ago

I tried using the attached file. And Id = Nothing My code:

image HU14_line.zip

Please try at you end. If it works please send me your proj9 database, I would like to find the source of the issue. If I would to deinstall your library and instal via nuget again and it works, then I would not know the cause of the issue.

dusko23 commented 12 months ago

I have uninstalled and installed the library via nuget. I have replicated your code. Problem persists. Could you confirm that you are testing with 100% the same library as currently available via the nuget package? If yes, could you please provide a c# solution that contains everything and I will test it here?

wfletcher commented 12 months ago

The 3rd file you uploaded does not identify unless you reduce the confidence theshold because it has a difference scale factor (try changing "PARAMETER["Scale_Factor",0.999601272]" to "PARAMETER["Scale_Factor",0.9996012717]").

If you add this line before you read the prj file it will identify the ID. The library uses a default value of 70 and I am not going to reduce this. EGIS.Projections.CoordinateReferenceSystemFactory.IdentificationConfidenceThreshold = 25;

Please read this link for more documentation on the confidence theshold. https://proj.org/en/9.3/development/reference/functions.html#c.proj_identify

Unfortunately you may find slight differences in shapefile prj WKT formats if you are downloading from different sources and there is really not much you can do about this. If you do encounter unknown projections and you are relying on having the ID identified you might be best to convert the shapefile into a WGS84 shapefile

dusko23 commented 12 months ago

Thanks for explaining the issue. Now I understand why are you using function isEquivalent. Currently, I do rely on ID. The thing is that isEquivalent function is very slow compared to if ID = XXX. I enable loading any EPSG type of Shapefile into my system. Converting them into WGS84 is not an option. To correctly detect the closest point, and record on mouse move I iterate through all shape files loaded and convert cursor coordinates to matching shape file EPSG (if not equal). The sole purpose is to enable editing Shapefiles whatever the SFmap object Coordinate system is. Currently, I allow editing only those Shapefiles where EPSG matches the SFmap object. The thing is that my system defines SFmap object EPSG using three options to define SFmap object EPSG (it is a must):

  1. Loading raster map background from the database.
  2. Loading raster maps from online resources as Spherical Mercator (Google Maps, OSM, ....)
  3. Defining blank raster map by selection of any of the projected ICRs.

Having ShapFiles loaded user can switch to any of the above raster map definitions. It is quite interesting to look at how Ordnance Survey Data UK data aligns when projected to Bolivia, for example.

I like to concept of EPSG. To me, it is just a definition of a "flat paper" where data is drawn.