Tuesday 21 August 2018

raster - CRS validity and PROJ versions


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.




No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...