1

I am trying to calculate the distance from several pixels of a raster to a boundary of some polygons in R as it is showed in the image. I have two problems:

  1. How to correctly define a boundary from a SpatialPolygon object. I currently have two options, but I can not decide for which is the best one. First, I can coerce the SpatialPolygon as a SpatialLines object and then to SpatialPoints. Second, I can raster::rasterize the SpatialPolygon. I suppose each option has its pros and cons.

  2. I only want make the distance calculation from one of the sides of the border. In the example, I only want the distance to the border from the pixels inside Colombia to the border between both countries. I tried to "filter" the border using a latitude parameter but the real border is quite irregular, so using a strict cut-off value can lead to error.

The following code can be used to get the a example shapefile:

#Get the world map and select two countries [Colombia and Venezuela]
library(maptools) #To get the shp data 
data(wrld_simpl)
colven <- c("Colombia", "Venezuela")
colven_map <- wrld_simpl[wrld_simpl@data$NAME %in% colven, ]

Raster data can be created with this code:

raster <- raster(nrow=100, ncol=100, crs=proj4string(wrld_simpl))
raster_extent <- extent(colven_map)
raster <- setExtent(raster, raster_extent)
raster[] <- 1:length(raster)
raster_colven <- mask(raster, coven)
raster_colven <- crop(raster, colven_map)
raster_colven <- mask(raster_colven, colven_map)

Example

2

Perhaps you've already found an answer but this is how I would do it.

  1. I would use rgeos::gIntersection for two adjacent polygons.
  2. I'm assuming you'd like the distance from centroid of each raster cell to the border. I basically use raster::rasterToPoints and rgeos::gDistance to find the distance from the shared border and the raster centroids.

I also use raster::getData to use GADM boundaries instead of maptools.

library(maptools)
library(raster)
library(rgeos)
colombia <- raster::getData(name = "GADM", country = "CO", level = 0)
venezuela <- raster::getData(name = "GADM", country = "VE", level = 0)
venezuela <- sp::spChFIDs(venezuela, "2") # Unique IDs for join
colven <- maptools::spRbind(colombia, venezuela)
border <- rgeos::gIntersection(colombia, venezuela)

Get the raster cell centers by country

raster <- raster(nrow=100, ncol=100, crs=proj4string(colven))
raster <- setExtent(raster, extent(colven))
raster_points <- raster::rasterToPoints(raster, spatial = TRUE)
colombia_points <- raster_points[colombia, ]
venezuela_points <- raster_points[venezuela, ]

Project and calculate distance with rgeos::gDistance. I'm using the equidistant cylindrical projection and kilometers.

crs.edc <- CRS("+proj=eqc +lat_ts=60 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=km +no_defs")

rgeos::gDistance(spTransform(colombia_points, crs.edc),
                 spTransform(border, crs.edc), byid = TRUE)
rgeos::gDistance(spTransform(venezuela_points, crs.edc),
                 spTransform(border, crs.edc), byid = TRUE)

Plot

pch <- "."
plot(colven)
plot(colombia_points, pch = pch, add = TRUE, col = "green")
plot(venezuela_points, pch = pch, add = TRUE, col = "blue")
plot(border, add = TRUE, col = "red")

Hopefully this helps

sample

  • Thank you sboysel, I used a similar approach. Nonetheless, now that you mention it, which command you will use if the polygons were not adjacent, but one on top of the other one? – topcat Feb 16 '16 at 2:16
  • gIntersection would simply return the intersection of two polygons. Depending on your objective, perhaps some use of gDifference will give you a workable border between to polygons. I've never done this personally. – sboysel Feb 16 '16 at 18:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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