[R-sig-Geo] Changing the units of shape file areas

[Roger Bivand]

On Sat, 29 Mar 2014, laanders wrote:

> Hello!
>
> I've been working with Philippine shape files I found at:
> http://www.philgis.org/freegisdata.htm
> Specifically, I've been looking at the barangay level shape files.  I'm
> trying to get the area of each barangay so I can bet the population density
> and plot a heat map.  Unfortunately, I don't understand the units of these
> areas.

Thanks for providing a reproducible example:

library(rgdal)
bg <- readOGR(".", "Barangays")
ZCB <- bg[bg$NAME_2=="Zamboanga City",]
proj4string(ZCB)
plot(ZCB, axes=TRUE)
bbox(ZCB)

The imported data is in geographical coordinates (WGS84 is the datum, it 
is not a projection), so the "areas" are in decimal degrees, which is not 
helpful. To calculate areas, you must project the features, and to do that 
you must find a relevant projection. Using the EPSG table:

EPSG <- make_EPSG()
EPSG[grep("Philippines", EPSG$note),]

and searching around, it appears that 3124 # PRS92 / Philippines zone 4 
may be appropriate - others might also be:

CRS("+init=epsg:3124")
ZCBtm <- spTransform(ZCB, CRS("+init=epsg:3124"))
plot(ZCBtm, axes=TRUE)

with units in metres. To get the areas:

library(rgeos)
AR <- gArea(ZCBtm, byid=TRUE)
summary(AR)

and AR (here in square metres) can be added as a variable into ZCB:

ZCB$area <- AR

if you want to continue using the original unprojected dataset.

If you choose another projection and datum, the areas may vary:

ZCButm <- spTransform(ZCB, CRS("+proj=utm +zone=51 +datum=WGS84"))
summary(abs(AR - gArea(ZCButm, byid=TRUE)))

Hope this helps,

Roger

> Here's the R code I'm using to get the areas: 
>
> ###########################
> BarangaysPlot <- readOGR("~/Barangays", layer = "Barangays”) 
> ZCB <- BarangaysPlot[BarangaysPlot$NAME_2=="Zamboanga City”,]  ## Subsetting
> the shape file to get ZamboCity
>
> sapply(slot(ZCB, "polygons"), slot, "area”)  ## Get the areas
>
> [1] 1.392389e-04 1.946191e-04 7.702360e-05 2.711882e-03 4.610627e-05
> 5.659337e-05
> [7] 6.926905e-05 5.678706e-05 8.203009e-04 1.639304e-03 1.336925e-03
> 4.909287e-03
> ###########################
>
> Since the areas are so small, I thought maybe there were either holes in the
> shape files or I am dealing with the degree coordinate system.  The
> coordinates are planar, so we're good there.  So to see if there were holes
> I used the code provided by Dr. Bivand in an earlier post, I found there
> were no holes.  According to the website, the current projection for these
> shape files is WGS 1984.  So I then looked at spTransform() using the code I
> mentioned earlier by Dr. Bivand.
>
> ###########################
> ## epsg for WGS 1984 is 4326
> ZCBareas <- spTransform(ZCB, CRS("+init=epsg:4326")) 
> ## Gives some “non-simple polygon objects"
> sapply(slot(ZCBareas, "polygons"), function(x) length(slot(x, "Polygons"))) 
>
> ## Gives any holes in the polygons; there are none
> sapply(slot(ZCBareas, "polygons"), function(x) { 
> xi <- slot(x, "Polygons") 
> any(sapply(xi, slot, "hole")) 
> }) 
>
> ## Gives the same areas as before
> sapply(slot(ZCBareas, "polygons"), slot, "area”) 
> ###########################
>
> Again, I'm just trying to find the areas in a unit I can understand to apply
> to population densities. 
> Any help would be much appreciated!! 
> Lew
>
>
>
> --
> View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Changing-the-units-of-shape-file-areas-tp7586083.html
> Sent from the R-sig-geo mailing list archive at Nabble.com.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no