Archiving SPC

Run backplanesGMT (which is in ORExSPCsupport) (or backplanesISIS)

Topography

Run backplanesGMT

  # Make the backplanes
name=FIRED1
~/bin/backplanesGMT
FIRED1 // give bigmap name
1737.4 // give datum, 1737.4 is for Luna

Output:

 Version:   1.50000000    
 Input map name (only 6 char no MAPFILES and .MAP)
 BIGMAP position
     V   120.38036346435547       -510.48150634765625       -54.750602722167969     
     Ux  -2.9186813160777092E-002  0.12163921445608139      -0.99214518070220947     
     Uy  0.97239929437637329       0.23332308232784271        0.0000000000000000     
     Uz  0.23149037361145020      -0.96476125717163086      -0.12509183585643768     
     QSZ          99
 Lon   279.612366       286.917664    
 Lat  -9.55961323      -2.35706925

<Deprecated> We are currently not using blockmean and sphinterpolate, but this section is maintained as another set of procedures.

We'll need the Lon and Lat range for "argR", and QSZ will give an idea of what you need for "argI".

  # Convert them into lon/lat/val
  paste $name-lonlat.txt $name-r.txt > tmp.ll
  var=tmp.ll
  argR="-R279.612366/286.917664/-9.55961323/-2.35706925"
  argI="-I199+n"

  # Run the rebinning/interpolating routines
  echo "Blockmean"
  gmt blockmean $var $argR $argI  > out
  echo "sphereinterpolate"
  cat out |  gmt sphinterpolate $argR $argI -Gtest.nc

Now you need some math to convert the Lon range number of nodes (i.e. "argI") to deg/px. (286.917664 - 279.612366) / 199 = 0.03671

  # Convert to geoTiff
  echo "gdal translate"
  degPx=0.0367
  gdal_translate -of ISIS3 -tr $degPx $degPx -r bilinear -b 1 -co TARGET_NAME=Tethys -co DATA_LOCATION=GEOTIFF -a_srs support/60300.prj NETCDF:test.nc $name-r-ISIS.lbl

If you don't want to make it readable by ISIS, and don't care about square pixels, you can run the following

     gdal_translate -of GTiff -b 1 -a_srs support/60300.prj  NETCDF:test.nc $name-r.tif

See above in "Grid to GeoTIF" for details, but a better program to use is gdalwarp so you can properly define the projection. Plus this allows the pixel size to be defined in meters, which is cleaner. The example below if for the Moon (Luna)

gdalwarp -of PDS4 -co IMAGE_FORMAT=GEOTIFF -r bilinear -tr $metPx $metPx -s_srs support/30100.prj  -t_srs support/30110.prj  NETCDF:tmp.nc ${bigmap}-${type}-v${ver}.xml

How to trim size

Use ArcMap to find the new lat and lon ranges
Remake geotiff with new lat and lon ranges
Measure the size of the lat range. Divide this number by the GSD from SPC to get number of pixels.
Divide number of pixels by lon range to get degPx
Update config.txt with these new values (i.e. lat and lon ranges, number of pixels across, degPx)
Make geotiff again

</Deprecated> End section on blockmean and sphinterpolate

Current procedures for generating Geotiff

JRW has a shell script that will make these products, and their PDS label. It is found in /usr/local/spc/bin/makeGeoTiffBigNneighbor.sh. The shell script isn't needed, but can make the process easier.

Make a config.txt file for the shell script that looks something like this...

<rand> INGV02$ cat config.txt 
BIGMAP=INGV02
ARGR=-R161.625107/161.997177/-35.9918938/-35.6902657
ARGI=-I1821+n
METPX=5
PROJECT=30110
BODY=Moon
RUNALL=No
VERSION=1

METPX is the GSD in meters for the bigmap. Project 30110 is for Equidistant_Cylindrical on the Moon. Project 30118 is North_Pole_Stereographic if you ever make a North Polar region.

# Make input for GMT
bigMapName=BSTDTM
paste $bigMapName-lonlat.txt $bigMapName-r.txt > tmp.ll
var=tmp.ll

gmt nearneighbor $var -Gtmp.nc $argI $argR -S0.001d -N16+m8


# Make normal GeoTiff
echo "gdal warp"
gdalwarp -of PDS4 -dstnodata -3.40282306073709653e+38 -co IMAGE_FORMAT=GEOTIFF -r bilinear -tr $metPx $metPx -s_srs support/30100.prj  -t_srs support/$proj.prj  NETCDF:tmp.nc ${bigmap}-${type}-v${ver}.xml

#Make GeoTiff readable by ISIS
gdalwarp -of ISIS3 -dstnodata -3.40282306073709653e+38 -tr $metPx $metPx -r bilinear -co TARGET_NAME=$body -co DATA_LOCATION=GEOTIFF -s_srs support/30100.prj  -t_srs support/$proj.prj NETCDF:tmp.nc ${bigmap}-${type}-ISIS-v${ver}.lbl

gmt nearneighbor only writes out a value if data points exist in the specific sectors, then makes an average of the nearest values in those sectors. -S.001d is the search radius of each sector (in the example it is .001 degrees) and -N16+m8 means there are 16 sectors and there must be values in 8 of those sectors to write out a value.

Project 30100 is the data in degrees. Note that if you've previous run ISIS by using "conda activate isis" in the window you are using, then gdalwarp will throw a project error. This is true even if you've run "conda deactivate", so best to just start with a fresh window.

The "-dstnodata" option in gdalwarp is to use a small number for "no data" instead of NaN. NaN is not compatible with PDS4.

