[Proj] Lambert Azimuthal Equal Area scale

Wed Jun 10 08:41:54 PDT 2015

Thanks for you help, and incidentally also for GeographicLib! I'm using 
it to calculate the pixel side lengths from a C++ program, and found it 
nicely designed and easy to use. But then I went on to check the output 
against that of invgeod to make sure I hadn't made any mistakes...

I'll check everything again tomorrow, this time with lon and lat in the 
right order for invgeod to understand it correctly. I agree that geod 
should be fixed, both geod and proj are libproj utilities, they should 
correctly interpret each other's default output format.

Thanks again,

Lourens

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

-- 
ir. Lourens Veen - PhD Student - UvA-IBED-CGE
    m: P.O. Box 94248, NL-1090 GE  Amsterdam
              t: +31 (0)20-5258293