OSM, rgeos, osmar, area calculation, does not add up
Asked Answered
A

2

7

I am trying to get the size of a polygon from OSM, using osmar to download the data. However, sanity check tells me the are is not right.

Below is an example of what I mean.

(1) Geographical area around Hyde Park in London. Extracting all ways and relations tagged as 'park'.

devtools::install_github('osmdatar/osmdata')
library(osmdata)
library(osmar)
library(sp)
library(sf)
library(rgeos)

osmO <- get_osm(center_bbox(-0.167919, 51.5072682, 2000, 2000))

ids_relations <- osmO$relations$tags[osmO$relations$tags$v=="park","id"]
ids_ways <- osmO$ways$tags[osmO$ways$tags$v=="park","id"]
ids_sub <- find_down(osmO, way(c(ids_relations, ids_ways)))
sp_sub_park <- as_sp(subset(osmO, ids = ids_sub), "polygons")

Now, I want to know the area of each of these 'parks' [Fig1] (the big one in the middle being Hyde Park).

spplot(sp_sub_park, c("version"), colorkey = FALSE, col.regions=c('green'))

Fig1

There are two ways:

1) Use the slot 'area' in the polygon itself.

a1 <- sapply(sp_sub_park@polygons, function(x) x@area)

2) Calculate the area with the specified projection.

bg_poly_t <- spTransform(sp_sub_park, CRS("+proj=longlat +datum=WGS84"))
a2 <- rgeos::gArea(bg_poly_t, byid=TRUE)

These two give me the same result [Fig2] (note the two biggest areas are Hyde Park, split in two by a road)

plot(a1*1000000, a2*1000000)

Fig2

However, the size is not what I would expect. The area is returned in square-km (plotted square-meters). According to that the two parts of Hyde Park add up to about 300 square-meters, the size of a big flat but not a park (Hyde Park ~1.420.000 square-meters).

Any ideas?

Avrilavrit answered 14/6, 2017 at 16:17 Comment(1)
rgeos::gArea gave you a warning, which you could have used as a clue.Three
S
4

Since your "map" is in lat/lon coordinates, to compute the areas you have to either convert it to a "metric" projection (as suggested by @Phil), or use spherical geometry (e.g., as implemented in package geosphere).

Luckily, function st_area of the sf package computes areas of polygons using geosphere in case the object is in geographical coordinates. Therefore, you can do simply:

sf_sub_park <- st_as_sf(sp_sub_park)
areas <- st_area(sf_sub_park)
sum(areas)

, giving:

2551269 m^2

, which is pretty close to @Phil results.

HTH

Subreption answered 14/6, 2017 at 20:40 Comment(1)
not only that, it also gives you the unit of measurement.Three
J
1

I've transformed the data into British National Grid which uses metres as its unit of length (WGS84 maybe uses degrees, IDK?). The area is then 2.5million sq m, which seems more plausible?

sp_sub_park <- spTransform(sp_sub_park, CRS("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs "))
gArea(sp_sub_park)
# [1] 2550387

This is higher than your estimated ~1.4million sq m, but is Hyde Park both large areas of your map (i.e. is Kensington Gardens separate?).

Jardine answered 14/6, 2017 at 17:0 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.