Q:

Lat/lon in km in GeoTIFF medata file not matching calculated distance

I am extracting information from a GeoTIFF file using gdalinfo based on THIS and I get the following:
Driver: GTiff/GeoTIFF
Files: myfile.tif
       myfile.tif.aux.xml
Size is 953, 2824
Coordinate System is:
PROJCS["Lambert_Conformal_Conic_2SP_World_Geodetic_System_1984",
    GEOGCS["GCS_World_Geodetic_System_1984",
        DATUM["D_World_Geodetic_System_1984",
            SPHEROID["WGS_1984",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",66],
    PARAMETER["standard_parallel_2",78],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-42],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-432678.541899999952875,9726108.134099999442697)
Pixel Size = (300.000000000000000,-300.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -432678.542, 9726108.134) ( 53d58'12.29"W, 70d52'36.63"N)
Lower Left  ( -432678.542, 8878908.134) ( 50d38' 9.28"W, 63d22'47.25"N)
Upper Right ( -146778.542, 9726108.134) ( 46d 6'31.44"W, 71d13'15.48"N)
Lower Right ( -146778.542, 8878908.134) ( 44d56'51.10"W, 63d37'31.84"N)
Center      ( -289728.542, 9302508.134) ( 48d45'16.75"W, 67d18'24.75"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
  Min=0.900 Max=1.100 
  Minimum=0.900, Maximum=1.100, Mean=0.993, StdDev=0.025
  NoData Value=-3.4028234663852886e+38
  Metadata:
    RepresentationType=ATHEMATIC
    STATISTICS_MAXIMUM=1.1000000238419
    STATISTICS_MEAN=0.99293445610505
    STATISTICS_MINIMUM=0.89999997615814
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=0.025099979310242

I am trying to reproduce the km distances of the 4 corners based on the lat/lon degree references. I can get the longitude, but somehow my latitudes are way off... Here is how I do it:
import geopy
from geopy import distance
from geopy.location import Point
p1 = Point('''70 52'36.63" N 53 58'12.29" W''')
p2 = Point('''0 N 53 58'12.29" W''') #Same longitude, but latitude set to 0
distance.distance(p1,p2).m # WGS84 distance in m  

Base don the gdalinfo output, I would expect to get back 9726108.134, but instead, I get out 7866807.760742959. I think it's off by too much to be a projection mistake. Any other suggestions of what I'm doing wrong?
Note that when I'm reproducing the km in THIS EXAMPLE, I can get the latitude in km right with my method...

A:

Your distance check makes two assumptions:

the point 0° N, 0° E is also at 0, 0 in the re-projection.
the x and y axis have the same alignment in both projections

However, both of these do not hold true. Here is a plot to illustrate the difference between the two projections:

The left plot shows the area which is depicted in your GeoTiff as a green polygon in WGS84 and the right plot shows the same thing in the re-projection. As you can see point 0° N, 0° E ist now at 7633606.089886177, 2779916.995372205 (while Point 0, 0 of the-re-projection would be at 0° N, -42° E in WGS84). Furthermore, the image area, which forms a perfect rectangle on the right is distorted in WGS84.
Long story short: the check you performed will not come up with your expected result, due to the differences inherent to the projections. As an approximation, you could instead use the coordinates of the four cornes to do your check and the the results will turn out to be (roughly) correct. Here's an example using the upper left and lower left corner:
p1 = geopy_Point('''70 52'36.63" N 53 58'12.29" W''')
p2 = geopy_Point('''63 22'47.25" N 50 38' 9.28" W''')

distance.distance(p1,p2).m

Which returns:
848195.724897564

This is very close to the expected distance (9726108.134 - 8878908.134 = 847200) and just slightly off due to rounding of the given WGS84 coordinates.
Hope this helps!
EDIT:
Here's a suggestion how you can transform from WGS84 to the re-projection:
import osr
from pyproj import Proj, transform

def reproject(x_in, y_in):
    proj_in = Proj(init='epsg:4326')
    osr_spat_ref = osr.SpatialReference()
    osr_spat_ref.ImportFromWkt('''PROJCS["Lambert_Conformal_Conic_2SP_World_Geodetic_System_1984",
    GEOGCS["GCS_World_Geodetic_System_1984",
        DATUM["D_World_Geodetic_System_1984",
            SPHEROID["WGS_1984",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",66],
    PARAMETER["standard_parallel_2",78],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-42],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]''')
    spat_ref_pyproj = osr_spat_ref.ExportToProj4()
    proj_out = Proj(spat_ref_pyproj)
    out = transform(proj_in, proj_out, x_in, y_in)
    return out

Just enter the the longitude and latitude coordinates as x_in and y_in respectively to obtained the re-projected coordinates. From these you can then figure out which tile of the image you are interested in by deducting the coordinates of the corners.

