Sunday, July 17, 2016

more HiRISE conversion tips (until labels are fixed)

The HiRISE Team is still planning to fix a few label issues so these work arounds may go away.

background explanation for HiRISE projection issues - from wayback. See this post for more on how ISIS calculates a local scaling radius for equirectangular (aka equidistant cylindrical).

First you need to get some excellent DTMs (aka DEMs) and their orthorectified GeoJpeg2000 images: http://hirise.lpl.arizona.edu/dtm/

HIRISE label issues
(1) MOLA-orthorectified RDRs
For images without partnered DTMs (typical MOLA-rectified images), they have an older map projection definition (where std par is defined as the center lat). There is a simple script called "fix_jp2_v2" to switch the values in-place in the Jpeg2000 header. This updated version "v2" (Jan. 2023) can safely be run on older and newer (already corrected) HiRISE labels. For newer labels, the program will safely skip over them since they are corrected.
original source code (now outdated. This version should not be used on newer HiRISE image as they are already corrected.): https://trac.osgeo.org/gdal/ticket/2706

> gcc -o fix_jp2_v2 fix_jp2_v2.cpp
> fix_jp2_v2 PSP_001918_1735_RED_C_01_ORTHO.JP2
Success, file updated.

(2) partnered DTM-orthorectified images
2.1. first run fix_jp2_v2 as in step 1.

2.2. So a more subtle issue (only in the partnered-DTM ortho-images *.JP2s), is that the radius is not set correctly in the geoJpeg2000 header. This doesn't impact the other map projection parameters (like cellsize or offsets) -- it is really just an incorrect value. To fix this for other GDAL environments you need to override the radius in the GeoJpeg2000 file (this generated vrt header also works for ArcMap).

using gdal tools you can create a corrected virtual header for the GeoJpeg2000 file using the projection defined in the DTM.
> gdalsrsinfo DTEEC_001918_1735_001984_1735_U01.IMG -o wkt > fixme.prj

>cat fixme.prj
PROJCS["EQUIRECTANGULAR MARS",GEOGCS["GCS_MARS",DATUM["D_MARS",SPHEROID["MARS_localRadius",3396040,0]],PRIMEM["Reference_Meridian",0],UNIT["degree",0.0174532925199433]],PROJECTION["Equirectangular"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",283.07],PARAMETER["standard_parallel_1",-5],PARAMETER["false_easting",0],PARAMETER["false_northing",0]]

>gdal_translate -of vrt -a_srs fixme.prj -a_nodata 0 PSP_001918_1735_RED_C_01_ORTHO.JP2 PSP_001918_1735_RED_C_01_ORTHO.vrt
Input file size is 6368, 10767

This process will work for the full-res "A" ortho files too, e.g. PSP_001918_1735_RED_A_01_ORTHO. I am also setting the image's nodata value in the command above too. Just helps with display.

------------------
BTW, if you don't like Jpeg2000 images you can use the HiRISE team JP2_PDS to convert to raw PDS format or the free Kakadu tools to convert to GeoTiff/GeoBigTiff. e.g.

> kdu_expand -i PSP_001918_1735_RED_C_01_ORTHO.JP2 -o PSP_001918_1735_RED_C_01_ORTHO.tif

-----------------

To map project these all to the same coordinate system (use the corrected vrt) and get the latest GDAL binaries and for each file run the lines below (for GDAL binaries, I like Anaconda Python environment [once installed type: 'conda install gdal'] or Kyng Chaos for Mac, or Osgeo4W - Windows). These steps should work for DTM *.IMG files also.

-only if you don't have GDAL with jpeg2000 support:
> kdu_expand -i PSP_001918_1735_RED_C_01_ORTHO.JP2 -o PSP_001918_1735_RED_C_01_ORTHO.tif
-- again fix radius for converted Tiff
> gdalsrsinfo DTEEC_001918_1735_001984_1735_U01.IMG -o wkt > fixme.prj
> gdal_translate -of vrt -a_srs fixme.prj PSP_001918_1735_RED_C_01_ORTHO.tif PSP_001918_1735_RED_C_01_ORTHO.vrt

in Degrees (equivalent ~1m/p):--download "Mars 2000 Sphere.prj"
> wget http://webgis.wr.usgs.gov/pigwad/ArcMap_prjs/Solar_System/Mars%202000%20Sphere.prj
> gdalwarp -t_srs "Mars 2000 Sphere.prj" -tr 1.70670654e-5 1.70670654e-5 -r bilinear PSP_001918_1735_RED_C_01_ORTHO.vrt PSP_001918_1735_RED_C_01_ORTHO_degrees.tif

Or in Meters (e.g. Simple Cylindrical or Equirectangular w/ center at 0, out cellsize set to 1m):
> echo 'PROJCS["Mars_Equidistant_Cylindrical",GEOGCS["Mars 2000",DATUM["D_Mars_2000",SPHEROID["Mars_2000_IAU_IAG",3396190.0,0]],PRIMEM["Reference_Meridian",0],UNIT["Decimal_Degree",0.0174532925199433]],PROJECTION["Equidistant_Cylindrical"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",0],PARAMETER["Standard_Parallel_1",0],UNIT["Meter",1]]' > Mars2000_Equi0.prj

> gdalwarp -t_srs Mars2000_Equi0.prj -tr 1 1 -r bilinear PSP_001918_1735_RED_C_01_ORTHO.vrt PSP_001918_1735_RED_C_01_ORTHO_1meter.tif

-------------------

Tip: To check out your handy work for the  DTM and fixed Ortho use TuiView or 3D Analyst. For 3D Analyst, simply load the images in and under the DTM and Ortho, set the Base Height Tab to use the DTM - almost instant 3D: https://pdsimage2.wr.usgs.gov/pub/pigpen/mars/tutorials/arcScene_Hirise.jpg