Differences between revisions 76 and 77
Revision 76 as of 2022-03-10 15:52:30
Size: 24442
Editor: JohnWeirich
Comment:
Revision 77 as of 2022-03-10 16:02:02
Size: 24689
Editor: JohnWeirich
Comment:
Deletions are marked like this. Additions are marked like this.
Line 672: Line 672:

The below was written prior to utilizing gdalwarp to make the geotiffs. See "Grid to GeoTIF" section for details, but if you use gdalwarp to make the geotiff you can also output a sample xml label that has most (all?) of the information below.

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

  • globalStat - Renamed version of coverage_p_low that does more (JRW not sure 29 Nov 2021)? Output test version of coverage and max res of images, as well as png of coverage


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

Note: gdal_translate doesn't define everything that is needed because the test.nc projection is not defined. Use the following instead. Below is an example for the Moon (Luna), where several important objectives are accomplished.

  1. The generated GeoTIFF is ready for ingestion into a PDS archive because it will use BandSequential format, and it will generate an example xml label which will have many hard-to-determine details like the header byte offset and the data type (i.e. LSB vs MSB).

  2. Trent Hare at USGS (personal communication on 02 Mar 2022) says in the planetary community usually uses square pixels, even though GeoTIFFs can handle rectangular pixels. Note that with the projections definitions (point 3.) the units of -tr will be meters. So -tr 1 1 will be 1 meter square vertices.
  3. gdalwarp defines the projection of the source file (i.e. -s_srs) and output/target file (i.e. -t_srs). I sent Trent Hare the first 10 lines of the input file to GMT and the tmp.nc and on 09 Mar 2022 he determined the source projection is a Cartographic system in degrees, so we can use 30100.prj. The example for 30110.prj is an Equidistant Cylindrical projection.

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


Also note that you can check the details of a GeoTIFF with gdalinfo and gdalsrsinfo.

Global

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

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

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

* 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* jweirich@hawk.psi.edu:/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:

#\!/bin/bash
#
#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="*" jweirich@hawk.psi.edu:/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

Note: See "Grid to GeoTIF" for a better program than gdal_translate, since gdal_translate doesn't always write the projection.


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 .

How to make DSK from ICQ and MAP

MAP to DSK

Pathway to make a DSK. Convert MAP to OBJ using AltWG tools, then convert OBJ to DSK using SPICE tools.

Example using TRALHL.MAP, to get obj (note: make sure you use the --local)

Maplet2FITS TRALHL.MAP tmp.plt
FITS2OBJ --local tmp.plt TRALHL.obj

Now turn OBJ into DSK (note: lbl file is listed below)

<rand> makeDSK$ /usr/local/toolkit/mice/exe/mkdsk 
 
MKDSK Program; Ver. 2.0.0, 28-FEB-2017; Toolkit Ver. N0066
 
SETUP FILE NAME> TRALHL-obj.lbl
Reading plate model input file...
...Done reading plate model input file.
 
Generating Spatial Index...
Segregating and closing DSK file...
DSK file was created.
 
All done.

Example of LBL file. Refer to https://naif.jpl.nasa.gov/pub/naif/utilities/MacIntel_OSX_64bit/mkdsk.ug for more details.

