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

Question

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...

Hoping @fmark has a suggestion on this issue? – Cynthia GS Nov 05 '19 at 15:47 — Cynthia GS, Nov 05 '19 at 15:47

bexi · Accepted Answer · 2019-11-05T16:46:13.167

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.

Thanks, this does help a little, but ultimately, I'm hoping to give my code degree lat and lon, and find the corresponding pixel in my GeoTIFF file. The way I'm doing it is by calculating the km lat and lon from "some" origin (so far been using (0,0)) and then finding the closest match in my lat and lon arrays for the GeoTIFF file. So I'm still confused about how to achieve that... — Cynthia GS, Nov 05 '19 at 16:05
@CynthiaGS: I added a little script to transform the coordinates, which should help you find the closest tile in the GeoTiff — bexi, Nov 05 '19 at 16:17
Thanks so much, this is way more than I was hoping for! Seems to work as is without having to deduct anything. — Cynthia GS, Nov 05 '19 at 19:32
Happy to help. I tend to convert the GeoTiffs to numpy arrays which then requires converting the coordinates to x/y index of the array when I want to obtain data for specific locations. Just out of curiosity: how do you access data at specific coordinates (i.e. what library/method)? — bexi, Nov 06 '19 at 09:29
Using gdal to convert to a numpy array as well: `import gdal src_ds_2 = gdal.Open("myfile.tif" ) Data = np.array(src_ds_2.ReadAsArray())` — Cynthia GS, Nov 06 '19 at 16:20

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

1 Answers1