Overview

Archiving SPC-focus products. There are two major categories (global and regional) with two subdivisions (topography and backplanes -- aka other data). The regional products are derived from bigmaps.

Tools

These tools are stored in GitHub as archivingTools.

• convertLatLon - Reads an X,Y,Z (or X,Y,Z,V) list and converts it into lat, lon, distance or lat, lon, value

• backplanesGMT - Reads a bigmap and generates a list of latlon, radius, albedo, slope. Use paste latlon radius > tmp for GMT.

• backplanesISIS - Reads a bigmap and generates 2D matrix for radius, albedo, slope, lat, lon

• phasei - Reads a bigmap and image (the SUMFILE) and generates 2D matrix for incidence, emission, phase angles

• coverage_p_low -- Reads all images and maplets. Creates a pgm that has every pixel marked for every image that it has, counter is only 1.

• map_coverage_p_low - Counts the number of images that cover each bigmap pixel. Creates MAPNM-cov.txt and MAPNM-cov.pgm

• bustGrid - Takes an ascii grid output and puts each entry on its own line. Needed to take grid and paste it with lat/lon to make grid

Grid to GeoTIF

var=radius.ll
new=big3
argR="-R280/287/-9/-3"
argR="-R0/360/-90/90"
argI="-I0.04"

date
echo "Blockmean"
gmt blockmean $var $argR $argI  > out

echo "sphereinterpolate"
cat out |  gmt sphinterpolate $argR $argI -Gtest.nc


echo "gdal translate"
gdal_translate -of GTiff -b 1 -a_srs support/60300.prj  NETCDF:test.nc $new.tif

Global

Global ICQ resolution to bin sizes [degrees]

Topography

• /bin/cp -f SHAPE.TXT tmp
  vi tmp
  1,%s/D/E/g
  ZZ
  convertLatLon tmp > radius.ll

Sigmas

• convertLatLon tmp 1 > sigma.ll 

• Use sigma.ll as the input to GMT

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 Resolution

• Run globalStat
• Use global-res as the input to GMT

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

Slope

Albedo

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

Regional

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

Topography

• Run backplanesGMT

  •   # Make the backplanes
      name=LEADEQ
      echo $name | ~/bin/backplanesGMT

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    

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

prj files can be found here

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

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

Photometric Data

This includes emission angle, incidence angle and phase angle.

Run phasei (located in ORExSPCsupport)

Run process similar to Number of Images (i.e. use "bustGrid" and "paste")

Running GMT on the PSI Hawk Cluster

To begin making projected geoTiffs for ArcGIS and ISIS, start with the txt files generated when making the photometric cubes. The scripts used produce the following txt files:

• IoverF-<bigmap>-<img>.txt (Ex. IoverF-LEADHL-N1563651588.txt)

• LEADHL-a.TXT
• LEADHL-e.TXT
• LEADHL-i.TXT

The above files need to be uploaded to the Hawk cluster, while keeping their parent directory structure.

Go to the photocube directory for each bigmap

rsync -hapvP --exclude=*.cub --exclude=*.prt N1* [email protected]:/home/jweirich/Tethys/geoTiffTethys/LEADHL/

Also copy the output from backplanesGMT and config.txt

Ex. config.txt

BIGMAP=LEADHL
ARGR=-R226.209259/261.399750/3.66536975/33.9945526
ARGI=-I801+n
DEGPX=0.0378641483770287
PROJECT=60300
BODY=Tethys
RUNALL=No
VERSION=1

The basic command being run is similar to the Topography step above:

echo "Blockmean"
gmt blockmean $var $argR $argI  > ${product}out
echo "sphereinterpolate"
cat ${product}out |  gmt sphinterpolate $argR $argI -G${product}.nc

... where $product is going to be IoverF, a, e, or i. Everything else is just a wrapper to deconflict all the files since we're running in parallel.

I have the GMT commands above in a script called makeaGeoBig.sh which reads most of the inputs from config.txt. makeGeoBig.sh is called by makeIoverFGeoBig.sh, makePhaseGeoBig.sh, makeEmissGeoBig.sh, and makeIncGeoBig.sh

The key lines in these scripts are: (Ex. below is for Incidence)

echo "getting ready for GMT ..."
~/bin/bustGrid $bigmap-i.TXT > bustedi-$bigmap-$img.txt
paste ../$bigmap-lonlat.txt bustedi-$bigmap-$img.txt > $var

echo "calling makeGeoTiffBig.sh"
sh ~/bin/makeGeoTiffBig.sh $var $product

$product is as before, while $var is IoverF.ll, phase.ll, emiss.ll, or inc.ll

The above four scripts are called by run4.sh, which is essentially a copy/paste script with different images. Note the use of & and "wait". & is to get multiple processors running on each machine, and "wait" is to keep the script running until all four processes complete.

The key lines from run4.sh look like this:

cd $img
sh ~/bin/makeIoverFGeoBig.sh $img &
sh ~/bin/makePhaseGeoBig.sh $img &
sh ~/bin/makeEmissGeoBig.sh $img &
sh ~/bin/makeIncGeoBig.sh $img &
cd ..

wait

Since there are 8 processors per Hawk node, we want to call 2 images per node. To do this we have a script that uses sbatch commands. Here's what one of those looks like:

#
#SBATCH -J Call1.sh              # Job name
#SBATCH -t 2-0:00 # time (D-HH:MM)
#SBATCH --partition=cpu
#SBATCH -o slurm.%N.%j.out # STDOUT
#SBATCH -e slurm.%N.%j.err # STDERR
 