prj files can be found here

IMPORTANT: After you've make the geoTiff and it's label, the label is still not ready to go because there will be placeholder values. To fix those, use a script Trent Hare sent JRW on 27 June 2022, which works for LRO data. PDS4_replaceStrings.py into the folder (note it will change ALL xml files in the directory and ALL subdirectories) and run it with "python PDS4_replaceStrings.py". The python file for LRO products can be found here

How to create bigmap products and labels for ISIS cubes

Make the cubes using backplanesISIS and ascii2isis. Now make the cubes Band sequential, which is something ascii2isis v4.1.1 cannot do, despite the ISIS manual saying it can. Instead use cubeatt. Also need to write the files out with a specific filename structure. This filename structure will be used by the shell scripts to make csv files readable by OLAF.

This will be done for the topography, latitude, longitude, and photocube products. I give an example of the conversion for the photocubes below.

bigmap=KAR543
dir=/Users/JW/Dropbox/PDS_LDAP16/products/isis/Karpinskiy/orig/$bigmap


while read img ; do
cubeatt from=img-${img}/IoverF-${bigmap}.cub to=${dir}/${bigmap}_${img}_IoverF.cub+BandSequential
cubeatt from=img-${img}/rawP-$bigmap.cub to=${dir}/${bigmap}_${img}_a.cub+BandSequential
cubeatt from=img-${img}/$bigmap-localI.cub to=${dir}/${bigmap}_${img}_i.cub+BandSequential
cubeatt from=img-${img}/$bigmap-localE.cub to=${dir}/${bigmap}_${img}_e.cub+BandSequential
done < lastFourList_$bigmap

Next is to make the label. We use a shell script to generate a csv file, and then upload the cube and csv file into OLAF. To use the shell script, create a config.txt that looks similar to this...

<rand> KAR543$ cat config.txt 
OBJECT=Moon
BIGMAP=KAR543
AUTHLIST="Weirich, J.R."
OBSSYSBOOK="LROC"
TARGETNAME="Moon"
TARGETTYPE=Satellite
STARTTIME="N/A"
STOPTIME="N/A"
REFKEY="N/A"
VER=v1
SAMPLES=1015
LINES=1015
VERTPXSCALE=5
HOZPXSCALE=5
WEST=167.284134
EAST=167.880630
SOUTH=73.4460144
NORTH=73.6149368
SPHNAME=Moon_2000
AAXIS=1737400
BAXIS=1737400
CAXIS=1737400
LOCALDESC="This file is simple in a local reference frame and thus not in a map projection.  Within cart:local_georeference_information 
we provide the needed vectors such that a transformation could be applied to this file."
GEOINFO=""Vector" is the radius vector that extends from the center of the object to the middle vertex of the DTM. It is in Cartesian space with the first number being in the direction of x, the second number the direction of y, and the third number the direction of z. 0 West Longitude defines the positive x-axis, while north pole is the positive z-axis. 270 West Longitude defines the positive y-axis.
The plane of the DTM is determined by the local slope. This plane of the DTM is defined by two horizontal unit vectors (i.e. Ux and Uy) and a vertical unit vector (i.e. Uz). Each of these unit vectors are defined in the same Cartesian space as the radius vector, where again the first number is in the direction of x, the second number the direction of y, and the third number the direction of z.
 # Vector   -480.19241333007812        105.74529266357422        1663.2192382812500
 # Ux  -0.94595688581466675       0.21812622249126434      -0.23997193574905396
 # Uy  -0.22469174861907959      -0.97442990541458130        0.0000000000000000
 # Uz  -0.23383583128452301        5.3919713944196701E-002  0.97077983617782593
"

Vector and Ux, Uy, and Uz come from dumpMapHeaders.

Next is the run shell scripts. For the topography, lat, and lon cubes, you'll use productInfoOLAF-isis-unprojected-topocubes.sh with a command variable giving the directory location of config.txt. For the photometric cubes use productInfoOLAF-isis-unprojected-photocubes.sh with command variables of the directory location and a file with the name of all the images.

Once you have the cube in PDS format and a csv file, upload these to OLAF. Use the "2-D Array Images" Product, and upload everything.

Slope

* See topography

Albedo

* See topography

Sigmas

Run backplanesGMT (this gets lon/lat only file)

Take the SIGMAS.pgm and convert it to lon/lat/val files

   name=LEADEQ
   std2isis from=SIGMAS.pgm to=tmp.cub
   isis2ascii from=tmp.cub to=tmp1.txt
   sed '1,2d' tmp1.txt > tmp2.txt
    ~/bin/bustGrid tmp2.txt > $name-sig.txt
   paste $name-lonlat.txt $name-sig.txt > sigma.ll
   var=sigma.ll

Continue with topography script

Number of Images

run map_coverage_p_low
RESLIM should be 3x best maplet GSD
Be sure to exclude limb images using coverage_p.in

Convert coverage out to lon/lat/value

   name=LEADEQ
   bustGrid $name-cov.txt > tmp1.txt
   paste $name-lonlat.txt tmp1.txt > imgNum.ll
   var=imgNum.ll

* Continue with topography script

Best Image Resolution

run map_coverage_p_low
RESLIM should be high enough that there are no holes in coverage
Be sure to exclude limb images using coverage_p.in

Convert resolution out to lon/lat/value

   name=LEADEQ
   bustGrid $name-res.txt > tmp1.txt
   paste $name-lonlat.txt tmp1.txt > imgRes.ll
   var=imgRes.ll