0

I have many point shapefiles which are not overlapping, which I would like to attribute to analagous rasters, also not overlapping. I would like to attribute these points with the raster data. For some of the raster data types with which I am doing this, I was able to merge the rasters first, then attribute.
However, my last sets of raster data do not have the same origin, so I was unable to merge/mosaic them. I am trying to instead attribute the points to rasters without merging the rasters. This would require me to use extract() on specific spatial point - raster pairs. I have named each spatial point file with a unique 4-letter name, which is also part of the raster name I would like the extract() to use.

I've created a reproducible example below that mimics my data and problem. Can anyone make suggestions on how I could code the looping for extract() to get the spatial point file to extract to the analagously-named raster?

Alternatively, if it is better/possible to combine all the spatial points and just loop through extracting all the rasters, then managing the data so that all extracted values are in one vector or column of a dataframe, perhaps that would be better.

I am using RStudio 1.2.1335

library(raster)
library(sp)   

#create point shapefiles
loc1 <- data.frame(x = c(-100,-90, -80, -70),
                  y = c(-100,-90, -80, -70))
loc2 <- data.frame(x = c(-100,-90, -80, -70),
                   y = c(100,90, 80, 70))
loc3 <- data.frame(x = c(100,90, 80, 70),
                   y = c(100,90, 80, 70))
coordinates(loc1) <- ~x+y
sp.loc1<- SpatialPoints(loc1,proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs"))
coordinates(loc2) <- ~x+y
sp.loc2<- SpatialPoints(loc2,proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs"))
coordinates(loc3) <- ~x+y
sp.loc3<- SpatialPoints(loc3,proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs"))
plot(sp.loc1)

#create rasters which have a common part of the naming convention of point shapefiles targeting for attribution
loc1_blablabla<- raster(xmn=-200, xmx=0, ymn=-200, ymx=0)
loc1_blablabla
ncell(loc1_blablabla)
#it has 64800 cells
values(loc1_blablabla)<-1:64800
plot(loc1_blablabla, add=TRUE)
plot(sp.loc1, add=TRUE)
loc2_blablabla<- raster(xmn=-200, xmx=0, ymn=0, ymx=200)
values(loc2_blablabla)<-1:64800
loc3_blablabla<- raster(xmn=0, xmx=200, ymn=0, ymx=200)
values(loc3_blablabla)<-1:64800

#plot all, but first create extent poly
#borrowed some code from here:https://gis.stackexchange.com/questions/206929/r-create-a-boundingbox-convert-to-polygon-class-and-plot/206952
library(rgeos)
ext.box<-rgeos::bbox2SP(n=250, s=-250, w=-250, e=250)
plot(ext.box)
plot(loc1_blablabla, add=TRUE)
plot(loc2_blablabla, add=TRUE)
plot(loc3_blablabla, add=TRUE)
plot(sp.loc1, add=TRUE)
plot(sp.loc2, add=TRUE)
plot(sp.loc3, add=TRUE)

#now attempt to extract raster values at points for multiple non-overlapping point and raster files ie. extract(loc1_blablabla, loc1)
#try lapply as in: https://stackoverflow.com/questions/59164538/create-a-loop-to-extract-data-from-multiple-raster
#create list as below, however, with my real data I would use list.files()
loc.list<-list(loc1,loc2,loc3)
rast.list<-list(loc1_blablabla,loc2_blablabla,loc3_blablabla)
attr.data<-lapply(rastlist.extract,loc.list)
#this doesn't work -- i think I need coordinates as the last term, not a list of spatial points

#also tried a for loop, but this gave error:  "Error in (function (classes, fdef, mtable): unable to find an inherited method for function ‘shapefile’ for signature ‘"list"’"
for (i in 1:length(loc.list)) {
  #this reads in each spatialpoint file and assigns each sp's name to the variable 'sp.name'
  sp.name<-shapefile(loc.list[i]) 
  attr.data<-data.frame(coordinates(sp.name),
                             extract(rast.list[grep(sp.name,rast.list)],sp.name))
  #eventually add this line to affix the new vector attr.data to the coordinates for each plot:  names(attr.data)<-c("x","y","raster.value")
}

1 Answer 1

1

Using data from your example (thanks for providing such a good example), the error is in the use of lapply. Here the simplest solution is to use mapply instead (or Map) like this:

attr.data<-mapply(raster::extract,rast.list,loc.list)

If you want to speed up computation (if your lists are very long) you could try the extract function from the terra package, which is the new version of the raster package (see documentation).

However, another possible option, as you said, is to merge the raster together and the points to do the extraction in one call with one of the functions that I mentioned above. To merge raster with different extent you can use rgdal and gdalUtils packages:

mergevrt <- function(listname,output.vrt,output.tif,sep=F){
#listname=character vector of raster path to merge
#output.vrt= character. Path of the output virtual raster
#output.tif= Character. Path of the output .tif
gdalbuildvrt(gdalfile=listname ,output.vrt=output.vrt,separate = sep,overwrite = T)
print("vrt done")
gdal_translate(output,dst_dataset = output.tif)
print("tif saved")
}

r <- raster(output.tif)
loc_merge <- rbind(loc1, loc2,loc3)
attr.data<-extract(r,loc_merge)

in terra you have to use the native class of the package, by converting the input like this:

library(terra)
t <- rast(r)
loc_terra <- vect(loc_merged)
attr.data<-terra::extract(r,loc_merge)
3
  • I have found that "on-the-fly" coercion works quite well with terra objects and does not slow things down much so, you can simply do: terra::extract(rast(r), vect(p)) Commented Apr 16, 2021 at 19:58
  • If you build a vrt, then you should be able to just use that, there is no need to copy all the data into a new tif. You can (maybe, not much tested yet) create a vrt with terra::vrt. Also, combining the rasters may not be the way to go if the rasters overlap or if they not have the same origin --- it would change the results. Commented Apr 16, 2021 at 20:52
  • Yes, just the vrt should be enough. I think that combine raster in this particular case should work well because of the non overlapping raster of the example
    – Elia
    Commented Apr 16, 2021 at 21:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.