TransWikia.com

Loop to extract specific raster/spatial point pairs

Geographic Information Systems Asked on June 3, 2021

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")
}

One Answer

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)

Answered by Elia on June 3, 2021

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