1

I would like to add the average, minimum and maximum value of, let’s say soil moisture (variable "value" in my case), to each NUTS3 region (= a NUTS3 region is a polygon of a certain small area in Europe).

However, soil moisture values are stored in little grids. Therefore, I need to find a way to calculate which grids overlap with a certain NUTS3 region in order to determine the average soil moisture value for that NUTS3 region. Visually, it can be compared to this picture where little grids approximate together the shape of US states: http://fieldguide.mt.gov/RangeMaps%5CESMapQUAD_9217.jpg

It is possible that one NUTS3 region consists of 2 full grid cells and 3 half grid cells, while another NUTS3 region consists of only a 0.75 grid. It is also possible that some grids are empty because I don’t have values there.


My raster and NUTS shapefile can be found through the following link: https://drive.google.com/folderview?id=0By9u5m3kxn9yfjZtdFZLcW82SWpzT1VwZXE1a3FtRGtSdEl1c1NvY205TGpack9xSFc2T2s&usp=sharing

  • library(rgdal)
  • grid <- readOGR(".", layer = "Filled_grid") # the variable "value" is the variable that I want to transfer to the respective NUTS3 regions.
  • nuts <- readOGR(".", layer = "NUTS_RG_60M_2006") # in the end, in this shapefile should be an additional variable with the average of "value". This information should come from overlapping the NUTS3 regions with the grids.

I would like to find a solution for this problem in R. For ArcGis, I already found numerous related questions (such as Extract Raster Value into Polygon Attribute ) that can help me.

In R, I found out that the rasterToPolygons command of the raster package, cannot help me because it takes forever due to the amount of data I have.

Another option in R could be this one: https://stat.ethz.ch/pipermail/r-sig-geo/2011-March/011093.html , however, it does not work for me and I am not sure how to indicate that I want the max, min and mean of the variable "value". I believe the mean reason why I have trouble with solutions proposed on the internet is the fact that I stored the raster data in a shapefile, and then want to overlap with another shapefile (with the nuts3 regions). But this shouldn't be a problem, or is it?

I am a beginner in ArcGis but have some basic knowledge in R.

  • If the raster-based calculation is taking forever, then the calculation using a shapefile conversion of the raster will take forever-squared. Go back to the raster and do the conversion in the other way: represent the NUTS regions in raster format. You will have to use a smaller cellsize to succeed. – whuber Jul 24 '15 at 19:00
2

Depending on whether you want a mean or a spatially weighted mean there are a couple things you can do using the rgeos and raster packages

library(rgeos)
proj4string(grid) <- proj4string(nuts) # I assumed these were the same projection?

First find out which grid cells intersect your NUTS polygons

grid_nuts <- gIntersects(grid,nuts,byid = TRUE)

Then you can use the apply() function to calculate the mean, min, and max of your value. However, for the mean this strategy will calculate the mean of any grid cell which touches that NUTS polygon.

nuts@data$average_value <- apply(grid_nuts,1,function(x) mean(grid@data$value[x]))
nuts@data$min_value <- apply(grid_nuts,1,function(x) min(grid@data$value[x]))
nuts@data$max_value <- apply(grid_nuts,1,function(x) max(grid@data$value[x]))

That might work for you, but if you're bothered by the fact that a grid cell that just overlaps NUTS a couple square km is worth the same as one that is totally overlapping, then you need a spatially weighted mean.

Using the intersect() function from the raster package you get the intersection of grid and nuts and it merges the dataframes. (which means the grid cells that are shared between 2 NUTS polygons are split into 2)

library(raster)
new_spdf <- intersect(grid,nuts)

Then as before, find out which of the new grid cells overlap

grid_nuts <- gIntersects(new_spdf,nuts,byid = TRUE)

There's probably a more 'elegant' way to arrange this, but with this you can take all the values in a particular NUTS polygon, multiply by each the areas of the corresponding new grid cells then divide by the total area of the NUTS polygon

EDIT: added the if statement to avoid non-overlapping NUTS polygons

for(i in 1:length(nuts)){
    if(any(grid_nuts[i,])){
        nuts@data$average_spatial_value[i] <- mean(new_spdf@data$value[grid_nuts[i,]]*
                                                       gArea(new_spdf[grid_nuts[i,],])/
                                                       gArea(nuts[i,]))
    } else {
        nuts@data$average_spatial_value[i] <- NA
    }

}
  • Hi @Remi ! Thank you very much for your answer. The first part (the unweighted part) works perfectly for me. However, when trying the second part (the for loop in particular), I got this error message: Error in mean(new_spdf@data$value[grid_nuts[i, ]] * gArea(new_spdf[grid_nuts[i, : error in evaluating the argument 'x' in selecting a method for function 'mean': Error in is.projected(spgeom) : error in evaluating the argument 'obj' in selecting a method for function 'is.projected': Error in if (is.numeric(i) && i < 0) { : missing value where TRUE/FALSE needed – user33125 Jul 30 '15 at 10:13
  • 1
    Sorry about that @user33125, I developed the code on a truncated dataset to save time. In the full dataset you have some cases where nuts and grid do not overlap which was causing an error since you can't take the mean or area of nothing! I edited the above to include an if statement that will assign NA to the average_spatial_value when there is no overlap – Remi Daigle Aug 1 '15 at 3:44
  • No need to apologize! It works now and I can't thank you enough for your help! Thank you for the well structured, patiently explained answer! – user33125 Aug 3 '15 at 7:26

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.