TransWikia.com

What is wrong with this proj4 string that the raster package in R will not accept?

Geographic Information Systems Asked by user168982 on December 27, 2020

I am trying to use a .tif file from USDA’s CropScape, which has the proj4 string

"+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0
+datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

Creating the original raster object succeeds and is plottable, but comes along with the error message

"Warning message: In .newCRS(value) : +proj=aea +lat_1=29.5
+lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 is not a valid PROJ.4 CRS string"

Likewise, trying to manually assign the CRS by typing

crs(cropscape) <- "proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

does not work and returns a similar error message, that it is not a valid PROJ.4 CRS string.

I’m no CRS or PROJ4 genius but I do not see what is wrong with this string. What is the problem here?

For reproducability: this is the code being used to generate the cropscape raster object:

cropscape <- GET(url="https://nassgeodata.gmu.edu",
                 path="axis2/services/CDLService/GetCDLFile?year=2014&fips=06") %>%
  content('text') %>%
  str_remove(".*<returnURL>") %>%
  str_remove("</returnURL>.*") %>%
  raster()

One Answer

Your CRS:

> crs(r) <- "proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
Warning message:
In .newCRS(value) :
  proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 is not a valid PROJ.4 CRS string

Works if you add a + sign at the start:

> crs(r) <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
> 

Interestingly the sf package seems to cope with it:

> pts = st_as_sf(data.frame(x=1,y=1),coords=1:2)
> st_crs(pts) <- "proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

but is silently setting the CRS to NA:

> pts
Simple feature collection with 1 feature and 0 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 1 ymin: 1 xmax: 1 ymax: 1
CRS:            NA

Some of this behaviour might be different on newer versions of GDAL with newer versions of PROJ - lots of CRS handling has changed.

Answered by Spacedman on December 27, 2020

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP