[R-sig-Geo] Match Coordinates to NUTS 2 ID

Fri May 27 07:39:31 CEST 2016

Hi Milu

Your original point data is not within the bounding box of your map_nuts2 data:

> bbox(point_data)
        min     max
LON -125.25 -124.25
LAT   42.75   49.75
> bbox(map_nuts2)
        min      max
x -61.74369 55.83663
y -21.37656 71.17308

So I changed the longitudes of your point data only for illustration. This is one way to do it.

> library(“rgdal”)
Loading required package: sp
rgdal: version: 0.9-2, (SVN revision 526)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
Path to GDAL shared files: c:/Programmer/R/R-3.1.3/library/rgdal/gdal
GDAL does not use iconv for recoding strings.
Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
Path to PROJ.4 shared files: c:/Programmer/R/R-3.1.3/library/rgdal/proj

> ## Download Administrative Level data from EuroStat
> temp <- tempfile(fileext = ".zip")
> download.file("http://ec.europa.eu/eurostat/cache/GISCO/geodatafiles/NUTS_2010_60M_SH.zip",
+               temp)
trying URL 'http://ec.europa.eu/eurostat/cache/GISCO/geodatafiles/NUTS_2010_60M_SH.zip'
Content type 'application/zip' length 6187118 bytes (5.9 MB)
opened URL
downloaded 5.9 MB

> unzip(temp)

> ## Read data
> EU_NUTS <- readOGR(dsn = "./NUTS_2010_60M_SH/data", layer = "NUTS_RG_60M_2010")
OGR data source with driver: ESRI Shapefile 
Source: "./NUTS_2010_60M_SH/data", layer: "NUTS_RG_60M_2010"
with 1920 features
It has 4 fields

> ## Subset NUTS 2 level data
> map_nuts2 <- subset(EU_NUTS, STAT_LEVL_ == 2)

> ## I also have data for a variable by coordinates, which looks like this:
> ## edited the longitudes from original example
> ## Only you know in which projection these coordinates are
> point_data <- structure(list(LON = c(15.25, 16.75, 17.25, 14.25, 14.25, 13.25),
+                              LAT = c(49.75, 49.25, 42.75, 43.25, 48.75, 49.25),
+                              yr = c(2.91457704560515, 9.94774197180345, -2.71956412885765,
+                                  -0.466213169185147, -36.6645659563374, 10.5168056769535)),
+                         .Names = c("LON", "LAT", "yr"), row.names = c(NA, 6L),
+                         class = "data.frame")

> ## Here is the projection for map_nuts data
> proj4string(map_nuts2)
[1] "+proj=longlat +ellps=GRS80 +no_defs"

> ## Make a spatial points dataframe from point_data
> ## Assuming the projection to be longlat on the WGS84 datum
> ## First define the coordinates
> coordinates(point_data) <- ~LON + LAT
> ## Then the projection
> proj4string(point_data) <- CRS("+proj=longlat +datum=WGS84")

> ## Now transform the projection of point_data to that of map_nuts2
> point_data <- spTransform(point_data, proj4string(map_nuts2))

> ## Check the bounding boxes
> bbox(map_nuts2)
        min      max
x -61.74369 55.83663
y -21.37656 71.17308
> bbox(point_data)
      min   max
LON 13.25 17.25
LAT 42.75 49.75

> ## Visualize data
> ## Two points seem to be in the ocean
> plot(map_nuts2)
> plot(point_data, add = TRUE, col = "red", pch = 19)

> ## Use the over() function to get the regions for the points
> ## No regions for 2 points in the ocean
> over(point_data, map_nuts2, by.id = TRUE)
  NUTS_ID STAT_LEVL_ SHAPE_Leng SHAPE_Area
1    CZ06          2   7.018153   1.688086
2    CZ06          2   7.018153   1.688086
3    <NA>         NA         NA         NA
4    <NA>         NA         NA         NA
5    CZ03          2   8.350082   2.222739
6    CZ03          2   8.350082   2.222739
>

Yours sincerely / Med venlig hilsen

Frede Aakmann Tøgersen
Specialist, M.Sc., Ph.D.
Plant Performance & Modeling

Technology & Service Solutions
T +45 9730 5135
M +45 2547 6050
frtog at vestas.com
http://www.vestas.com

Company reg. name: Vestas Wind Systems A/S
This e-mail is subject to our e-mail disclaimer statement.
Please refer to www.vestas.com/legal/notice
If you have received this e-mail in error please contact the sender. 

-----Original Message-----
From: R-sig-Geo [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of Miluji Sb
Sent: 26. maj 2016 23:31
To: R-sig-geo mailing list
Subject: [R-sig-Geo] Match Coordinates to NUTS 2 ID

Dear all,

I have downloaded the NUTS 2 level data from
library(“rgdal”)
library(“RColorBrewer”)
library(“classInt”)
#library(“SmarterPoland”)
library(fields)

# Download Administrative Level data from EuroStat
temp <- tempfile(fileext = ".zip")
download.file("
http://ec.europa.eu/eurostat/cache/GISCO/geodatafiles/NUTS_2010_60M_SH.zip
",
              temp)
unzip(temp)

# Read data
EU_NUTS <- readOGR(dsn = "./NUTS_2010_60M_SH/data", layer =
"NUTS_RG_60M_2010")

# Subset NUTS 2 level data
map_nuts2 <- subset(EU_NUTS, STAT_LEVL_ == 2)

I also have data for a variable by coordinates, which looks like this:

structure(list(LON = c(-125.25, -124.75, -124.25, -124.25, -124.25,
-124.25), LAT = c(49.75, 49.25, 42.75, 43.25, 48.75, 49.25),
    yr = c(2.91457704560515, 9.94774197180345, -2.71956412885765,
    -0.466213169185147, -36.6645659563374, 10.5168056769535)), .Names =
c("LON",
"LAT", "yr"), row.names = c(NA, 6L), class = "data.frame")

I would like to match the coordinates to their corresponding NUTS 2 region
- is this possible? Any help will be high appreciated. Thank you!

Sincerely,

Milu

	[[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo