[R-sig-Geo] Match coordinates with NUTS3 regions

Bede-Fazekas Ákos b|@|ev||@t @end|ng |rom gm@||@com
Thu Apr 4 07:53:34 CEST 2019


Dear Michele,
be sure that the order of the coordinates is  the same for the two objects.
proj4string(EU_NUTS29) is starting with "+proj=longlat" (long, then 
lat), while coordinates(awards_with_coordinates_first2500) is "~g_lat + 
g_lon" (lat, then lon).
HTH,
Ákos Bede-Fazekas
Hungarian Academy of Sciences

2019.04.04. 1:47 keltezéssel, Michele Molteni írta:
> Hi everyone,
>
> I have a .dta file with coordinates (WGS 84) and I would like to match them
> to NUTS3 regions. I followed indications from a previous post but I get no
> match even though the coordinates fall within NUTS3 bboxes. Could you
> please be so kind to help me understand what I am doing wrong? This is what
> I have so far:
>
>> library(rgdal)
>>
>> ## Download Administrative Level data from EuroStat
>>
>> temp <- tempfile(fileext = ".zip")
>> download.file("
> https://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/nuts/download/ref-nuts-2013-60m.shp.zip",
> temp)
> trying URL '
> https://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/nuts/download/ref-nuts-2013-60m.shp.zip
> '
> Content type 'application/zip' length 2330060 bytes (2.2 MB)
> downloaded 2.2 MB
>> unzip(temp)
>>
>> ##Read data after having unzipped it
>> EU_NUTS <- readOGR(dsn = "C:/Users/User/Documents/TED_project/Replicating
> paper/Shapefile", layer="NUTS_RG_60M_2013_4326_LEVL_3")
> OGR data source with driver: ESRI Shapefile
> Source: "C:\Users\User\Documents\TED_project\Replicating paper\Shapefile",
> layer: "NUTS_RG_60M_2013_4326_LEVL_3"
> with 1480 features
> It has 5 fields
>> ##Consider only EU29 countries (1361 NUTS3)
>> EU_NUTS29 <- EU_NUTS[which(EU_NUTS using data$CNTR_CODE %in% c('BE', 'BG',
> 'CZ', 'DK', 'DE', 'EE', 'IE', 'EL', 'ES', 'FR', 'HR', 'IT', 'CY', 'LV',
> 'LT', 'LU', 'HU', 'MT', 'NL', 'AT', 'PL', 'PT', 'RO', 'SI', 'SK', 'FI',
> 'SE', 'UK', 'NO')), ]
>> ##Import data contaning lat and lon and consider first 2500 observations
>> awards_with_coordinates_first2500 <- read_dta("~/TED_project/Replicating
> paper/awards_with_coordinates_first2500.dta")
>> awards_with_coordinates_first2500 <- awards_with_coordinates_first2
> 500[(1:2500),]
>> ## Here is the projection for map_nuts data
>> proj4string(EU_NUTS29)
> [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
>> ##Transform g_lat and g_lon in numeric
>> awards_with_coordinates_first2500$g_lat <- as.numeric(as.character(awards
> _with_coordinates_first2500$g_lat))
>> awards_with_coordinates_first2500$g_lon <- as.numeric(as.character(awards
> _with_coordinates_first2500$g_lon))
>> ## Make a spatial points dataframe from point_data
>> ## Assuming the projection to be longlat on the WGS84 datum
>> ## First define the coordinates
>> coordinates(awards_with_coordinates_first2500) <- ~g_lat + g_lon
>>
>> ## Then the projection
>> proj4string(awards_with_coordinates_first2500) <- CRS("+proj=longlat
> +datum=WGS84")
>> ## Now transform the projection of awards_with_coordinates_first2500 to
> that of EU_NUTS29
>> awards_with_coordinates_first2500 <- spTransform(awards_with_coordinates_first2500,
> proj4string(EU_NUTS29))
>> ## Use the over() function to get the regions for the points
>> over(awards_with_coordinates_first2500,EU_NUTS29,by.id=)
>      NUTS_ID LEVL_CODE CNTR_CODE NUTS_NAME  FID
> 1      <NA>        NA      <NA>      <NA> <NA>
> 2      <NA>        NA      <NA>      <NA> <NA>
>
> 200    <NA>        NA      <NA>      <NA> <NA>
>   [ reached 'max' / getOption("max.print") -- omitted 2300 rows ]
>> bbox(awards_with_coordinates_first2500)
>              min      max
> g_lat 34.794891 51.46780
> g_lon -3.018564 33.63785
>> bbox(EU_NUTS29)
>        min    max
> x -61.806 55.826
> y -21.350 71.127
>> ##The points seem to be all in Africa
>> plot(EU_NUTS29)
>> plot(awards_with_coordinates_first2500, add = TRUE, col = "red", pch = 19)
> Thanks in advance.
>
> Cheers,
> Michele.
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>



More information about the R-sig-Geo mailing list