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