GeoscienceAustralia / ptha

Probabilistic Tsunami Hazard Assessment from major earthquake source-zones
BSD 3-Clause "New" or "Revised" License
20 stars 9 forks source link

Changes in CRS("+init=epsg:4326") #23

Closed gareth-d-ga closed 3 months ago

gareth-d-ga commented 1 year ago

Recent updates to dependencies of the R spatial ecosystem can potentially break code that refers to WGS84 as +init=epsg:4326. We get a warning like:

In CPL_crs_from_input(x) :
  GDAL Message 1: +init=epsg:XXXX syntax is deprecated. It might return a CRS with a non-EPSG compliant axis order.

The updates affect axis order. Historically R would always treat this case as lon/lat ordered, but apparently the EPSG standard implies lat/lon ordering.

Need to figure out whether workarounds are required. At one point I thought I hit a bug related to this, but it seems to have been something else.

See also here: https://cran.r-project.org/web/packages/rgdal/vignettes/CRS_projections_transformations.html

gareth-d-ga commented 1 year ago

The aforementioned website says "It will be useful to know that in general OGC:CRS84 should be used instead of EPSG:4326, because the latter presupposes that Latitude is always ordered before Longitude. OGC:CRS84 is the standard representation used by GeoJSON, with coordinates always ordered Longitude before Latitude. This is also the standard GIS and visualization order:"

gareth-d-ga commented 1 year ago

See also this blogpost under the heading "Axis order".

R-spatial packages have, for the past 25 years, pretty much assumed that two-dimensional data are XY-ordered, or longitude-latitude. Geodesists, on the other hand, typically use (ϕ,λ), or latitude-longitude, as coordinate pairs....

Further down:

Package sf by default uses a switch in GDAL that brings everything in the old, longitude-latitude order, but data may come in [in another ordering](https://github.com/r-spatial/sf/issues/1245).

Workflows using sp/rgdal should expect “GIS style” axis order to be preserved

So it sounds like there should not be any problems?....

gareth-d-ga commented 1 year ago

Beware that newer versions of R+packages will accept the notation CRS('epsg:4326') but older versions (e.g. like I'm using on NCI) will not.

However it seems both the newer and older versions accept the equivalent in capitals: CRS("EPSG:4326")

gareth-d-ga commented 3 months ago

This seems to have gone smoothly so closing.