sh ~/bin/run4.sh N1563651588 &
sh ~/bin/run4.sh N1563651648 &
 
wait

Rather than type each by hand, I instead have a script (makeScript.sh) to make a bunch of these:

# 29 Sep 2021 - John R. Weirich
# Script to call different photo scripts

imgList=$1

if [ -z $imgList ]; then
        echo "Please select a list of images"
        echo "Usage: <program> <Image List File>"
        exit
fi

list=`cat $imgList`

top=1
count=0

rm -f Call*.sh

for i in $list
do
 if [ $top = "1" ]
 then
  let count=$count+1
  out="Call$count.sh" 
  echo "#!/bin/bash" > $out
  echo "#" >> $out
  echo "#SBATCH -J $out              # Job name" >> $out
  echo "#SBATCH -t 2-0:00 # time (D-HH:MM)" >> $out
  echo "#SBATCH --partition=cpu" >> $out
  echo "#SBATCH -o slurm.%N.%j.out # STDOUT" >> $out
  echo "#SBATCH -e slurm.%N.%j.err # STDERR" >> $out
  echo " " >> $out
  echo "sh ~/bin/run4.sh $i &" >> $out

  top=2
 else
  echo "sh ~/bin/run4.sh $i &" >> $out
  echo " " >> $out
  echo "wait" >> $out

  top=1
 fi

done

if [ $top = "2" ]
then
 echo " " >> $out
 echo "wait" >> $out
fi

To fire everything off, you'll make a copy/paste that will look like this:

sbatch Call1.sh
sbatch Call2.sh
etc.

Once all those processes are finished, you'll have converted each *.txt file into a *.nc file. Now pull that back down to your local machine to run through GDAL.

rsync -hapvP --prune-empty-dirs --include="*/" --include="*.nc" --exclude="*" [email protected]:/home/jweirich/Tethys/geoTiffTethys/LEADHL/ .

Once downloaded, run

sh ../bin/wrapperAfterClusterImage.sh imgListLEADHL 

wrapperAfterClusterImages.sh is this:

# 30 Sep 2021 - John R. Weirich
# Wrapper to turn all the *.nc into geoTIFFS

imgList=$1

if [ -z $imgList ]; then
        echo "Please select a list of images"
        echo "Usage: <program> <Image List File>"
        exit
fi

list=`cat $imgList`

for i in $list
do
 cd $i
 sh /usr/local/spc/bin/afterClusterImage.sh IoverF $i
 sh /usr/local/spc/bin/afterClusterImage.sh a $i
 sh /usr/local/spc/bin/afterClusterImage.sh e $i
 sh /usr/local/spc/bin/afterClusterImage.sh i $i
 cd ..
done

And afterClusterImage.sh is this:

# 28 Sep 2021 - John R. Weirich
# Make the various geoTiffs using GDAL
# Run from geoTiff[directory]/
# 30 Sep 2021: Began modifying, then undid changes (hopefully they are undone correctly!)

# Usage in geoTiff directory : sh /usr/local/spc/bin/afterClusterBigmap.sh <product/type name> <image>

type=$1
img=$2

if [ -z $type ]; then
        echo "Please select the type"
        echo "Usage: <program> <type> <Image>"
        exit
fi

if [ -z $img ]; then
        echo "Please select the image"
        echo "Usage: <program> <type> <Image>"
        exit
fi



 if [ ! -e ../config.txt ]
 then
  echo "Make config.txt!"
  exit
 fi

 source ../config.txt

 if [ ! -e ./${type}.nc ]
 then
  echo "Make ${type}.nc!"
  exit
 fi



bigmap=$BIGMAP
argR=$ARGR
argI=$ARGI
degPx=$DEGPX
proj=$PROJECT
body=$BODY
runAll=$RUNALL
ver=$VERSION

echo "$bigmap"
echo "$argR"
echo "$argI"
echo "$degPx"
echo "$proj"
echo "$body"
echo "$runAll"


# Make normal GeoTiff

echo "gdal translate"

gdal_translate -of GTiff -b 1 -a_srs ../support/$proj.prj  NETCDF:${type}.nc ${bigmap}-${img}-${type}-v${ver}.tif

#Make GeoTiff readable by ISIS
gdal_translate -of ISIS3 -tr $degPx $degPx -r bilinear -b 1 -co TARGET_NAME=$body -co DATA_LOCATION=GEOTIFF -a_srs ../support/$proj.prj NETCDF:${type}.nc ${bigmap}-${img}-${type}-ISIS-v${ver}.lbl

Make Pretties

gmt begin GMT_cont
gmt set GMT_THEME cookbook
gmt grdcontour test.nc
gmt end show
#gmt grdcontour test.nc -C10 -A50

gmt begin GMT_img
gmt set GMT_THEME cookbook
#gmt makecpt -Crainbow
gmt grdimage test.nc  -JM6i -B -BWSnE
gmt colorbar -DJTC -Bxa -By+lm
gmt end show


gmt begin GMT_img
gmt makecpt -Crainbow
gmt set GMT_THEME cookbook
gmt grdimage test.nc  $argR  -JM6i -B -BWSnE
#gmt colorbar -DJTC -I0.4 -Bxa -By+lm
gmt colorbar -DJTC -Bxa -By+lm
gmt end show

Setup

Use Dropbox spcShare directory

• cd /opt/local/spc/
  relink ~/Dropbox/spcShare/Tethys .

Then in the working directory

• relink.sh /opt/local/spc/Tethys/support .
  relink.sh /opt/local/spc/Tethys/mapConfig .