<rand> makeDSK$ cat TRALHL-obj.lbl 
\begindata
 
      INPUT_SHAPE_FILE    = 'TRALHL.obj'
      OUTPUT_DSK_FILE     = 'TRALHL.bds'
      SURFACE_NAME        = 'TRALHL Tethys'                 ### This can be user defined; see NAIF_SURFACE_NAME below
      CENTER_NAME         = 'TETHYS'                        ### Must be spice compatible name
      REF_FRAME_NAME      = 'IAU_TETHYS'                    ### Must be spice compatible frame
      START_TIME          = '1950-JAN-1/00:00:00'
      STOP_TIME           = '2050-JAN-1/00:00:00'
      DATA_CLASS          = 1                               ### Data class 1 is for shapes with a single radii for each lat/lon
      INPUT_DATA_UNITS    = ( 'ANGLES    = DEGREES'
                              'DISTANCES = KILOMETERS' )
      COORDINATE_SYSTEM   = 'LATITUDINAL'                   ### Haven't yet experimented with other systems
      DATA_TYPE           = 2                                ### As of 22 Nov '21, all DATA TYPES are 2
      PLATE_TYPE          = 3                               ### For an OBJ input use a Plate Type of 3

      KERNELS_TO_LOAD     = ( 'naif0012.tls' )              ### Not sure if this is needed


      NAIF_SURFACE_NAME   += 'TRALHL Tethys'                ### Here is where you define user surface
      NAIF_SURFACE_CODE   += 1                              ### The above SPICE link doesn't explain this code well; 1 seems to work
      NAIF_SURFACE_BODY   += 603                            ### This is the SPICE code for Tethys
 

      MINIMUM_LATITUDE    = 11.715                          ### I've played around changing the min/max lat/lon and it doesn't seem to matter much
      MAXIMUM_LATITUDE    = 29.104                          ### At least it doesn't seem to change anything once converted back to an OBJ
      MINIMUM_LONGITUDE   = 102.452
      MAXIMUM_LONGITUDE   = 122.322

      \begintext

ICQ to DSK

Pathway to make a DSK. Convert ICQ to OBJ using AltWG tools, then convert OBJ to DSK using SPICE tools.

Example using Tethys.txt, to get obj

ICQ2PLT Tethys.txt tmp.plt
PLT2OBJ tmp.plt Tethys.obj

To get OBJ to DSK, use the same technique above, except the min/max lat/lon should be -90 to 90 and -180 to 180.

DSK to OBJ

-dsk gives the input dsk, -text gives the output file, -format vertex-facet determines the output format to be in an obj

/usr/local/toolkit/mice/exe/dskexp -dsk TRALHL.bds -text TRALHL_fromDSK.obj -format vertex-facet

Note: To read the obj generated by dskexp in Meshlab, first trim the spaces in front of the "v"'s and "f"'s of the output file.

The SPICE commands DSKBRIEF and COMMNT may also be useful to check the dsk generated.

Notes on extracting information from geotiff. i.e. Entries for OLAF

The below was written prior to utilizing gdalwarp to make the geotiffs. See "Grid to GeoTIF" section for details, but if you use gdalwarp to make the geotiff you can also output a sample xml label that has most (all?) of the information below.

To get the pixel size of a geotiff, use gdalinfo. Command and output of a rectangular pixel size are shown below. Number of pixels is shown by the "Size is 501, 100" where 501 is px and 100 is line. Pixel size (in deg) is given by "Pixel Size = (0.036126000000000,-0.181000000000000)". Note you can also get this number by using the "Corner Coordinates:" of (292.297 - 274.198) / 501 = 0.036126, or (3.115 - -14.984) / 100 = 0.181. GDAL outputs a negative sign in front of the 0.181, I am not sure why.

<rand> test$ gdalinfo -stats LEADEQ-radius-v1.tif
Driver: GTiff/GeoTIFF
Files: LEADEQ-radius-v1.tif
       LEADEQ-radius-v1.tif.aux.xml
