2

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() 
7
  • Where does GET and the %>% operator come from? Commented Aug 23, 2020 at 21:45
  • Sorry, GET is from package httr and is being used to interface with the cropscape API. %>% is a piping operator from dplyr. Commented Aug 23, 2020 at 21:59
  • and str_remove? If you are going for reproducibility you do need to tell us all the add-on packages. Commented Aug 23, 2020 at 22:17
  • from stringr, also a part of tidyverse. Sorry for not including this off the bat, I was under the impression that tidyverse was universal enough that these would be recognized. GET pulls the info from the API, content strips everything but the text, which is a string from which I want just a url that is surrounded by <returnURL> and </returnURL> tags Commented Aug 23, 2020 at 22:35
  • No, its not universal! And anyway, doing library(tidyverse) is never a good idea - best practice is to always use as few packages as possible and get each one with a separate library call. Commented Aug 24, 2020 at 7:15

1 Answer 1

1

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.

1
  • I think I made a mistake transcribing this to stackexchange. My original code does have a leading plus sign and still generates the error ``` 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" 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 Commented Aug 23, 2020 at 21:59

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.