[R-sig-Geo] Match Coordinates to NUTS 2 ID
Miluji Sb
milujisb at gmail.com
Fri May 27 10:42:15 CEST 2016
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
>
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list