0

I have a raster stack that I want to extract from:

raster resolution: 0.0417, 0.0417  (x, y). 

My raster is coarser than my coordinates for extraction. I have finer coordinates whose resolution is ~ .000217, .000217. (My coordinates are coarser because I have other rasters with finer resolution). My coordinates aren't quite a rectangular grid (only am using coordinates in the form of a US state) but they increase incrementally in both LAT/LONG direction.

I've been using fast_extract from the library prioritizr & velox after converting my coordinates into SpatialPoints but I have too many coordinates and it's taking too long to go through each coordinate and extract from the raster.

I feel like the most effective way would be to extract at the resolution of the raster, and use that value for all MY coordinates that fall within that raster's pixel.

for example, let's take one pixel, which is from -74.0417 to -74.0834 (LONG) and 40 to 40.0417 (LAT). The value of this pixel would be for all my coordinates that fall in this range, ie. 40.00027 x 74.0417 / 40.00054 x 74.0417 / and so on.

I think I might want to use something like the aggregate function but I don't want to lose the number of my LatLon coordinates by making it coarser---I ultimately want to keep the same number of LatLon coordinates but reduce the number of extractions.

Here is some sample data

LatLon coordinates to extract:

temp_LAT<-seq(44.86,45.10,.0027) #89
temp_LONG<-seq(-80,-71.5,.0027) #3149
temp_LL<-data.frame(matrix(0,nrow = 280261, ncol =2))
colnames(temp_LL)<-c("LAT","LONG")
temp_LL[,"LAT"]<-rep(temp_LAT,each=length(temp_LONG))
temp_LL[,"LONG"]<-temp_LONG
pts<-SpatialPoints(temp_LL[,c("LONG","LAT")])
crs_prj<-"+proj=utm +zone=48 +datum=WGS84"
proj4string(pts)<-CRS(crs_prj)

raster to extract from:

library(raster)
x<-raster()
x <- raster(ncol=35, nrow=150, xmn= - 81.5, xmx= -71, ymn= 40, ymx=45)
res(x)<-.03333
res(x)
projection(x) <- "+proj=utm +zone=48 +datum=WGS84"
values(x) <- 1:ncell(x)
  • 1
    It sounds like you want to call cellFromXY(raster, xy) to get the cell number for every point, then extract raster values at the unique set of cell numbers, then assign the values back to the points by cell number. If you can set up some sample data for a test case I'll answer with full code. – Spacedman Dec 10 '19 at 6:29
  • Hi @Spacedman Thanks for your help, esp with following up with the code. I've edited my code to create a raster + LatLong points (sorry for the delay--learning how to make a fake raster). What you're describing sounds like what I want--do you have suggestions on how to quickly assign values by cell number rather than row by row? With all the points I have, I'm a little concerned about the efficiency of that part. Thanks! – Tammy Dec 10 '19 at 16:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.