1

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.

5
  • would a similar solution using library(sf) be valuable for you?
    – Nate
    Nov 30, 2018 at 18:37
  • 1
    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 Nov 30, 2018 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, 2018 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, 2018 at 21:00
  • @Nate i tried using sf() to convert the spatiallinesdataframe to a spatialpolygonsdataframe. I followed the instructions from here: stackoverflow.com/questions/47147242/…, 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, 2018 at 21:03

1 Answer 1

2

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.

enter image description here

1
  • 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 3, 2018 at 15:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.