Size is 501, 100
Coordinate System is:
GEOGCRS["Tethys 2000",
    DATUM["D_Tethys_2000",
        ELLIPSOID["Tethys_2000_IAU_IAG",535600,54.6530612244898,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433,
            ID["EPSG",9122]]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]
Data axis to CRS axis mapping: 2,1
Origin = (274.197937000000024,3.115500000000000)
Pixel Size = (0.036126000000000,-0.181000000000000)
Metadata:
  AREA_OR_POINT=Area
  lat#actual_range={-14.894,3.025}
  lat#axis=Y
  lat#long_name=latitude
  lat#standard_name=latitude
  lat#units=degrees_north
  lon#actual_range={274.216,292.279}
  lon#axis=X
  lon#long_name=longitude
  lon#standard_name=longitude
  lon#units=degrees_east
  NC_GLOBAL#Conventions=CF-1.7
  NC_GLOBAL#GMT_version=6.2.0 [64-bit]
  NC_GLOBAL#history=sphinterpolate -R274.216/292.279/-14.894/3.025 -I501+n/100+n -Gtmp.nc
  z#actual_range={525616.3125,530976.5625}
  z#long_name=z
  z#_FillValue=nan
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (     274.198,       3.115) (274d11'52.57"E,  3d 6'55.80"N)
Lower Left  (     274.198,     -14.984) (274d11'52.57"E, 14d59' 4.20"S)
Upper Right (     292.297,       3.115) (292d17'49.43"E,  3d 6'55.80"N)
Lower Right (     292.297,     -14.984) (292d17'49.43"E, 14d59' 4.20"S)
Center      (     283.248,      -5.934) (283d14'51.00"E,  5d56' 4.20"S)
Band 1 Block=501x4 Type=Float32, ColorInterp=Gray
  Minimum=525616.312, Maximum=530976.562, Mean=528625.394, StdDev=890.336
  NoData Value=nan
  Metadata:
    actual_range={525616.3125,530976.5625}
    long_name=z
    NETCDF_VARNAME=z
    STATISTICS_MAXIMUM=530976.5625
    STATISTICS_MEAN=528625.39353917
    STATISTICS_MINIMUM=525616.3125
    STATISTICS_STDDEV=890.33610412375
    STATISTICS_VALID_PERCENT=100
    _FillValue=nan
  • Bits per Pixel is 32
  • Sample Direction is RIGHT
  • Line Direction is DOWN
  • Axis Order is Last_Index_fastest
  • Data Type is IEEE_REAL_LSB (or at least the ones I make with GDAL are LSB)
  • Header size can be found using Hex Fiend. Notice the red line a bit after the word "Greenwich". That is where you need to click, and then the header size is given at the bottom of the window. Click to see larger image if needed.

attachment:HexFiend2.png

  • To get Vertical Pixel scale take the second Pixel Size value (0.181000000000000), and radius of the object. V Px Scale (meters) = 2*PI*(Radius in meters)/360*(0.181000000000000)
  • To get Horizontal Pixel scale take the first Pixel Size value (0.036126000000000), radius of the object, and Latitude. H Px Scale (meters) = 2*PI(Radius in meters)*cos(Lat in radians)/360*(0.036126000000000)

Notes on extracting information from cube. i.e. Entries for OLAF

Most data can be gleaned from open the cube in the program "vi". Note that to get the cube into BandSequential you may need to run the ISIS program "cubeatt from=<filename1>.cub to=<filename2>.cub+BandSequential"

Object = IsisCube
  Object = Core
    StartByte = 65537
    Format    = BandSequential

    Group = Dimensions
      Samples = 1821
      Lines   = 1821
      Bands   = 1
    End_Group

    Group = Pixels
      Type       = Real
      ByteOrder  = Lsb
      Base       = 0.0
      Multiplier = 1.0
    End_Group
  End_Object
End_Object

Object = Label
  Bytes = 65536
End_Object

Object = History
  Name      = IsisCube
  StartByte = 13329701
  Bytes     = 1069
End_Object
End
  • Bits per Pixel is 32
  • Sample Direction is RIGHT
  • Line Direction is DOWN
  • Axis Order is Last_Index_fastest
  • Data Type is IEEE_REAL_LSB (or at least the ones I make with GDAL are LSB)
  • Header size is one less than the StartByte, so in the above case the StartByte is 65536.

Archiving SPC (last edited 2022-03-10 16:02:02 by JohnWeirich)