[Proj] Lambert Azimuthal Equal Area scale

[Charles Karney]

Bletch...  It seems that geod expects latitude before longitude ignoring
the strong hint you gave it with the W and N.  At least with proj 4.9.0
(which incorporates a new backend), you get nan's in this case to
give you a hint that's something is amiss.  The GeographicLib utility
for geodesic calculations, GeodSolve, does pay attention to the
hemisphere designators.  This would have told you that your pixel is
approximately a rectangle 11km x 9km.  geod/invgeod should be fixed!

   --Charles

On 06/10/15 09:42, Charles Karney wrote:
> 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
>>