mikejohnson51
commented
2 years ago @rmcd-mscb thanks! I think this comes down to how the datum is represented as these are both the "same" projection (which is shown at the end). First though, I'll show were they come from:
Seeming Inconsistency
Looking at the "global attributes" of the file here, we would expect EPSG:4326 as you did:
However when looking at the CRS variable of the same file, there is an inverse flattening and semi-major axis provided:
The +f flag in the PROJ string is the "flattening dimension" and "un-inverting" the diminsion gives:
1/298.257223563 = 0.00335281
Equally, the +a flag is populated with the semi-major axis provided in the CRS. The code I use to detect the proj string sees these values and applies them accordingly.
Kinda the same though...
Looking at the "official" WGS84
We see the title indicates "EPSG:4326"
The proj4string is
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
But the OGC WKT is:
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.01745329251994328,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]]Note the AUTHORITY["EPSG","4326"]] and SPHEROID["WGS 84",6378137,298.257223563,
PROJ should handle this...
With the R bindings to PROJ, this is handled "behind the scenes":
library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.4.0, PROJ 8.1.1; sf_use_s2() is TRUE
library(AOI)
crs = "+proj=longlat +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs"
# look at the CRS argument (contains flattening and semi-major axis)
(sf_version = aoi_get(state = "CO") |>
st_transform(crs))
#> Simple feature collection with 1 feature and 14 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: -109.0602 ymin: 36.99246 xmax: -102.0415 ymax: 41.00342
#> CRS: +proj=longlat +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs
#> state_region state_division feature_code state_name state_abbr name
#> 1 4 8 1779779 Colorado CO Colorado
#> fip_class tiger_class combined_area_code metropolitan_area_code
#> 1 <NA> G4000 NA <NA>
#> functional_status land_area water_area fip_code
#> 1 A 268418796417 1185716938 08
#> geometry
#> 1 MULTIPOLYGON (((-105.155 36...
# the WKT CRS uses these to define the elipsoid!
st_crs(sf_version)
#> Coordinate Reference System:
#> User input: +proj=longlat +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs
#> wkt:
#> GEOGCRS["unknown",
#> DATUM["unknown",
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1,
#> ID["EPSG",9001]]]],
#> PRIMEM["unknown",0,
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]],
#> CS[ellipsoidal,2],
#> AXIS["longitude",east,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]],
#> AXIS["latitude",north,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]]]
# But is dropped in the proj4string representation:
st_crs(sf_version)$proj4string
#> [1] "+proj=longlat +ellps=WGS84 +no_defs"Created on 2022-04-05 by the reprex package (v2.0.1)
So, long story short I am surprised pyproj is choking on that PROJ string as it believe it is valid, and "correct".
I am not sure how to avoid the inclusion of CRS elements in the automated catalog building and hope to leave that to the spatial backends used.
Hopefully we can find a way to make this work with python!
rmcd-mscb
Mike I think the "proj" grid value is wrong for gridmet. If i try to use it to reproject a shapefile it fails. So I checked it as follows against the stated epsg code 4326:
import pyproj test = pyproj.CRS.from_user_input(4326) test.to_proj4()resulting in '+proj=longlat +datum=WGS84 +no_defs +type=crs'
as compared to:
+proj=longlat +a=6378137 +f=0.0033528106647474...
when I use the proj value in the grid file if doesn't like the "f" value.