wjbenfold
commented
2 years ago Closed wjbenfold closed 2 years ago
wjbenfold
commented
2 years ago 
snowman2
commented
2 years ago Looking at this from the perspective of creating a Datum for a CRS may help: https://pyproj4.github.io/pyproj/stable/build_crs.html
When you specify +datum=WGS84 it could be a datum or an ensemble of datums. Each datum implicitly has a prime meridian & ellipsoid already defined:
>>> from pyproj.crs import Datum
>>> Datum.from_string("WGS84")
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)",
ID["EPSG",1166]],
MEMBER["World Geodetic System 1984 (G730)",
ID["EPSG",1152]],
MEMBER["World Geodetic System 1984 (G873)",
ID["EPSG",1153]],
MEMBER["World Geodetic System 1984 (G1150)",
ID["EPSG",1154]],
MEMBER["World Geodetic System 1984 (G1674)",
ID["EPSG",1155]],
MEMBER["World Geodetic System 1984 (G1762)",
ID["EPSG",1156]],
MEMBER["World Geodetic System 1984 (G2139)",
ID["EPSG",1309]],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1],
ID["EPSG",7030]],
ENSEMBLEACCURACY[2.0],
ID["EPSG",6326]]
>>> Datum.from_string("WGS84").name
'World Geodetic System 1984 ensemble'
>>> Datum.from_epsg(1309).ellipsoid
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1],
ID["EPSG",7030]]
>>> Datum.from_epsg(1309).prime_meridian
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]When you are specifying the +datum as well as ellipsoid parameters, it has to pick either the datum or the ellipsoid parameters to create the datum as it cannot override the existing datum presets. In this scenario, it is using the ellipsoid parameters:
>>> from pyproj.crs.datum import CustomDatum, CustomEllipsoid
>>> ell = CustomEllipsoid(semi_major_axis=6378137, inverse_flattening=298.257223563)
>>> ell
ELLIPSOID["undefined",6378137,298.257223563,
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]
>>> datum = CustomDatum(ellipsoid=ell)
>>> datum
DATUM["undefined",
ELLIPSOID["undefined",6378137,298.257223563,
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]However, it has figured out the WGS84 ellipse in your scenario:
>>> from pyproj import CRS
>>> cc = CRS("+proj=merc +a=6378137.0 +rf=298.257223563 +datum=WGS84 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +units=m +no_defs +type=crs")
>>> cc.datum
DATUM["unknown",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
The PROJ string parsing logic is handled by PROJ (https://github.com/OSGeo/PROJ/). I recommend raising an issue upstream if you have further questions or think there is an issue that needs to be addressed.
wjbenfold
commented
2 years ago Thanks @snowman2, I'll take it up with PROJ. I think there's an issue to be addressed in the way that information is silently dropped without indication that the provided string was invalid, but I've not been working with these strings for very long so would be interested to know if you think differently?
snowman2
commented
2 years ago I think there's an issue to be addressed in the way that information is silently dropped without indication that the provided string was invalid
I think that would something to discuss with the PROJ devs (https://github.com/OSGeo/PROJ/). They have more background on the history and how that might impact other applications using PROJ.
snowman2
commented
2 years ago This is also a helpful reference (note how the +datum= is dropped):
>>> from pyproj import CRS
>>> cc = CRS("+proj=merc +a=6378137.0 +rf=298.257223563 +datum=WGS84 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +units=m +no_defs +type=crs")
>>> cc.to_proj4()
UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
return self._crs.to_proj4(version=version)
'+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs'
Code Sample, a copy-pastable example if possible
Problem description
Currently, passing the datum along with ellipsoid parameters (even accurate ones) when creating a CRS causes the datum to be silently ignored by pyproj. The datum is still visible when printing the CRS, but specifically interrogating the datum shows that it is not attached (and indeed transformations using the CRS that has been created will act as if the datum is not present).
Expected Output
I would be happy for the two examples above to be interpreted as equivalent, or for the user to be warned that the datum is being ignored (along with the datum not being visible when the CRS is printed). This would make it easier to work out what has gone wrong when this instance is encountered. There could alternatively be an error raised if this definition is considered invalid because a datum specifies the ellipse, or the ellipse could be modified using the parameters given (in the same way that providing parameters after an ellipse specification works).
Environment Information
Installation method
Conda
Conda environment information (if you installed with conda):
Environment (
conda list):Details about
condaand system (conda info):