search

conda-forge / rasterio-feedstock

A conda-smithy repository for rasterio.
BSD 3-Clause "New" or "Revised" License
42 stars 22 forks source link

to_epsg() returns None for some WKT strings #181

Closed forkozi closed 3 years ago

forkozi commented 3 years ago

Issue: The following code returns None for some WKT strings (of common coordinate reference systems)

from rasterio.crs import CRS
CRS.from_wkt(wkt_str).to_epsg()


This WKT string results in an EPSG code...:

PROJCS["NAD83(2011) / UTM zone 18N",GEOGCS["NAD83(2011)",DATUM["NAD83_National_Spatial_Reference_System_2011",SPHEROID["GRS 1980",6378137,298.257222101004,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","1116"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-75],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","6347"]]


This WKT string results in None:

PROJCRS["NAD83(2011) / UTM zone 17N",BASEGEOGCRS["NAD83(2011)",DATUM["NAD83 (National Spatial Reference System 2011)",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",6319]],CONVERSION["UTM zone 17N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-81,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["meter",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["meter",1],ID["EPSG",8807]],ID["EPSG",16017]],CS[Cartesian,3],AXIS["x",east,ORDER[1],LENGTHUNIT["meter",1,ID["EPSG",9001]]],AXIS["y",north,ORDER[2],LENGTHUNIT["meter",1,ID["EPSG",9001]]],AXIS["up",up,ORDER[3],LENGTHUNIT["meter",1,ID["EPSG",9001]]]]

Does rasterio.crs support v2.0.6 of the OGC Geographic information — Well-known text representation of coordinate reference systems specifications (here)? Could it be the GEOGCS-versus-BASEGEOGCRS difference that's causing the issue?


Environment (conda list):

``` $ conda list # Name Version Build Channel affine 2.3.0 py_0 conda-forge attrs 20.3.0 pyhd3deb0d_0 conda-forge boost-cpp 1.74.0 hd4e6614_0 conda-forge bzip2 1.0.8 he774522_3 conda-forge ca-certificates 2020.11.8 h5b45459_0 conda-forge certifi 2020.11.8 py37h03978a9_0 conda-forge cfitsio 3.470 h0af3d06_7 conda-forge click 7.1.2 pyh9f0ad1d_0 conda-forge click-plugins 1.1.1 py_0 conda-forge cligj 0.7.0 py_0 conda-forge curl 7.71.1 h4b64cdc_8 conda-forge expat 2.2.9 h33f27b4_2 conda-forge freetype 2.10.1 vc14_0 [vc14] esri freexl 1.0.5 hd288d7e_1002 conda-forge geos 3.8.1 he025d50_0 conda-forge geotiff 1.6.0 h8884d1a_3 conda-forge gettext 0.19.8.1 hfbb10ce_1004 conda-forge glib 2.66.2 h0e60522_0 conda-forge hdf4 4.2.13 hf8e6fe8_1003 conda-forge hdf5 1.10.6 nompi_h89124ea_1110 conda-forge icu 67.1 h33f27b4_0 conda-forge intel-openmp 2020.3 h57928b3_311 conda-forge jpeg 9d h8ffe710_0 conda-forge kealib 1.4.13 h3b59ab9_1 conda-forge krb5 1.17.2 hbae68bd_0 conda-forge libblas 3.8.0 21_mkl conda-forge libcblas 3.8.0 21_mkl conda-forge libcurl 7.71.1 h4b64cdc_8 conda-forge libffi 3.2.1 ha925a31_1007 conda-forge libgdal 3.1.4 h0e5aa5a_0 conda-forge libglib 2.66.2 h35efcdc_0 conda-forge libiconv 1.16 he774522_0 conda-forge libkml 1.3.0 he9e54da_1012 conda-forge liblapack 3.8.0 21_mkl conda-forge libnetcdf 4.7.4 nompi_h2ee746f_106 conda-forge libpng 1.6.37 h1d00b33_2 conda-forge libpq 12.3 hd9aa61d_2 conda-forge libspatialite 5.0.0 hf693123_0 conda-forge libssh2 1.9.0 hb06d900_5 conda-forge libtiff 4.1.0 hc10be44_6 conda-forge libwebp-base 1.1.0 h8ffe710_3 conda-forge libxml2 2.9.10 h1006b36_2 conda-forge lz4-c 1.9.2 h62dcd97_2 conda-forge m2w64-gcc-libgfortran 5.3.0 6 m2w64-gcc-libs 5.3.0 7 m2w64-gcc-libs-core 5.3.0 7 m2w64-gmp 6.1.0 2 m2w64-libwinpthread-git 5.0.0.4634.697f757 2 mkl 2020.4 hb70f87d_311 conda-forge msys2-conda-epoch 20160418 1 numpy 1.19.4 py37hd20adf4_1 conda-forge openjpeg 2.3.1 h57dd2e7_3 conda-forge openssl 1.1.1h he774522_0 conda-forge pcre 8.44 ha925a31_0 conda-forge pip 20.2.4 py37haa95532_0 poppler 0.89.0 h0cd1227_0 conda-forge poppler-data 0.4.10 0 conda-forge postgresql 12.3 he14cc48_2 conda-forge proj 7.1.1 h7d85306_3 conda-forge pyparsing 2.4.7 pyh9f0ad1d_0 conda-forge python 3.7.9 h60c2a47_0 python_abi 3.7 1_cp37m conda-forge rasterio 1.1.8 py37hc4b0cd6_0 conda-forge setuptools 50.3.1 py37haa95532_1 snuggs 1.4.7 py_0 conda-forge sqlite 3.33.0 h2a8f88b_0 tiledb 2.1.2 h968eb34_1 conda-forge tk 8.6.10 he774522_1 conda-forge vc 14.1 h0510ff6_4 vs2015_runtime 14.16.27012 hf0eaf9b_0 esri wheel 0.35.1 pyhd3eb1b0_0 wincertstore 0.2 py37_0 xerces-c 3.2.3 ha925a31_1 conda-forge xz 5.2.5 h62dcd97_1 conda-forge zlib 1.2.11 h62dcd97_4 zstd 1.4.5 h1f3a1b7_2 conda-forge ```


Details about conda and system ( conda info ):

``` $ conda info active environment : rasterio_test active env location : C:\Users\Nick.Forfinski-Sarko\Anaconda3\envs\rasterio_test shell level : 2 user config file : C:\Users\Nick.Forfinski-Sarko\.condarc populated config files : C:\Users\Nick.Forfinski-Sarko\.condarc conda version : 4.9.2 conda-build version : 3.18.11 python version : 3.8.3.final.0 virtual packages : __cuda=10.2=0 __win=0=0 __archspec=1=x86_64 base environment : C:\Users\Nick.Forfinski-Sarko\Anaconda3 (writable) channel URLs : https://conda.anaconda.org/esri/win-64 https://conda.anaconda.org/esri/noarch https://repo.anaconda.com/pkgs/main/win-64 https://repo.anaconda.com/pkgs/main/noarch https://repo.anaconda.com/pkgs/r/win-64 https://repo.anaconda.com/pkgs/r/noarch https://repo.anaconda.com/pkgs/msys2/win-64 https://repo.anaconda.com/pkgs/msys2/noarch https://conda.anaconda.org/conda-forge/win-64 https://conda.anaconda.org/conda-forge/noarch package cache : C:\Users\Nick.Forfinski-Sarko\Anaconda3\pkgs C:\Users\Nick.Forfinski-Sarko\.conda\pkgs C:\Users\Nick.Forfinski-Sarko\AppData\Local\conda\conda\pkgs envs directories : C:\Users\Nick.Forfinski-Sarko\Anaconda3\envs C:\Users\Nick.Forfinski-Sarko\.conda\envs C:\Users\Nick.Forfinski-Sarko\AppData\Local\conda\conda\envs platform : win-64 user-agent : conda/4.9.2 requests/2.24.0 CPython/3.8.3 Windows/10 Windows/10.0.17763 administrator : False netrc file : None offline mode : False ```
snowman2 commented 3 years ago

There are many WKT strings out there that don't have a a corresponding EPSG code. The EPSG code is just a specific CRS that has been standardized.

snowman2 commented 3 years ago

Actually, that one works with pyproj:

>>> from pyproj import CRS
>>> cc = CRS('PROJCS["NAD83(2011) / UTM zone 18N",GEOGCS["NAD83(2011)",DATUM["NAD83_National_Spatial_Reference_System_2011",SPHEROID["GRS 1980",6378137,298.257222101004,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","1116"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-75],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","6347"]]')
>>> cc
<Projected CRS: EPSG:6347>
Name: NAD83(2011) / UTM zone 18N
Axis Info [cartesian]:
- [east]: Easting (metre)
- [north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: UTM zone 18N
- method: Transverse Mercator
Datum: NAD83 (National Spatial Reference System 2011)
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

>>> cc.to_epsg()
6347
snowman2 commented 3 years ago

Works for me with rasterio:

>>> from rasterio.crs import CRS
>>> cc = CRS.from_wkt('PROJCS["NAD83(2011) / UTM zone 18N",GEOGCS["NAD83(2011)",DATUM["NAD83_National_Spatial_Reference_System_2011",SPHEROID["GRS 1980",6378137,298.257222101004,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","1116"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-75],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","6347"]]')
>>> cc.to_epsg()
6347
snowman2 commented 3 years ago

Ah, it is the second WKT:

>>> cc = pyproj.CRS('PROJCRS["NAD83(2011) / UTM zone 17N",BASEGEOGCRS["NAD83(2011)",DATUM["NAD83 (National Spatial Reference System 2011)",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",6319]],CONVERSION["UTM zone 17N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-81,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["meter",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["meter",1],ID["EPSG",8807]],ID["EPSG",16017]],CS[Cartesian,3],AXIS["x",east,ORDER[1],LENGTHUNIT["meter",1,ID["EPSG",9001]]],AXIS["y",north,ORDER[2],LENGTHUNIT["meter",1,ID["EPSG",9001]]],AXIS["up",up,ORDER[3],LENGTHUNIT["meter",1,ID["EPSG",9001]]]]')
>>> cc.to_epsg()
>>> cc.to_epsg(0)
6346
snowman2 commented 3 years ago

https://gis.stackexchange.com/questions/326690/explaining-pyproj-to-epsg-min-confidence-parameter/326919

forkozi commented 3 years ago

Terrific, thanks!