Join SpatialPointsDataFrame and SpatialLinesDataFrame using over(), R

Question

I am stuck with something that it seems should be quite straightforward. Apologies, I am new to using spatial data in R.

I am trying to map city data, onto a map of the world's coastlines. I have taken the coastlines from the natural earth data set (https://www.naturalearthdata.com/downloads/) 1:110m data and generated the spatial lines dataframe:

coast_rough_sldf
class       : SpatialLinesDataFrame 
features    : 134 
extent      : -180, 180, -85.60904, 83.64513  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 3
names       : scalerank, featurecla, min_zoom 
min values  :         0,  Coastline,      0.0 
max values  :         1,    Country,      1.5

I further have a dataset of cities, a sample of which looks as follows:

city_coast <- data.frame(Latitude = c(-34.60842, -34.47083, -34.55848, -34.76200, -34.79658, -34.66850), 
              Longitude = c(-58.37316, -58.52861, -58.73540, -58.21130, -58.27601, -58.72825), 
              Name1 = c("Buenos Aires", "San Isidro", "San Miguel", "Berazategui", "Florencio Varela", "Merlo"), 
              distance = c(7970.091,  5313.518, 26156.700, 11670.274, 18409.738, 33880.259))
city_coast

Latitude Longitude            Name1  distance
1 -34.60842 -58.37316     Buenos Aires  7970.091
2 -34.47083 -58.52861       San Isidro  5313.518
3 -34.55848 -58.73540       San Miguel 26156.700
4 -34.76200 -58.21130      Berazategui 11670.274
5 -34.79658 -58.27601 Florencio Varela 18409.738
6 -34.66850 -58.72825            Merlo 33880.259

I then successfully create the spatial points dataframe:

city_spdf <- SpatialPointsDataFrame(coords = select(city_coast, c("Longitude", "Latitude")),
                                    proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84"),
                                    data = select(city_coast, c("Name1", "distance")))

city_spdf

class       : SpatialPointsDataFrame 
features    : 6 
extent      : -58.7354, -58.2113, -34.79658, -34.47083  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 2
names       :       Name1,  distance 
min values  : Berazategui,  5313.518 
max values  :  San Miguel, 33880.259

Now i want to join the city_spdf with the coast_sldf, so that i can plot them using tmap. Looking at tutorials it seems that i should use over():

city_coast_shp <- over(coast_rough_sldf, city_spdf)

city_coast_shp

Name1 distance
1  <NA>       NA

Which is clearly wrong. Switching the order of the objects changes things but still doesn't give me what i need.

Can anyone tell me what i am not getting right with this over function? Every example i have seen simply has people joining the two spatial objects. Apologies if i am missing something extremely simple.

would a similar solution using `library(sf)` be valuable for you? — Nate, Nov 30 '18 at 18:37
Are your points and lines actually overlapping? While Buenos Aires is a coastal city it doesn't mean the point would intersect with your line feature. Can you get a polygon of the city boundaries and try that? Or maybe do something like a nearest neighbor to assign a coast — elmuertefurioso, Nov 30 '18 at 18:46
@elmuertefurioso sorry, perhaps i have misunderstood something. I don't need them to intersect. I simply want to plot them and color them based on a categorical variable i generate from the distance measure. My understanding was that i needed them to be in the same spatial object to use in tmap. Do i have this wrong? — MorrisseyJ, Nov 30 '18 at 18:53
I just don't know `spatial` as well as I do `sf` but if no one gets back tonight, I will help tomorrow with and `sf` solution — Nate, Nov 30 '18 at 21:00
@Nate i tried using sf() to convert the spatiallinesdataframe to a spatialpolygonsdataframe. I followed the instructions from here: https://stackoverflow.com/questions/47147242/convert-spatial-lines-to-spatial-polygons, but generated an error: Error in SpatialPolygonsDataFrame(x, y, match.ID = match.ID, ...) : Object length mismatch. Since i am not really sure what an sf object is, i am struggling to fix this. — MorrisseyJ, Nov 30 '18 at 21:03

score 2 · Accepted Answer · answered Dec 01 '18 at 15:52

Like @elmuertefurioso point out in the comments, I think a one reason this isn't working how you expect is because of confusion of types of geometries.

Since the coastline data is lines, and not polygons like data(World) from tmap, you are restricted a bit in the calculations and comparisons you can make with cities, which is points.

Reading in the data the sf way:

library(sf)

# downloaded from https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/physical/ne_110m_coastline.zip
coastline <- read_sf("~/Downloads/ne_110m_coastline/ne_110m_coastline.shp")

cities <- data.frame(
  Latitude = c(-34.60842, -34.47083, -34.55848, -34.76200, -34.79658, -34.66850), 
  Longitude = c(-58.37316, -58.52861, -58.73540, -58.21130, -58.27601, -58.72825), 
  Name1 = c("Buenos Aires", "San Isidro", "San Miguel", "Berazategui", "Florencio Varela", "Merlo"), 
  distance = c(7970.091,  5313.518, 26156.700, 11670.274, 18409.738, 33880.259)
  )

In order to do any comparisons between sf objects they must have the same Coordinate Reference System. So as we read in cities we will set the CRS to be that of coastline.

cities <- st_as_sf(
  cities,
  coords = c("Longitude", "Latitude"), # must be x, y order
  crs = st_crs(coastline) # must be equivilant between objects
  )

Now you can make comparisons using the st_{comparison}() family of functions.

The function over() and its sf counterpart st_intersects() would work on a set of points and polygons, but we don't have that here. We can use distance functions like st_nearest_feature() with points and lines, to get the closest geometry from coastline for each city.

st_nearest_feature(cities, coastline)

It returns the row index for the nearest geometry in coastlines which happens to be the same for all the cities here because they are all in Argentina. The order matters in the function because it defines the question being asked If we flipped it to st_nearest_feature(coastline, cities)it would return the closest city for each geometry in coastline, so the return would have 134 elements.

All that to say you don't actually have to do any joining or comparisons to plot your points together on the same tmap.

library(tmap)

tmap_mode("view")

tm_shape(coastline) +
  tm_lines() +
tm_shape(cities) +
  tm_bubbles("distance")

I'm not a tmap user but I just zoomed in snapped this screen shot to show its working.

thanks so much for this. It's an extremely thorough answer. I was mis-understanding the operation of `over()`, and thanks also for the `sf()` syntax. The explanation of line and point spatial objects makes sense. Further for `tmap()`, i had things wrong. I was trying: `tm_shape(coastline) + tm_lines() + tm_bubbles(data = cities, size = "distance")`. Following the ggplot() syntax. Apologies for the simple error - glad to have this all working now. Will do more exploring of `sf()`. — MorrisseyJ, Dec 03 '18 at 15:37

Join SpatialPointsDataFrame and SpatialLinesDataFrame using over(), R

1 Answers1