2

I have the following snippet of code:

tracts<-readOGR('data/.')
proj4string(tracts)
"+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"

tracts_projected <- spTransform(tracts, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
proj4string(tracts_projected)
"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"

If I take the min area of both, the values are different, even though both spatialpolygondataframes are the same except for the difference in projection.

min(area(tracts)/1000000)
 1.221094e-05

min(area(tracts_projected)/1000000)
9.505325e-06

Why is this so? And more importantly which is the better way to get area

  • 1
    Mercator "+proj=merc" is no good for area, it's designed for lines of constant bearing and not area calculations, I would compare the result from tracts_projected (which is calculated using an approximation in longitude/latitude) with area(spTransform(tracts, "+proj=laea +datum=WGS84")) (an equal area projection)- but, it might be better to transform to a local centre - so please print(tracts) and show us the output, and also plot(tracts); axis(1); axis(2) to show us how the data looks. – mdsumner Aug 9 at 8:21

Your Answer

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

Browse other questions tagged or ask your own question.