[R-sig-Geo] Faster way to extract raster values using multiple polygons?
Dear all, I am trying to extract a time series from a raster brick using multiple polygons. The raster brick is a global temperature netcdf file with 1200 layers (months) and the polygons are each of the 5567 municipalities in Brazil. I intend to extract the temperature time series for each municipality. As a final result, I would like to have a data frame containing: -the name and coordinates of the municipality, -the date tag from the raster, and -the extracted temperature value. My initial, VERY SLOW attempt was something like this: library(raster) # Create some sample raster data idx <- seq(as.Date("1960/1/1"), as.Date("1990/12/01"), by = "month") r <- raster(ncol=360, nrow=180) b <- brick(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r) b <- setZ(b, idx) # Import polygon data br <- getData("GADM", country="Brazil", level=2) # about 7MB to download # Subset the SMALLEST state - only 75 municipalities br_sub <- br[br$NAME_1=="Sergipe" & !is.na(br$NAME_1),] plot(br_sub) # Now let's extract the data for each municipality beginCluster() e <- extract(b, br_sub, fun=mean, df=T) endCluster() Even using the smallest state possible, this example takes about 20 minutes to run on a dual-core 2.5GHz Macbook Pro using four threads. As a reference, there are states with over 850 municipalities. And remember, in total there are over 5500 municipalities I need to extract the data for. I have the feeling that my code is not very optimized. Has anybody ever dealt with this amount of data in this kind of raster operation? Is there any fastest way to achieve my expected result? Thanks in advance, -- Thiago V. dos Santos PhD student Land and Atmospheric Science University of Minnesota ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Faster way to extract raster values using multiple polygons?
Hi, The resolution of your raster data (1 degree) is much more coarse than what your polygons represent. Could you short-circuit the process by assuming that the temp at the centroid of each polygon would suitably represent the mean temperature across each polygon? Unless you have some much bigger polygons, I can't imagine it will be very far off. If so, then you could pretty quickly extract the values for each layer in the raster at each centroid. Perhaps like this? cents <- coordinates(br_sub) v <- extract(b, cents) Is that close enough? Cheers, Ben > On Apr 6, 2016, at 8:00 AM, Thiago V. dos Santos via R-sig-Geo > wrote: > > Dear all, > > I am trying to extract a time series from a raster brick using multiple > polygons. > > The raster brick is a global temperature netcdf file with 1200 layers > (months) and the polygons are each of the 5567 municipalities in Brazil. I > intend to extract the temperature time series for each municipality. > > As a final result, I would like to have a data frame containing: -the name > and coordinates of the municipality, -the date tag from the raster, and -the > extracted temperature value. > > My initial, VERY SLOW attempt was something like this: > > > library(raster) > > # Create some sample raster data > idx <- seq(as.Date("1960/1/1"), as.Date("1990/12/01"), by = "month") > r <- raster(ncol=360, nrow=180) > b <- brick(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r) > b <- setZ(b, idx) > > # Import polygon data > br <- getData("GADM", country="Brazil", level=2) # about 7MB to download > > # Subset the SMALLEST state - only 75 municipalities > br_sub <- br[br$NAME_1=="Sergipe" & !is.na(br$NAME_1),] > plot(br_sub) > > # Now let's extract the data for each municipality > beginCluster() > e <- extract(b, br_sub, fun=mean, df=T) > endCluster() > > > Even using the smallest state possible, this example takes about 20 minutes > to run on a dual-core 2.5GHz Macbook Pro using four threads. As a reference, > there are states with over 850 municipalities. And remember, in total there > are over 5500 municipalities I need to extract the data for. > > I have the feeling that my code is not very optimized. > > Has anybody ever dealt with this amount of data in this kind of raster > operation? Is there any fastest way to achieve my expected result? > Thanks in advance, > -- Thiago V. dos Santos > > PhD student > Land and Atmospheric Science > > University of Minnesota > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo Ben Tupper Bigelow Laboratory for Ocean Sciences 60 Bigelow Drive, P.O. Box 380 East Boothbay, Maine 04544 http://www.bigelow.org ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Faster way to extract raster values using multiple polygons?
Hi Ben, Thank you so much for the suggestion. It is pretty fast, and I guess it's a good assumption for temperature. I'll have to do the same with precipitation, which is way more spatially variable than temperature. In that case, I think I'll have to wait the 10 hours it takes to run "extract" using a mean value of all cells covering a polygon. It is a one-time processing anyway, so I should not worry too much about the waiting time. Greetings, -- Thiago V. dos Santos PhD student Land and Atmospheric Science University of Minnesota On Wednesday, April 6, 2016 10:48 AM, Ben Tupper wrote: Hi, The resolution of your raster data (1 degree) is much more coarse than what your polygons represent. Could you short-circuit the process by assuming that the temp at the centroid of each polygon would suitably represent the mean temperature across each polygon? Unless you have some much bigger polygons, I can't imagine it will be very far off. If so, then you could pretty quickly extract the values for each layer in the raster at each centroid. Perhaps like this? cents <- coordinates(br_sub) v <- extract(b, cents) Is that close enough? Cheers, Ben > On Apr 6, 2016, at 8:00 AM, Thiago V. dos Santos via R-sig-Geo > wrote: > > Dear all, > > I am trying to extract a time series from a raster brick using multiple > polygons. > > The raster brick is a global temperature netcdf file with 1200 layers > (months) and the polygons are each of the 5567 municipalities in Brazil. I > intend to extract the temperature time series for each municipality. > > As a final result, I would like to have a data frame containing: -the name > and coordinates of the municipality, -the date tag from the raster, and -the > extracted temperature value. > > My initial, VERY SLOW attempt was something like this: > > > library(raster) > > # Create some sample raster data > idx <- seq(as.Date("1960/1/1"), as.Date("1990/12/01"), by = "month") > r <- raster(ncol=360, nrow=180) > b <- brick(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r) > b <- setZ(b, idx) > > # Import polygon data > br <- getData("GADM", country="Brazil", level=2) # about 7MB to download > > # Subset the SMALLEST state - only 75 municipalities > br_sub <- br[br$NAME_1=="Sergipe" & !is.na(br$NAME_1),] > plot(br_sub) > > # Now let's extract the data for each municipality > beginCluster() > e <- extract(b, br_sub, fun=mean, df=T) > endCluster() > > > Even using the smallest state possible, this example takes about 20 minutes > to run on a dual-core 2.5GHz Macbook Pro using four threads. As a reference, > there are states with over 850 municipalities. And remember, in total there > are over 5500 municipalities I need to extract the data for. > > I have the feeling that my code is not very optimized. > > Has anybody ever dealt with this amount of data in this kind of raster > operation? Is there any fastest way to achieve my expected result? > Thanks in advance, > -- Thiago V. dos Santos > > PhD student > Land and Atmospheric Science > > University of Minnesota > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo Ben Tupper Bigelow Laboratory for Ocean Sciences 60 Bigelow Drive, P.O. Box 380 East Boothbay, Maine 04544 http://www.bigelow.org ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo