[R-sig-Geo] readOGR default unit?

[Roger Bivand]

On Fri, 2 Mar 2012, Christian Jansson wrote:

> Hi,
>
> I read some MapInfo TAB-files (or MID MIF, same result) like this:
> shp_0 <- rgdal::readOGR(dsn=fnP, layer=xR)
>
> When I look at summary(shp_0) I get:
> -------------------------------------------
> Is projected: FALSE
> proj4string :
> [+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0]
> -------------------------------------------
>
> Note: "No units="!
>
> The coordinates I read are lat long and show up correctly when I plot them
> on maps. But the problem is I need to use the Area-info for each polygon,
> and the values associated to each polygon are extremely low, like
> 0.000323445. The largest area-number I've ever seen is 0.04.. The polygons
> should be in the tens or hundreds of square kilometers.

If the representation of the coordinate reference system in the source 
file did not include units, they will not be present in the object read. 
The units here are - of course - decimal degrees, so why you would expect 
measurements of area in a metric is a misunderstanding on your part. You 
can try areaPolygon() in geosphere:

library(maptools)
xx <- readShapeSpatial(system.file("shapes/sids.shp",
  package="maptools")[1], IDvar="FIPSNO",
  proj4string=CRS("+proj=longlat +ellps=clrk66"))
library(geosphere)
res <- lapply(slot(xx, "polygons"),
  function(x) sapply(slot(x, "Polygons"),
  function(y) areaPolygon(slot(y, "coords"), r=6378.137)))
summary(unlist(res))
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#  41.78  799.00 1182.00 1179.00 1464.00 2456.00
# in square km for this r.

Or project with spTransform in rgdal:

library(rgdal)
proj4string(xx) <- CRS("+proj=longlat +datum=NAD27")
xxlcc <- spTransform(xx, CRS("+init=epsg:3631"))
res1 <- lapply(slot(xxlcc, "polygons"), function(x)
  sapply(slot(x, "Polygons"), slot, "area")/1000000)
summary(unlist(res1))
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#  41.68  797.20 1179.00 1176.00 1461.00 2450.00


for this datum transformation and projection (and choice of sphere 
parameters in areaPolygon).

Hope this helps,

Roger

>
> I suspect the readOGR is using another default-unit than meters. How can I
> know what it is or how to set it, so that readOGR calculate areas I can use?
>
> By the way, if I look in the Mapinfo-MIF file, it has "CoordSys Earth
> Projection 1,104". Not sure if that is relevant here.
>
> Have a great weekend everyone.
>
> /Chris
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> 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, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no