Archiving SPC

Global ICQ resolution to bin sizes [degrees]

ICQ Q Value	Vertices	Bin Size [deg]
512	1,579,014	0.041
256	396,294	0.164
128	99,846	0.649
64	25,350	2.556
32	6,534	9.917

Topography

/bin/cp -f SHAPE.TXT tmp
vi tmp
1,%s/D/E/g  # Remove first line of shape if it has the Q size, replace D with E
ZZ
convertLatLon tmp > radius.ll # If you are not using the "1" option, "tmp" needs to be 3 columns, instead of 4.

Commands developed while talking with Trent Hare 22 Aug 2022.

GMT pads with half a pixel, so for a product with 1 deg bins, need to trim 0.5 deg. 0.3 deg bins would trim 0.15 deg. Example below is for a product for Phoebe with 1 deg bins where the output from covertLatLon is radius128.ll

var=radius128.ll
argR=-R0.5/359.5/-89.5/89.5 # trim per above; gives a warning when there are more than two decimal places, but the end product is still good. 
argI=-I1.0 # 1 deg bins

gmt blockmean $var $argR $argI  > out # This step may not be needed, but keep for now
cat out |  gmt sphinterpolate $argR $argI -GtestRadius128Trim.nc

Make a cube with Long from 0 to 360. Phoebe is 106.5 km in radius. Circumference divided by 2 gives the distance in meters from center to western/eastern end. Circumference divided by 4 gives distance in meters from equator to pole. -a_ullr is upper left pixel in -long lat and lower right pixel in lon -lat. BOUNDING_DEGREES is what is put in the label (remember that GMT pads by half a pixel, so we are not generating data when there is none). May need to contact Trent Hare to get the correct prj file

gdal_translate -of ISIS3 -co TARGET_NAME=Phoebe -co FORCE_360=True -a_srs support/60915FromTrent.prj -a_ullr -334579.61760731 167289.808803656 334579.61760731 -167289.808803656 -co BOUNDING_DEGREES=0,-90,360,90 NETCDF:testRadius128Trim.nc out_testRadius128Trim.cub

The above command should have GDAL write the label correctly, but apparently there is a bug in the code. My version is more recent than Trent's version, so it shouldn't be a version problem. So now we use ISIS to make a correct ISIS label.

First set up the MAP parameters. Note that for 0 to 360 clon=180, this is because ISIS isn't using clon (it's a misnomer), but is really setting x=0. -180 to 180 would use clon=0. Resolution is pixels/deg, so for 0.3 deg bins you would enter 3.3333 here.

maptemplate map=smp180.map projection=simplecylindrical clon=180 targopt=user targetname=Phoebe eqradius=106500 polradius=106500 resopt=ppd resolution=1

The below will overwrite the given file. Note this is mapping a single pixel (sample and line are always 0.5, 0.5 since that refers to the top left corner of the top left pixel, and maps to lat 90 lon 0.)

maplab from=out_testRadius128Trim.cub map=smp180.map sample=0.5 line=0.5 coordinates=latlon lat=90 lon=0

The above gives a working cube, now to turn it into a geotiff. We use gdal_translate again. To make the label a bit cleaner, could use -tr 1858.77565 1858.77565 (which comes from gdalinfo on the cube and truncating the Pixel Size).

gdal_translate -of PDS4 -co IMAGE_FORMAT=geotiff out_testRadius128Trim.cub out.lbl

Follow up with using vi to edit the <cart:pixel_scale_x unit="pixel/deg">0.999999999999978906</cart:pixel_scale_x> value to be 1 so the label looks nicer.

Sigmas

```
convertLatLon tmp 1 > sigma.ll 
```

Use sigma.ll as the input to GMT

After you have the cube, change values of "-1" (which represent 1 or 0 maplets) to "9999"

fx f1=phoebe_sigma_tmp.cub to=phoebe_sigma_c.cub equation="(9999-f1)*(f1<0) + f1" // Use GUI for this since command line will fail

Number of Images

Run globalStat

Use global-cov.ll as the input to GMT

var=global-cov.ll
new=tCov
argR="-R0/360/-90/90"
argI="-I1"

Best Image Resolution

Run globalStat

Use global-res as the input to GMT

var=global-res.ll
new=tRes
argR="-R0/360/-90/90"
argI="-I1"

Best Maplet Resolution

Use coverage and choose min/max ranges that bracket all the different maplets resolutions. So if you have 8 different resolutions, you'll have 8 different pgm files. Note that while coverage/SPC works in West longitude, when writing out to an image the presentation doesn't matter if it's West or East longitude since they will look the same. West longitude just starts with 0 in the East and increases West, while East Longitude starts with 0 in the West and increases to the East.

For each pgm file, use USGS ISIS to turn it into a cube where each pixel has the value of the resolution in meters, and the rest of the pixels have a 9999 meters.

Change pgm to cube
Turn all DN values of 0 to (9999-resolution), and all non-zero DN to 0. The "==" logic comparison returns 1 if true and 0 if false
Add maplet resolution to all DN values

std2isis fr=Phoebe0p00-0p13Maplets.pgm to=Phoebe0p00-0p13Maplets.cub
fx f1=Phoebe0p00-0p13Maplets.cub to=Phoebe0p00-0p13flip.cub equation=(9999-125)*(f1==0) // Here the maplet resolution is 125 ; command line will error out, need to use GUI
fx f1=Phoebe0p00-0p13flip.cub to=Phoebe0p00-0p13res.cub equation=f1+125
...
std2isis fr=Phoebe1p10-2p00Maplets.pgm to=Phoebe1p10-2p00Maplets.cub
fx f1=Phoebe1p10-2p00Maplets.cub to=Phoebe1p10-2p00flip.cub equation=(9999-1500)*(f1==0) // Here the maplet resolution is 1500
fx f1=Phoebe1p10-2p00flip.cub to=Phoebe1p10-2p00res.cub equation=f1+1500

Now take the minimum value from all the different cubes

fx f1=Phoebe0p00-0p13res.cub f2=Phoebe0p13-0p29res.cub to=Phoebe0p00-0p29.cub equation=min(f1,f2) // ; command line will error out, need to use GUI
...
fx f1=Phoebe0p00-1p10.cub f2=Phoebe1p10-2p00res.cub to=Phoebe0p00-2p00.cub equation=min(f1,f2)

Slope

Albedo

/bin/cp -f SHAPEA.TXT tmp
vi tmp
1,$s/D/E/g
ZZ
convertLatLon tmp 1 > albedo.ll