Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I need to calculate how much (in percent) each of thousands of spatial polygons shrink after some events.

I use R and the very useful sp-library.

Example:

  • Projection: "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

  • poly.1 is a Spatial Polygon. It has an area-slot with a value
    area.1

  • poly.2 is the poly.1 after shrinking. It has an area-slot with a
    value area.2.

  • I need to calculate area.2/area.1 = percent after shrink.

I know that to get the square-meters etc. of the polygons, they must first be transformed to the appropriate projection. But I hope to avoid to do the transformation, to save calculation-time.

To calculate the percentage-difference, can I use the projection above? I’m not sure what unit the areas are in, square-degrees? Each polygon is less than 200x200 meters.

share|improve this question

3 Answers 3

up vote 1 down vote accepted

As @radouxju wrote, you really should use an equal area projection. Interesting is the question of the ratio of the same polygon with shrinked area. I've made test and shift a rudimentray 200x200m patch and the shrinked version 100x100m along the latitude axis. As you see in the results the ration will be something near 25% BUT THERE ARE SOME SERIOUS PROBLEMS

  1. The area in LonLat for objects < 200mx200m is very small and you have to deal with rounding errors?

  2. The area not a real area because you underestimate it by missing the great circle calculation.

#!/usr/bin/perl -w
use strict;
use Geo::Proj4;

# Initialize a transverse mercator
my $proj = Geo::Proj4->new(
    proj => "tmerc", 
    ellps => "WGS84"
);

# Shifting the patch in latitude and check the 
# differences in size and ratio metric versus lonlat
for my $offs (-4..4) {

    # calculate the offset im Meters
    my $loffs = $offs*2000000;

    # getting the center latitude
    my ($la, $lo) = $proj->inverse(0, $loffs);  

    # create a rudimentary patch 200 by 200 in lon lat
    my ($x1,$y1) = (-100, $loffs-100);  
    my ($x2,$y2) = ( 100, $loffs+100);  
    my ($la1, $lo1) = $proj->inverse($x1,$y1);
    my ($la2, $lo2) = $proj->inverse($x2,$y2);

    # create a rudimentary patch 100 by 100 lon lat
    my ($x3,$y3) = (-50, $loffs-50);    
    my ($x4,$y4) = ( 50, $loffs+50);    
    my ($la3, $lo3) = $proj->inverse($x3, $y3); 
    my ($la4, $lo4) = $proj->inverse($x4, $y4); 

    # calculate the area for both
    # This is not exact because we have great circles
    my  $la1 = ($lo2-$lo1)*($la2-$la1);
    my  $la2 = ($lo4-$lo3)*($la4-$la3);
    my  $lr  = $la2/$la1*100.0; 

    # calculate the areas in tmerc
    my $ma1 = ($x2-$x1)*($y2-$y1);
    my $ma2 = ($x4-$x3)*($y4-$y3);

    # calculate the ratio
    my $mr  = $ma2/$ma1*100.0; 

            # result
    print "OFFS: $loffs  METRIC A1: $ma1 A2: $ma2 RATIO: $mr \n";
    print " LAT: $la  METRIC A1: $la1 A2: $la2 RATIO: $lr  \n\n";
}

# EOF -------------------------------------------------------

Results:

OFFS: -8000000  METRIC A1: 40000 A2: 10000 RATIO: 25 
 LAT: -72.0705339311268  METRIC A1: 1.04284937849616e-05 A2: 2.60712344564335e-06 RATIO: 24.9999999942749  

OFFS: -6000000  METRIC A1: 40000 A2: 10000 RATIO: 25 
 LAT: -54.126533597883  METRIC A1: 5.49683618861433e-06 A2: 1.37420904709556e-06 RATIO: 24.9999999989444  

OFFS: -4000000  METRIC A1: 40000 A2: 10000 RATIO: 25 
 LAT: -36.130292919723  METRIC A1: 4.00472249368946e-06 A2: 1.00118062344132e-06 RATIO: 25.0000000004732  

OFFS: -2000000  METRIC A1: 40000 A2: 10000 RATIO: 25 
 LAT: -18.0814780948804  METRIC A1: 3.41404470186791e-06 A2: 8.53511175486969e-07 RATIO: 25.0000000005856  

OFFS: 0  METRIC A1: 40000 A2: 10000 RATIO: 25 
 LAT: 0  METRIC A1: 3.2496356947654e-06 A2: 8.1240892371699e-07 RATIO: 25.000000000789  

OFFS: 2000000  METRIC A1: 40000 A2: 10000 RATIO: 25 
 LAT: 18.0814780948804  METRIC A1: 3.41404470186791e-06 A2: 8.53511175486969e-07 RATIO: 25.0000000005856  

OFFS: 4000000  METRIC A1: 40000 A2: 10000 RATIO: 25 
 LAT: 36.130292919723  METRIC A1: 4.00472249368946e-06 A2: 1.00118062344132e-06 RATIO: 25.0000000004732  

OFFS: 6000000  METRIC A1: 40000 A2: 10000 RATIO: 25 
 LAT: 54.126533597883  METRIC A1: 5.49683618861433e-06 A2: 1.37420904709556e-06 RATIO: 24.9999999989444  

OFFS: 8000000  METRIC A1: 40000 A2: 10000 RATIO: 25 
 LAT: 72.0705339311268  METRIC A1: 1.04284937849616e-05 A2: 2.60712344564335e-06 RATIO: 24.9999999942749  
share|improve this answer

It is not possible to convert from an area in degrees to an area in meters. You can indeed have two polygons with the same area in degree but different area once projected (it depends on the size and on the orientation).

If you want consistent results, you should really start projecting your polygons in an "equal area" projection (e.g. sinusoidal, cylindrical equal area...).

EDIT : Within a 200 by 200 m area, what you can assume is that the ratio between X and Y distances for one degree is constant. The area ratios will thus be approximately the same, as W. Huber said.

share|improve this answer
    
Thanks for quick response. Can one assume that within the same 200x200m on the ground, that the proportion between decimal-polygon-to-decimal-polygon is same as if first converted to m? –  Chris Feb 20 at 21:00
1  
Chris, that is correct provided you are more than a kilometer or so away from one of the poles :-). –  whuber Feb 20 at 21:20

The PolygonArea class and the Planimeter utility in GeographicLib (written by me) allow you to compute the area of a polygon whose vertices are given in terms of latitudes and longitudes and whose edges are ellipsoidal geodesics. The accuracy is around 0.001 m2. These routines are available in C++, C, Fortran, Java, Python, MATLAB/Octave, Javascript, IDL, C#, but not, unfortunately, R. There's an online version of Planimeter here, if you'd like to try it out.

There are two alternative methods each of which entails an approximation:

(1) Transform the coordinates to an equal-area projection. This entails an error in that a geodesic is typically mapped to a curved line in the projection. (This error is small if the edges of your polygon are short.)

(2) Transform your polygon to the authalic sphere (geographic latitude goes to authalic latitude) and compute the area of the resulting spherical polygon. Just be sure not to use L'Huilier's formula for the area of spherical triangles, because it's badly conditioned. Better is to use the formula for a spherical quadrilateral, the last formula of this Wikipedia article. This method results in a small error because a geodesic is typically not mapped to a great circle.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.