4

I have species distributions (in vector format) and I would like to calculate how much of the area for a given species distribution is within certain lines of latitude, such as temperate and tropic zones.

Using the Wikipedia definition: The north temperate zone extends from the Tropic of Cancer (approximately 23.5° north latitude) to the Arctic Circle (approximately 66.5° north latitude). The south temperate zone extends from the Tropic of Capricorn (approximately 23.5° south latitude) to the Antarctic Circle (at approximately 66.5° south latitude).

Thus, the tropics would be between 23.5° south latitude and 23.5° north latitude.

For example, using this shapefile of the Atlantic Ocean (choose shapefile in the drop-down menue) plotted on a worldmap, one could easily calculate the total area of the Atlantic Ocean;

require(sp)
require(ggplot2)
require(rgeos)
require(rgdal)
require(maps)

setwd("~/test/iho")
ao <- readOGR(getwd(), layer="iho")
aof <- fortify(world, region="name")

# Not necessary for the calculation per se, but still nice. Although not the best looking map
world <- map("worldHires", col="gray90", fill=T)

# Plot 
pp <- ggplot(data = world, aes(x = long, y = lat, group = group)) +
geom_polygon(fill = "grey50") +
geom_polygon(data = AO, aes(x = long, y = lat, group = group), fill = alpha("cyan", 0.5)) +
coord_equal()

# Total area
gArea(SpatialPolygons(ao@polygons))
[1] 7512.821

But how could we restrict the calculation of the area to specified lines of latitude, e.g. tropic and temperate zones as defined above?

Plot

enter image description here

  • In what format are your distributional data? Vector or raster? Counts, densities, ranges, or something else? – whuber Apr 4 '14 at 16:29
  • Thanks for commenting. I have used merged shapefiles, so the Area distributions are in SpatialPolygons class (as the example above). See related thread gis.stackexchange.com/questions/82667/… – jO. Apr 4 '14 at 16:33
  • 1
    Would this be a correct restatement of your need? "Compute the area of a given polygon (in vector format) lying between specified lines of latitude." – whuber Apr 4 '14 at 16:37
  • Cheers- yeah, that sounds better! – jO. Apr 4 '14 at 16:38
4

You can use functions from the rgeos package to extract such regions (e.g. gIntersection, gDifference). I use gDifference in this example, because gIntersection returns a SpatialCollections object here:

# define rectangular region 
y_lim <- c(-1, 1)*23.5
rect_lim <- cbind(c(rep(bbox(ao)["x", ], each=2), bbox(ao)["x", 1]),
                  c(y_lim, rev(y_lim), y_lim[1]))

rect <- SpatialPolygons(list(Polygons(list(Polygon(rect_lim)), ID="1")))
proj4string(rect) <- proj4string(ao)

# compute difference between the two geometries
res <- gDifference(ao, rect)

plot(ao, axes=TRUE)
plot(res, border="red", lwd=2, add=TRUE)
plot(rect, col="#00FF0010", add=TRUE)

plot

# area between 23.5° south latitude and 23.5° north latitude
gArea(SpatialPolygons(ao@polygons)) - gArea(SpatialPolygons(res@polygons))
# [1] 2355.3448

It is important to note that this area is in square degrees. You have to use spTransfrom with an appropriate projection (see @WHuber's comment below). gArea returns also a warning (you didn't fixed it, you used a workaround):

R> gArea(ao)
[1] 7513
Warning message:
In RGEOSMiscFunc(spgeom, byid, "rgeos_area") :
  Spatial object is not projected; GEOS expects planar coordinates
  • 3
    What are the units in which that area is reported? If seems like square degrees, which would be inappropriate for data that extend so far in latitude. It is crucial to perform the computation with data projected in a cylindrical equal area projection. – whuber Apr 5 '14 at 0:16
  • 1
    yes, it's square degrees, I've just reused the code from the question; gArea returns also a warning: gArea(ao): Spatial object is not projected; GEOS expects planar coordinates – rcs Apr 5 '14 at 5:42

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.