2

I'm not sure where the right place to raise this is, so:

The PROJ strings for EPSG codes 4151 and 4283 are identical:

+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 

which is a hassle because one is Switzerland and one is Australia. This causes an issue when writing geospatial data to file from R with sf, raster or sp and then trying to open it in other software like QGIS. 4151 is usually autodetected as the file CRS, presumably because the EPSG code shows up first in a search.

With PROJ 4.9.1 on Windows, I thought I'd found a way around this by including +init=EPSG:4283 in the string:

+init=EPSG:4283 +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 

However, with PROJ 5.2.0 on Ubuntu 18.04 LTS, this string is considered invalid and an error, 'no arguments in initialisation list' is returned. This is encountered both with sp::CRS() and sf::st_crs() and persists even when +init is used alone, e.g.

sf::st_crs('+init=EPSG:4283')

Validity appears to be case-sensitive - If I run

sf::st_crs('+init=epsg:4283')

no error is returned in Ubuntu. However, if I do that under Windows, I get inconsistent results:

> st_crs('+init=EPSG:4283')
Coordinate Reference System:
  EPSG: 4283 
  proj4string: "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"
> st_crs('+init=epsg:4283')
Coordinate Reference System:
  EPSG: 4283 
  proj4string: "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"
 > st_crs('+init=EPSG:4283 +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs')
Coordinate Reference System:
  EPSG: 4283 
  proj4string: "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"
> st_crs('+init=epsg:4283 +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs')
Coordinate Reference System:
  No EPSG code
  proj4string: "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"

Note the missing EPSG code.

Is my best option as an R package developer to force the init tag to lower-case if a Linux-based OS is detected? I don't think I have an easy way to detect the PROJ version in use.

  • 1
    Use sf_extSoftVersion() to get versions of stuff sf links to. – Spacedman Jul 6 at 13:41
  • Thanks, I think that's preferable to OS detection. – obrl_soil Jul 7 at 2:49
1

I wonder what is the real problem if the EPSG number is interpreted wrong because both CRS definitions mean the same and there will be no error in coordinate operations, at least if you do not need to take the time dimension and movements of the continents into account. However, GDAL 3.0 which is using Proj v.6 means an improvement because coordinate system definitions are using OGC WKT:2018. If you are R package developer I believe you should consentrate on upgrading to Proj v.6.just as GDAL did https://gdalbarn.com/

gdalsrsinfo epsg:4151

PROJ.4 : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs

OGC WKT2:2018 :
BOUNDCRS[
    SOURCECRS[
        GEOGCRS["CHTRF95",
            DATUM["Swiss Terrestrial Reference Frame 1995",
                ELLIPSOID["GRS 1980",6378137,298.257222101,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["geodetic latitude (Lat)",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["geodetic longitude (Lon)",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            USAGE[
                SCOPE["unknown"],
                AREA["Europe - Liechtenstein and Switzerland"],
                BBOX[45.82,5.96,47.81,10.49]],
            ID["EPSG",4151]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["latitude",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["longitude",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    ABRIDGEDTRANSFORMATION["CHTRF95 to WGS 84 (1)",
        METHOD["Geocentric translations (geog2D domain)",
            ID["EPSG",9603]],
        PARAMETER["X-axis translation",0,
            ID["EPSG",8605]],
        PARAMETER["Y-axis translation",0,
            ID["EPSG",8606]],
        PARAMETER["Z-axis translation",0,
            ID["EPSG",8607]],
        USAGE[
            SCOPE["unknown"],
            AREA["Europe - Liechtenstein and Switzerland"],
            BBOX[45.82,5.96,47.81,10.49]],
        ID["EPSG",1511]]]


gdalsrsinfo epsg:4283

PROJ.4 : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs

OGC WKT2:2018 :
BOUNDCRS[
    SOURCECRS[
        GEOGCRS["GDA94",
            DATUM["Geocentric Datum of Australia 1994",
                ELLIPSOID["GRS 1980",6378137,298.257222101,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["geodetic latitude (Lat)",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["geodetic longitude (Lon)",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            USAGE[
                SCOPE["unknown"],
                AREA["Australia - GDA"],
                BBOX[-60.56,93.41,-8.47,173.35]],
            ID["EPSG",4283]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["latitude",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["longitude",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    ABRIDGEDTRANSFORMATION["GDA94 to WGS 84 (1)",
        METHOD["Geocentric translations (geog2D domain)",
            ID["EPSG",9603]],
        PARAMETER["X-axis translation",0,
            ID["EPSG",8605]],
        PARAMETER["Y-axis translation",0,
            ID["EPSG",8606]],
        PARAMETER["Z-axis translation",0,
            ID["EPSG",8607]],
        USAGE[
            SCOPE["unknown"],
            AREA["Australia - GDA"],
            BBOX[-60.56,93.41,-8.47,173.35]],
        ID["EPSG",1150]]]
  • Yes the ideal is that r-spatial packages move to the same stack cross-platform (at least in terms of major releases), but this is beyond my control. – obrl_soil Jul 7 at 2:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.