Shapefile – Extracting NUTS-3 Shapefile of Germany from Eurostat’s Shapefile in R

europerrgdalshapefile

I have downloaded shapefile from Eurostat website (download site for downloading Eurostat shapefile and direct download link of Eurostat shapefile that I am interested in). To be short, I need shapefile of Germany' NUTS3 region only (general concept for NUTS3). To do so, I used rgdal package. Here is the R script down below:

library(rgdal)
library(sp)
library(maptools)

url = "http://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/nuts/download/ref-nuts-2013-03m.shp.zip"
download.file(url, basename(url))
gunzip(basename(url))
rname <- list.files(getwd(), "shp$")

but I got following error anyway:

> eu <- readOGR(dsn = getwd(), layer ="NUTS_RG_03M_2013")
Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv = use_iconv,  : 
  Cannot open layer

don't know why does this happen. Any quick solution? How can I fix this error?

If I fix this problem, I want to retrieve shapefile of only Germany's NUTS3 region (please take a look overview of NUTS3). To do so, I come up with following solution:

nuts2 <- subset(eu, STAT_LEVL_ == 2)
proj4string(nuts2)
coordinates(point_data) <- ~LON + LAT
proj4string(point_data) <- CRS("+proj=longlat +datum=WGS84")
point_data <- spTransform(point_data, proj4string(map_nuts2))
over(point_data, nuts2, by.id = TRUE)

My ultimate goal is to get Germany' NUTS3 level shapefile. How can I correct my above solution? How can I refine my above approach? Is there any more efficient and straightforward way to get certain countries' NUTS3 level shapefile? Any idea?

Best Answer

Quick solution? Check if filename exists!

library(rgdal)
library(sp)
library(maptools)

url = "http://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/nuts/download/ref-nuts-2013-03m.shp.zip"
download.file(url, basename(url))
unzip(basename(url))
rname <- list.files(getwd(), "shp$")

# your method
eu <- readOGR(dsn = getwd(), layer ="NUTS_RG_03M_2013")
## Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv = use_iconv,  : 
##                    Cannot open layer

# another method
eu <- shapefile("NUTS_RG_03M_2013.shp")
## Error: file.exists(extension(x, ".shp")) is not TRUE

You have a lot of files with the same pattern name:

list.files(pattern = 'NUTS_RG_03M_2013.*.shp$')
## [1] "NUTS_RG_03M_2013_3035_LEVL_0.shp" "NUTS_RG_03M_2013_3035_LEVL_1.shp"
## [3] "NUTS_RG_03M_2013_3035_LEVL_2.shp" "NUTS_RG_03M_2013_3035_LEVL_3.shp"
## [5] "NUTS_RG_03M_2013_3035.shp"        "NUTS_RG_03M_2013_3857_LEVL_0.shp"
## [7] "NUTS_RG_03M_2013_3857_LEVL_1.shp" "NUTS_RG_03M_2013_3857_LEVL_2.shp"
## [9] "NUTS_RG_03M_2013_3857_LEVL_3.shp" "NUTS_RG_03M_2013_3857.shp"       
## [11] "NUTS_RG_03M_2013_4258_LEVL_0.shp" "NUTS_RG_03M_2013_4258_LEVL_1.shp"
## [13] "NUTS_RG_03M_2013_4258_LEVL_2.shp" "NUTS_RG_03M_2013_4258_LEVL_3.shp"
## [15] "NUTS_RG_03M_2013_4258.shp"        "NUTS_RG_03M_2013_4326_LEVL_0.shp"
## [17] "NUTS_RG_03M_2013_4326_LEVL_1.shp" "NUTS_RG_03M_2013_4326_LEVL_2.shp"
## [19] "NUTS_RG_03M_2013_4326_LEVL_3.shp" "NUTS_RG_03M_2013_4326.shp"

One could be (for WGS 84 - EPSG 4326):

eu <- readOGR(dsn = getwd(), layer ="NUTS_RG_03M_2013_4326_LEVL_2")
Germany <- eu[eu@data$CNTR_CODE == "DE",]
plot(Germany)

Related Solutions

[GIS] Overlaying grids with a polygon to determine the average value of the overlapping grids for that polygon

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
    }

}

[GIS] Need to extract rectangle coordinates from a shapefile

You can create a rectangle with the boundary box of your polygon. I use rgdal package to make a reproducible example, but you could use only raster package with shapefile() function instead of readOGR():

library(rgdal)
library(raster)

dsn <- system.file("vectors", package = "rgdal")[1]
ogrListLayers(dsn)
polygon <- readOGR(dsn=dsn, layer="scot_BNG")

Creates an extent object using bbox and convert it to polygon:

e <- extent(polygon@bbox)

e <- as(e,"SpatialPolygons")

See the results:

plot(polygon)
plot(e, add= T)

Best Answer

Related Solutions

[GIS] Overlaying grids with a polygon to determine the average value of the overlapping grids for that polygon

[GIS] Need to extract rectangle coordinates from a shapefile

Related Question