0

I have downloaded the following data as a GeoTIFF and I'm interested in determining the value from this raster at specific locations.

raster = readGDAL("filename.TIFF")
summary(raster)
Object of class SpatialGridDataFrame
Coordinates:
   min max
x -180 180
y  -90  90
Is projected: FALSE 
proj4string :
[+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
Grid attributes:
  cellcentre.offset cellsize cells.dim
x            -179.5        1       360
y             -89.5        1       180
Data attributes:
     band1      
 Min.   :  1.0  
 1st Qu.:255.0  
 Median :255.0  
 Mean   :226.9  
 3rd Qu.:255.0  
 Max.   :255.0  
 NA's   :245    

head(data[7:8]) #Xp = longitude, Yp = latitude
              Xp        Yp
    7  -50.41777   1.33138
    8  -59.25075 -19.16780
    9  -67.25043 -27.83412
    10 -54.91759 -15.00130
    12 -67.58375 -28.00078
    13 -65.08385 -34.33386

However, when I run the following, I get an error message:

sp <- SpatialPoints(data[7:8])
summary(sp)
Object of class SpatialPoints
Coordinates:
         min       max
Xp -76.08341 -35.25171
Yp -38.66702   1.33138
Is projected: NA 
proj4string : [NA]
Number of points: 184

data$ndvi<-extract(raster, sp)
    Error in .subset(x, j) : invalid subscript type 'S4'

If I use simulated data for the raster, I don't seem to have the same issue:

simr <- raster(ncol=36, nrow=18)
simr[] <- 1:ncell(simr)
simndvi<-extract(raster, sp)
head(simndvi)
[1] 301 373 408 373 408 444

How can I figure out why this isn't working for the GeoTIFF?

  • Show all your working, including how you read the tiff in, what the summary function says on the object you've read the tiff into, and the sp object. – Spacedman Nov 16 '18 at 19:54
  • Have you got other packages loaded, such as tdiyr, that mask the extract function? Start with a fresh session and make a reproducible example. – Spacedman Nov 16 '18 at 19:55
  • @Spacedman I have included the summaries you requested, and no, I am not using other packages that mask extract (and yes, I did start with a fresh session to double check). – tnt Nov 16 '18 at 20:05
  • It might help if you provide a CRS for your spatial points (WGS84 I assume). – AndreJ Nov 16 '18 at 20:14
  • @AndreJ how do I do that? – tnt Nov 16 '18 at 20:16
1

You are using readGDAL which returns a SpatialGridDataFrame - instead use the raster function from the raster package to create a Raster object. Then extract should work.

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.