Extracting point data from raster files

Question

I have many raster data in tiff format. I need to extract the values for some coordinates and save it as a csv file. I have a code to extract data from some stations across the raster flies, and the code works as long as I extract just one station. Once, I put other coordinates to extract their values it gives errors. so I have to run the code for each of stations separately (I have 8 stations).

Can you help me modify the code?

library(raster)
library(sp) # used to create a SpatialPoint object
setwd("c:/PET")


# LOAD RASTERS INTO A LIST OBJECT. 
tmp <- lapply(list.files("c:/PET", pattern = ".tif$", full.names = TRUE), raster)


coords <- data.frame(

# **1st station**
"lat" = c(498957.72),"lon" = c(4516020.57)


)

pts <- SpatialPoints(coords = coords, 
              proj4string = CRS("+proj=utm +zone=16  +datum=WGS84 +units=m"))


outt=as.data.frame(cbind(coords, do.call("rbind", lapply(tmp, extract, pts))))

write.csv(outt, file="myFileName.csv", row.names = FALSE)

Answer 1

I think your penultimate line of code needs modification. Right now, the # objects in your list determines the # rows in your output data.frame. You are then receiving an error because the station coordinates most likely have a different # rows.

In order to fix this, replace your rbind() call inside do.call() with a cbind(), and add the coordinates afterwards. Like that, you make sure that the extracted value matrix and station coordinates have identical # rows. Optionally, you may then transform the output data.frame from wide (ie. high # columns) into long format (ie. high # rows).

library(remote)
library(reshape2)

## example raster data
data(vdendool)
tmp = unstack(vdendool)

## create random points
set.seed(123)
pts = as(sampleRandom(vdendool, size = 5, sp = TRUE), "SpatialPoints")

coords = coordinates(pts)


### option #1: extract using *apply loop ----
xtr = do.call("cbind", lapply(tmp, extract, pts))
colnames(xtr) = names(vdendool)

out = data.frame(coords, xtr)

## if required, transform 'data.frame' from wide into long format
out = melt(out, id.vars = 1:2, variable.name = "layer")

head(out, 10)
#   x    y   layer      value
# 1  -175 67.5 layer.1 -24.620316
# 2  -165 32.5 layer.1  50.426976
# 3    85 62.5 layer.1 -34.658461
# 4   -45 27.5 layer.1 -17.551138
# 5  -125 22.5 layer.1 -12.083841
# 6  -175 67.5 layer.2  40.380589
# 7  -165 32.5 layer.2  25.252061
# 8    85 62.5 layer.2   4.944444
# 9   -45 27.5 layer.2   2.678345
# 10 -125 22.5 layer.2  -9.226877

Please note that in case your Raster* objects have the same properties (CRS, extent, resolution), you may as well stack() the single layers currently stored in the list and run extract() directly on the resulting RasterStack, thus avoiding the need for a loop structure altogether.

### option #2: use stacked Raster* objects -----

rst = stack(tmp)
xtr = extract(rst, pts)
out = data.frame(coords, xtr)

## same here, use reshape2::melt to bring 'data.frame' into long format