[gdal-dev] Shapefile: automatic retrieval of EPSG code

Tue Sep 26 16:34:21 PDT 2017

This is incredibly useful! Hopefully I will get around to recompiling my
docker image to test this out soon!

On Mon, Sep 25, 2017 at 12:30 PM, Even Rouault <even.rouault at spatialys.com>
wrote:

> Hi,
>
>
>
> As part of a AWEOC (*), following a recent email thread about year-long
> issues when importing shapefiles, I've committed a first version of a
> OSRFindMatches() function that iterates through the EPSG database to check
> for one or several matches between the input SRS and the catalog.
>
>
>
> The shapefile driver now uses this function, instead of the more simple
> AutoIdentifyEPSG(), which only worked on a tiny subset of SRS. When a
> single match is found by OSRFindMatches() with an excellent match quality
> (meaning IsSame() returning TRUE), it is automatically picked up by the
> Shapefile driver.
>
>
>
> OSRFindMatches() returns both a list of candidate SRS and a confidence
> value: 0=bad, 100=perfect . Currently if several candidates pass IsSame(),
> they will be returned with 90% confidence. If a SRS doesn't pass IsSame()
> but is the single one to have a PROJCS/GEOGCS name matching the input SRS
> name (name matching is done with a few hardcoded approximation rules), it
> is returned with 50% confidence.
>
>
>
> OSRFindMatches() execution time is :
>
> - first execution ever: around 500 ms. It will then create files in
>
> ~/.gdal/$(GDAL_VERSION_MAJOR).$(GDAL_VERSION_MINOR)/srs_cache to speed-up
> later executions
>
> - first execution in a process, using the above cache: around 100ms
>
> - later executions in the process: around 10 ms per SRS (55s to identify
> the 5216 entries of GEOGCS and PROJCS from the EPSG database that are valid
> for GDAL)
>
>
>
> So I feel this is good enough to be enabled by default in the shapefile
> driver.
>
>
>
> I've tested & tuned the function:
>
> - the 5216 above mentionned SRS exported from OGC WKT to ESRI WKT, and
> reimported from ESRI WKT to OGC WKT. 100% confidence single match is
> obtained for most SRS. Except occurences using Polar_Stereographic where
> ESRI WKT lose the scale_factor and other using Hotine_Oblique_Mercator
> where the rectified_grid_angle parametr is lost. But those are probably
> issues with morphToESRI() / morphFromESRI()
>
> - the SRS definitions at
>
> http://help.arcgis.com/en/arcgisserver/10.0/apis/rest/pcs.html . Results
> are less good there. Sometimes for good reason: no entry in the EPSG code,
> different projection or datum parameters, etc.
>
>
>
> It could certainly be improved. For example by using instead of the
> back&white IsSame() method a (to-be-created) GetSimilarityLevel() method
> that would compare how different 2 SRS are. More relaxed SRS name
> comparisons would also help (as the 'tokens' making SRS names are separated
> with underscore once 'massaged', you could split on that, and then compare
> how many tokens are in common (or have partial match), and sum the number
> of identical characters, etc..).
>
>
>
>
>
> The gdalsrsinfo utility is also enhanced to use the new OSRFindMatches()
> function if you use the -e switch (for epsg).
>
>
>
> For example, with some.prj containing: PROJCS["NAD_1983_StatePlane_Te
> xas_Central_FIPS_4203_Feet",GEOGCS["GCS_North_American_1983"
> ,DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.
> 0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",
> 0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],
> PARAMETER["False_Easting",2296583.333333333],PARAMETER["
> False_Northing",9842500.0],PARAMETER["Central_Meridian",-
> 100.3333333333333],PARAMETER["Standard_Parallel_1",30.
> 11666666666667],PARAMETER["Standard_Parallel_2",31.
> 88333333333333],PARAMETER["Latitude_Of_Origin",29.
> 66666666666667],UNIT["Foot_US",0.3048006096012192]]
>
>
>
> $ gdalsrsinfo -e ESRI::some.prj
>
>
>
> EPSG:2277
>
>
>
> PROJ.4 : +proj=lcc +lat_1=31.88333333333333 +lat_2=30.11666666666667
> +lat_0=29.66666666666667 +lon_0=-100.3333333333333 +x_0=699999.9998983998
> +y_0=3000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs
>
>
>
> OGC WKT :
>
> PROJCS["NAD83 / Texas Central (ftUS)",
>
> GEOGCS["NAD83",
>
> DATUM["North_American_Datum_1983",
>
> SPHEROID["GRS 1980",6378137,298.257222101,
>
> AUTHORITY["EPSG","7019"]],
>
> TOWGS84[0,0,0,0,0,0,0],
>
> AUTHORITY["EPSG","6269"]],
>
> PRIMEM["Greenwich",0,
>
> AUTHORITY["EPSG","8901"]],
>
> UNIT["degree",0.0174532925199433,
>
> AUTHORITY["EPSG","9122"]],
>
> AUTHORITY["EPSG","4269"]],
>
> PROJECTION["Lambert_Conformal_Conic_2SP"],
>
> PARAMETER["standard_parallel_1",31.88333333333333],
>
> PARAMETER["standard_parallel_2",30.11666666666667],
>
> PARAMETER["latitude_of_origin",29.66666666666667],
>
> PARAMETER["central_meridian",-100.3333333333333],
>
> PARAMETER["false_easting",2296583.333],
>
> PARAMETER["false_northing",9842500.000000002],
>
> UNIT["US survey foot",0.3048006096012192,
>
> AUTHORITY["EPSG","9003"]],
>
> AXIS["X",EAST],
>
> AXIS["Y",NORTH],
>
> AUTHORITY["EPSG","2277"]]
>
>
>
>
>
> It can also return several occurences if there isn't a single candidate
>
>
>
> For example:
>
> $ gdalsrsinfo -e 'PROJCS["undefined",GEOGCS["ET
> RS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS
> 1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["de
> gree",0.0174532925199433]],PROJECTION["Transverse_Mercator"]
> ,PARAMETER["latitude_of_origin",0],PARAMETER["central_
> meridian",27],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",
> 500000],PARAMETER["false_northing",0],UNIT["metre",1]]'
>
>
>
> Confidence in this match: 90 %
>
>
>
> EPSG:3047
>
>
>
> PROJ.4 : +proj=utm +zone=35 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
> +no_defs
>
>
>
> OGC WKT :
>
> PROJCS["ETRS89 / UTM zone 35N (N-E)",
>
> GEOGCS["ETRS89",
>
> DATUM["European_Terrestrial_Reference_System_1989",
>
> SPHEROID["GRS 1980",6378137,298.257222101,
>
> AUTHORITY["EPSG","7019"]],
>
> TOWGS84[0,0,0,0,0,0,0],
>
> AUTHORITY["EPSG","6258"]],
>
> PRIMEM["Greenwich",0,
>
> AUTHORITY["EPSG","8901"]],
>
> UNIT["degree",0.0174532925199433,
>
> AUTHORITY["EPSG","9122"]],
>
> AUTHORITY["EPSG","4258"]],
>
> PROJECTION["Transverse_Mercator"],
>
> PARAMETER["latitude_of_origin",0],
>
> PARAMETER["central_meridian",27],
>
> PARAMETER["scale_factor",0.9996],
>
> PARAMETER["false_easting",500000],
>
> PARAMETER["false_northing",0],
>
> UNIT["metre",1,
>
> AUTHORITY["EPSG","9001"]],
>
> AUTHORITY["EPSG","3047"]]
>
>
>
> Confidence in this match: 90 %
>
>
>
> EPSG:3067
>
>
>
> PROJ.4 : +proj=utm +zone=35 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
> +no_defs
>
>
>
> OGC WKT :
>
> PROJCS["ETRS89 / TM35FIN(E,N)",
>
> GEOGCS["ETRS89",
>
> DATUM["European_Terrestrial_Reference_System_1989",
>
> SPHEROID["GRS 1980",6378137,298.257222101,
>
> AUTHORITY["EPSG","7019"]],
>
> TOWGS84[0,0,0,0,0,0,0],
>
> AUTHORITY["EPSG","6258"]],
>
> PRIMEM["Greenwich",0,
>
> AUTHORITY["EPSG","8901"]],
>
> UNIT["degree",0.0174532925199433,
>
> AUTHORITY["EPSG","9122"]],
>
> AUTHORITY["EPSG","4258"]],
>
> PROJECTION["Transverse_Mercator"],
>
> PARAMETER["latitude_of_origin",0],
>
> PARAMETER["central_meridian",27],
>
> PARAMETER["scale_factor",0.9996],
>
> PARAMETER["false_easting",500000],
>
> PARAMETER["false_northing",0],
>
> UNIT["metre",1,
>
> AUTHORITY["EPSG","9001"]],
>
> AXIS["Easting",EAST],
>
> AXIS["Northing",NORTH],
>
> AUTHORITY["EPSG","3067"]]
>
>
>
> Confidence in this match: 90 %
>
>
>
> EPSG:5048
>
>
>
> PROJ.4 : +proj=utm +zone=35 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
> +no_defs
>
>
>
> OGC WKT :
>
> PROJCS["ETRS89 / TM35FIN(N,E)",
>
> GEOGCS["ETRS89",
>
> DATUM["European_Terrestrial_Reference_System_1989",
>
> SPHEROID["GRS 1980",6378137,298.257222101,
>
> AUTHORITY["EPSG","7019"]],
>
> TOWGS84[0,0,0,0,0,0,0],
>
> AUTHORITY["EPSG","6258"]],
>
> PRIMEM["Greenwich",0,
>
> AUTHORITY["EPSG","8901"]],
>
> UNIT["degree",0.0174532925199433,
>
> AUTHORITY["EPSG","9122"]],
>
> AUTHORITY["EPSG","4258"]],
>
> PROJECTION["Transverse_Mercator"],
>
> PARAMETER["latitude_of_origin",0],
>
> PARAMETER["central_meridian",27],
>
> PARAMETER["scale_factor",0.9996],
>
> PARAMETER["false_easting",500000],
>
> PARAMETER["false_northing",0],
>
> UNIT["metre",1,
>
> AUTHORITY["EPSG","9001"]],
>
> AUTHORITY["EPSG","5048"]]
>
>
>
> Confidence in this match: 90 %
>
>
>
> EPSG:25835
>
>
>
> PROJ.4 : +proj=utm +zone=35 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
> +no_defs
>
>
>
> OGC WKT :
>
> PROJCS["ETRS89 / UTM zone 35N",
>
> GEOGCS["ETRS89",
>
> DATUM["European_Terrestrial_Reference_System_1989",
>
> SPHEROID["GRS 1980",6378137,298.257222101,
>
> AUTHORITY["EPSG","7019"]],
>
> TOWGS84[0,0,0,0,0,0,0],
>
> AUTHORITY["EPSG","6258"]],
>
> PRIMEM["Greenwich",0,
>
> AUTHORITY["EPSG","8901"]],
>
> UNIT["degree",0.0174532925199433,
>
> AUTHORITY["EPSG","9122"]],
>
> AUTHORITY["EPSG","4258"]],
>
> PROJECTION["Transverse_Mercator"],
>
> PARAMETER["latitude_of_origin",0],
>
> PARAMETER["central_meridian",27],
>
> PARAMETER["scale_factor",0.9996],
>
> PARAMETER["false_easting",500000],
>
> PARAMETER["false_northing",0],
>
> UNIT["metre",1,
>
> AUTHORITY["EPSG","9001"]],
>
> AXIS["Easting",EAST],
>
> AXIS["Northing",NORTH],
>
> AUTHORITY["EPSG","25835"]]
>
>
>
> Regards,
>
>
>
> Even
>
>
>
> (*) Autumn Week-End Of Code
>
>
>
> --
>
> Spatialys - Geospatial professional services
>
> http://www.spatialys.com
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20170926/c345c2d7/attachment-0001.html>