[Proj] Lambert Azimuthal Equal Area scale

[Charles Karney]

I can confirm that the area of your pixel is close to 100 km^2.  You can
too.  Visit

http://geographiclib.sourceforge.net/cgi-bin/Planimeter

Paste the coordinates of the corners of your pixel

164d11'36.762"W	54d52'47.599"N
164d1'53.146"W	54d50'37.423"N
163d55'19.442"W	54d54'20.055"N
164d5'3.703"W	54d56'30.694"N

into the text box.  Hit "submit" and the area comes back as

   99999985.2 m^2

This is the true area of the geodesic quadrilateral on the ellipsoid.

On 06/10/15 09:23, Lourens Veen wrote:
> Dear all,
>
> I have a GeoTIFF file that uses a custom Lambert azimuthal equal area
> projection with proj.4 string
>
> +proj=laea +lat_0=15 +lon_0=-80 +x_0=0 +y_0=0 +datum=WGS84 +units=m
> +no_defs +ellps=WGS84 +towgs84=0,0,0
>
> The pixels are sized 10000x10000m, and so in the projected system have a
> surface area of 100km^2.
>
> I am now using the invproj tool (from proj 4.8.0) to un-project four
> corners of a pixel to WGS84 coordinates:
>
> $ invproj +proj=laea +lat_0=15 +lon_0=-80 +x_0=0 +y_0=0 +datum=WGS84
> +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 -S
> -4600000 6190000
> 164d11'36.762"W	54d52'47.599"N	<1.25548 0.798216 1 26.0103 1.25729 0.795362>
> -4600000 6180000
> 164d1'53.146"W	54d50'37.423"N	<1.25483 0.798516 1 25.943 1.25653 0.795842>
> -4590000 6180000
> 163d55'19.442"W	54d54'20.055"N	<1.25437 0.798719 1 25.8934 1.25597 0.796196>
> -4590000 6190000
> 164d5'3.703"W	54d56'30.694"N	<1.25502 0.798414 1 25.9607 1.25673 0.795716>
>
> Next, I'm calculating the distance-along-the-ellipsoid between these
> points using invgeod:
>
> $ invgeod +ellps=WGS84 +units=m
>
> 164d11'36.762"W	54d52'47.599"N 164d1'53.146"W	54d50'37.423"N
> 12d6'11.586"	-167d53'12.778"	18457.232
> 164d1'53.146"W	54d50'37.423"N 163d55'19.442"W	54d54'20.055"N
> -28d30'57.832"	151d28'0.716"	13856.109
> 163d55'19.442"W	54d54'20.055"N 164d5'3.703"W	54d56'30.694"N
> -167d51'52.105"	12d7'31.893"	18479.179
> 164d1'53.146"W	54d50'37.423"N 164d5'3.703"W	54d56'30.694"N
> -119d16'41.52"	60d41'41.449"	12043.634
> 164d5'3.703"W	54d56'30.694"N 164d11'36.762"W	54d52'47.599"N
> 151d21'28.799"	-28d37'30.228"	13847.743
> 164d1'53.146"W	54d50'37.423"N 164d5'3.703"W	54d56'30.694"N
> -119d16'41.52"	60d41'41.449"	12043.634
> 163d55'19.442"W	54d54'20.055"N 164d11'36.762"W	54d52'47.599"N
> 174d47'56.288"	-5d11'38.318"	30345.538
>
> The LAEA projection doesn't preserve angles or distances, so it's no
> surprise that my square pixel is warped. From the above lengths of the
> sides and diagonals, it's approximately a parallelogram, with top and
> bottom sides of 13.8 km, left and right sides 18.4 km, and the top-right
> corner almost directly above the bottom-left corner (i.e. it's leaning
> left about 45 degrees).
>
> What is surprising, at least to me, is that this parallelogram has a
> surface area of about 165 km^2! Of course in reality the parallelogram
> is embedded in the ellipsoid, and therefore not flat, but surely over a
> distance of 30 km the curvature can't be enough to give such a
> deviation? Besides, that should result in a surface area less than the
> actual 100 km^2 I think.
>
> Is there something fundamentally wrong with this calculation? Am I using
> invproj and invgeod incorrectly? What am I missing?
>
> Thanks in advance,
>
> Lourens
>
> _______________________________________________
> Proj mailing list
> Proj at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj
>

-- 
Charles Karney <charles.karney at sri.com>
SRI International, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662