Skip to content

Using GDAL for Conversions Between Lon Lat _ Sample Line _ X Y meters

Jump to bottom
Trent Hare edited this page Mar 2, 2022 · 1 revision

First, the best method to get GDAL v3.4 binaries is to install Anaconda (any OS) and then "conda install gdal". see: https://opensourceoptions.com/blog/how-to-install-gdal-with-anaconda/ 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.

Anyway, once GDAL is install there is a routine called gdallocationinfo which has methods to request line/sample or lon/lat or X/Y meters to sample/line for 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.tif 120 45

-- 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 example is using X= -10,000,000.0 Y=0.0 (meters).

gdallocationinfo -geoloc Venus_Magellan_Topography_Global_4641m_v02.tif -10000000 0

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.

--This will be 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:

To Check CART is good this might help. It should spit out a large text block as a PROJ4 string and Well-known Text (WKT) projection:

gdalsrsinfo input_PDS4.xml

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

gdalsrsinfo IAU_2015:29900

To see sinusoidal in action

gdalsrsinfo IAU_2015:29920

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 (pixel offsets) to PDS4 (X,Y meter offsets):

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

Clone this wiki locally