0

I have a set of coordinates that most fall within or near the Iberian Peninsula.

I'm interested in finding out the country name of a large set of coordinates and the province name. However, there are some problematic coordinates on the boundaries of the countries (coordinates near the coastal area) that are no detected within the country or province. First I have tried to find out the province name and after the country name and I find this issue with these two different approaches.

I'm aware that some coordinates could be typos but I have manually checked that some fall within the country limits.

Any advice of how to overcome this?

Here create here a minimal reproducible example:

First, I create a dataframe with the problematic coordinates. The last coordinate (36.76353, -4.425162) is a "correct one" in order to show that the function works.

lat <- c(36.81973, 43.69739, 36.51103, 38.50000, 42.25167, 42.25177, 39.31316, 36.76353)
long <- c(-2.411557, -5.919138, -4.635675, -0.100000, -8.804174, -8.790300, 2.995276, -4.425162)
coords <-data.frame(long, lat)

Then by readapting the answer of a colleague with the addition of library(mapSpain) Latitude Longitude Coordinates to State Code in R I create this function that should tell the provinces where the different points are.

#I create a function that should tell the provinces of Spain from coordinates
library(sf)
library(spData)
library(mapSpain)

## pointsDF: A data.frame whose first column contains longitudes and
##           whose second column contains latitudes.
##
## states:   An sf MULTIPOLYGON object with 50 states plus DC.
##
## name_col: Name of a column in `states` that supplies the states'
##           names.
lonlat_to_state <- function(pointsDF,
                            states = mapSpain::esp_get_prov(),
                            name_col = "ine.prov.name") {
  ## Convert points data.frame to an sf POINTS object
  pts <- st_as_sf(pointsDF, coords = 1:2, crs = 4326)
  
  ## Transform spatial data to some planar coordinate system
  ## (e.g. Web Mercator) as required for geometric operations
  states <- st_transform(states, crs = 3857)
  pts <- st_transform(pts, crs = 3857)
  
  ## Find names of state (if any) intersected by each point
  state_names <- states[[name_col]]
  ii <- as.integer(st_intersects(pts, states))
  state_names[ii]
}

lonlat_to_state(coords)

[1] NA       NA       NA       NA       NA       NA       NA       "Málaga"

And again I have the same issue if I use

library(maps)
map.where(x = coords$long, y = coords$lat)
[1] NA      NA      NA      NA      NA      NA      NA      "Spain"

Any advice would be more than welcome, thanks for your time!

J. Lan
  • 51
  • 4
  • 2
    maybe this: https://stackoverflow.com/a/39009509/2761575 – dww Feb 15 '22 at 21:35
  • I have tried with countriesHigh(resolution="high") that uses the package rworldextra but still no luck (I get the same output). – J. Lan Feb 15 '22 at 22:13

1 Answers1

1

st_intersects is working as expected, because the points (except one) lie outside the polygons. We can see this by plotting them (zooming in on each)

library(ggplot2)
for (i in seq_len(nrow(pts))) {
  xlim = unlist(st_geometry(pts[i,]))[1] + c(-5000, 5000)
  ylim = unlist(st_geometry(pts[i,]))[2] + c(-5000, 5000)
  g = ggplot() +
    geom_sf(data=states) +
    geom_sf(data=pts) +
    coord_sf(xlim = xlim, ylim = ylim)
  print(g)
}

A few examples of the plots:

enter image description here enter image description here enter image description here

dww
  • 30,425
  • 5
  • 68
  • 111
  • Thanks! I think that more than one coordinate falls within the Iberian peninsula. Still think that I need a high resolution map to solve this. I have been playing around with this https://www.gps-coordinates.net/ and I can see that some are on land. In fact I have thousands of points like this and is over my capacity to check them all. I tried to exemplify this with a minimal example and I do apologize if this question wasn't complete. – J. Lan Feb 15 '22 at 23:03