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()
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
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP