1

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)
4

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
2
  • Thank you for your answer, I am just not sure how to add my coordinates to your code? what does pts does in your script as I think it's different with mine
    – Seji
    Jul 24 '18 at 16:18
  • Since you didn't provide a reproducible example, I had to select points myself; that's what sampleRandom() is doing. I converted 'pts' to SpatialPoints, which is similar to what you are dealing with when calling SpatiaPoints() on your 'coords' variable. The rest remains the same.
    – fdetsch
    Jul 25 '18 at 5:18

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.