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
Match Coordinates to NUTS 2 ID
3 messages · Miluji Sb, Frede Aakmann Tøgersen
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 _______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Dear Dr. T?gersen, Thank you very much for your reply. I really appreciate it - I am following your suggestions but I get "NA" after the match here: over(point_data, map_nuts2, by.id = TRUE) I am trying to find what the problem is at my end. Thanks again! Sincerely, Milu On Fri, May 27, 2016 at 7:39 AM, Frede Aakmann T?gersen <frtog at vestas.com> wrote:
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