[Proj] Lambert Azimuthal Equal Area scale

[Lourens Veen]

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