3

Basically, I want to know how many times a 0.25°lat x 0.25°lon patch is fitting in various polygons located around the world. The latter have sizes around 3°lat x 10°lon or 2°lat x 4°lon.

I have the lat/lon values of the corners of the polygons and I am calculating their area like this:

from pyproj import Proj
from shapely.geometry import shape

def getArea(coords):
    c = {"type": "Polygon",
    "coordinates": [[ (coords[0], coords[2]), (coords[1], coords[2]),
                      (coords[0], coords[3]), (coords[1], coords[3]) ]]}
    lon, lat = zip(*c['coordinates'][0])
    pro = Proj("+proj=aea")
    x, y = pro(lon, lat)
    poly = {"type": "Polygon", "coordinates": [zip(x, y)]}
    return shape(cop).area

I took the approach from here: How to calculate the area of a polygon on the earth's surface using python?

Now the question is, which equal area projection shall I choose in order to have comparable area sizes for the polygons. The area of the small patch is always the same, regardless of where it is located on the globe in such a projection.

Taking the Albers Equal Area Projection (aea) results in these areas of three polygons:

  1. 240993868.90978813
  2. 699931593.1047173
  3. 212092562.5238676

Taking the Lambert Azimuthal Equal Area Projection (laea) results in these areas of the same polygons:

  1. 148709452.69292444
  2. 409253749.5468254
  3. 106218747.36092758

Why is the relation between the areas in the two projections different? First 1:3 = 0.344; Second 1:3 = 0.363; They should be the same since both are equal area projections?!

That makes me wonder whether it is legitimate to compare the small patch to the areas of the polygons in either projection. Do you have any advice?

Community
  • 1
  • 1
HyperCube
  • 3,870
  • 9
  • 41
  • 53
  • 1
    2°lat x 3°lon is pretty big, any errors will be proportionately big. And the reason there are different projections is that no projection is perfect, it simply can't be. – Mark Ransom Oct 30 '12 at 21:54
  • I updated the sizes of the patch/polygons just for the record. Hmm, so your saying the error is getting to large with this method? What else can I try? – HyperCube Oct 30 '12 at 22:15
  • 1
    The only way to do it perfectly is to calculate the surface area of the cut-out of a sphere without projecting it. Otherwise you'll have to decide on a projection based on the properties that are important to you. – Mark Ransom Oct 30 '12 at 22:41
  • 2
    "The area of the 2x3 patch is always the same" - not true. Meridianal degree is approximately 111 km as is one degree of latitude near equator. But the latter becomes smaller as the location moves from equator to a pole. – Igor Oct 31 '12 at 13:57
  • 1
    Maybe you can find something useful from [this question](http://gis.stackexchange.com/questions/711/how-can-i-measure-area-from-geographic-coordinates) – nadya Nov 01 '12 at 06:16

1 Answers1

2

If you care about actual relative area, counting tiles won't give you the right answer. If you need actual surface area, use ellipsoidal or spherical geometry, or the idea below.

Brute force vector approach: generate world-wide .25 deg grid as polygons, intersect with polygons, count results. Each polygon can be a Cartesian square or the actual ellipsoidal square. Create a custom ideal projection for each tile, then calculate and store the area as an attribute for each tile. You only have to do this once :)

Chris
  • 21
  • 1