Using GDAL for conversions between Lon Lat <> Sample Line <> X Y meters

First, the best method to get GDAL v3.4 binaries is to install Anaconda (any OS) and then:

conda create --name gdal
conda activate gdal
conda config --env --add channels conda-forge
conda install -c conda-forge gdal

for more conda information.

Note GDAL v3.4 or higher is required to support the new IAU_2015 codes. You don’t have to have the codes but makes it easier – see Data Workshop abstract (note abstract references IAU:2018 in error).

Anyway, once GDAL is installed there is a routine called gdallocationinfo which has methods to request line/sample or lon/lat or X/Y meters to sample/line from an image file. This program is really for reporting the pixel value, but does report what pixel (sample/line) the value was attained from.

Now to go from lon/lat to X/Y meters or back there is a routine called gdaltransform to do that.


gdallocationinfo Examples:

for Venus (code 29900 is Venus in degrees) – this example goes from Lon Lat to Sample Line:

> gdallocationinfo -l_srs IAU_2015:29900 Venus_Magellan_Topography_Global_4641m_v02.cub 120 45
Report:
  Location: (6826P,1024L)
  Band 1:
    Value: 646

– note this code is only available in the latest version of GDAL (v3.4). You can test using:

> gdallocationinfo --version

for the files internal map projection X,Y meters to Sample Line.
This next example is using X= -10,000,000.0, Y=0.0 (meters)

> gdallocationinfo -geoloc Venus_Magellan_Topography_Global_4641m_v02.cub -10000000 0
Report:
  Location: (1941P,2048L)
  Band 1:
    Value: 318

for Sample Line to Sample Line (why? it reports the pixel value at that location).

> gdallocationinfo Venus_Magellan_Topography_Global_4641m_v02.tif 2000 1000
Report:
  Location: (2000P,1000L)
  Band 1:
    Value: -424

gdaltransform Examples:

gdaltransform actually fires up an interactive session if you don’t supply it with a list of coordinates. You can also “echo X Y” in geotransform on the command-line (see below). -s_srs is the “source” spatial reference system and -t_srs is the target. Note if you don’t have GDAL v3.4, the code IAU_2015:29900 = “+proj=longlat +R=6051800 +no_defs” (means Venus in degrees). It will also try to report X,Y,Z (and even time) but here Z will always be reported as 0.

– Shown below, it will project from degrees to Sinusoidal meters (center meridian=35)

> gdaltransform -s_srs IAU_2015:29900 -t_srs "+proj=sinu +lon_0=35 +R=6051800"
Enter X Y [Z [T]] values separated by space, and press Return.
180 0  //IN Lon Lat
15315456.172468 0 0  //OUT in meters
180 45   //IN Lon Lat
10829662.9165175 4753072.60524868 0  //OUT in meters

OR without GDAL 3.4 IAU_2015 codes (use short proj string):

> gdaltransform -s_srs "+proj=longlat +R=6051800" -t_srs "+proj=sinu +lon_0=35 +R=6051800"

This is Sinusoidal meters to degrees:

> gdaltransform -s_srs "+proj=sinu +lon_0=35 +R=6051800" -t_srs IAU_2015:29900
Enter X Y [Z [T]] values separated by space, and press Return.
10829662.9165175 4753072.60524868  //IN X Y meters 
180.000000000001 45 0  //OUT in degrees

using echo at the command-line:

echo 10829662.9165175 4753072.60524868 0 | gdaltransform -s_srs "+proj=sinu +lon_0=35 +R=6051800" -t_srs IAU_2015:29900
180.000000000001 45 0

gdalsrsinfo Examples:

BTW, to check out those IAU_2015 codes in degrees you can run

> gdalsrsinfo IAU_2015:29900

PROJ.4 : +proj=longlat +R=6051800 +no_defs

OGC WKT2:2018 :
GEOGCRS["Venus (2015) - Sphere / Ocentric",
    DATUM["Venus (2015) - Sphere",
        ELLIPSOID["Venus (2015) - Sphere",6051800,0,
            LENGTHUNIT["metre",1]],
        ANCHOR["Ariadne: 0.0"]],
    PRIMEM["Reference Meridian",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["IAU",29900,2015],
    REMARK["Source of IAU Coordinate systems: doi://10.1007/s10569-017-9805-5"]]

To see a Sinusoidal projection try:

> gdalsrsinfo IAU_2015:29920
PROJ.4 : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6051800 +units=m +no_defs

OGC WKT2:2018 :
PROJCRS["Venus (2015) - Sphere / Ocentric / Sinusoidal, clon = 0",
    BASEGEOGCRS["Venus (2015) - Sphere / Ocentric",
        DATUM["Venus (2015) - Sphere",
            ELLIPSOID["Venus (2015) - Sphere",6051800,0,
                LENGTHUNIT["metre",1]],
            ANCHOR["Ariadne: 0.0"]],
        PRIMEM["Reference Meridian",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["IAU",29900,2015]],
    CONVERSION["Sinusoidal, clon = 0",
        METHOD["Sinusoidal",
            ID["PROJ","SINUSOIDAL"]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    ID["IAU",29920,2015]]

PDS4 tips

For map projections, PDS4 requires the CART (Cartography) Discipline Dictionary. GDAL can help with the generation of a “map projected” PDS4 label (which is held in the CART section). You don’t need to manually calculate any X/Y map projection offset in meters If the input file has a projection set (and is understood by GDAL). But if not, you can help the conversion along.

So to use that Sinusoidal code in GDAL (again requires the input file has a projection or geospatial transform – cellsize and map projection offsets).

gdal_translate -of PDS4 -a_srs IAU_2015:29920 input_pds3.img output_pds4.xml
– basic example which should require a input template (hence the warnings listed).

recall to go from PDS3 to PDS4:

cart:upperleft_corner_x = (SAMPLE_PROJECTION_OFFSET + 0.5) * (MAP_SCALE * 1000) * -1
cart:upperleft_corner_y = ( LINE_PROJECTION_OFFSET  + 0.5) * (MAP_SCALE * 1000)

where SAMPLE_PROJECTION_OFFSET, LINE_PROJECTION_OFFSET, and MAP_SCALE (in km) are
from the original IMAGE_MAP_PROJECTION group of the PDS3 label. If MAP_SCALE is
not in km the convert to meters accordingly.
from (user doc): https://github.com/pds-data-dictionaries/ldd-cart

To check the PDS4 CART is good, GDAL can also try to parse the PDS4 label. It should spit out a large text block as a PROJ4 string and Well-known Text (WKT) projection:

> gdalsrsinfo input_PDS4.xml

1 Like