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.
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
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
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.
> 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.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.
------------------
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
-----------------
> 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
Thank you Trent for the helpful instruction and script for fixing the coordinate issue of HiRISE jp2s. May I ask if there is a Polar_Stereographic version of the fix_jp2 available? For example, the current fix_jp2 wouldn't work for this scene: https://hirise.lpl.arizona.edu/PSP_002081_1055
ReplyDeletewhich is not in Equirectangular projection. Many thanks!!
HiRISE polar files don't need any updates. They are good as is.
ReplyDeleteHi again Trent,
ReplyDeleteIs there a simple way to identify which JP2 images need to use fix_jp2 and which ones not?
For example, if I batch process (using fix_jp2) a group of HiRISE orthoimages *_RED_A_01_ORTHO.JP2, most of the JP2 images get corrected correctly w.r.t the DTMs (e.g., https://www.uahirise.org/dtm/dtm.php?ID=ESP_019604_2360), however, occasionally, I get a few of them being incorrectly corrected by fix_jp2 (e.g., https://www.uahirise.org/dtm/dtm.php?ID=ESP_069665_2055).
This means the original JP2 images (e.g., https://www.uahirise.org/dtm/dtm.php?ID=ESP_069665_2055) have no problem and can be displayed (e.g., in QGIS) correctly w.r.t the DTMs WITHOUT using fix_jp2, but they cannot be correctly displayed any longer after using fix_jp2 (an error is generated by fix_jp2 to the JP2 image).
My guess is that the HiRISE team has now somehow fixed the JP2 issue for the "new" (?since when?) HiRISE releases and continuing to apply fix_jp2 to the "newly" released JP2 images leads to a newly generated geometric error.
Do you have any idea on this? If my guess is right, do you happen to know which JP2 images I should apply fix_jp2 (released until when?) and which JP2 images I should not apply fix_jp2 (released since when?) ?
Many thanks and best regards,
Yu
I have updated a new version, now called fix_jp2_v2 which doesn't touch the newer already HiRISE-team corrected files. This routine will eventually no longer be needed when ALL the HiRISE images are re-released, but that is a ways off.
ReplyDelete