Gridding with GMT

Online Resources

https://youtu.be/D7do_whEkI0
https://developers.google.com/kml/articles/raster
https://www.geologie.ens.fr/~ecalais/teaching/maps-and-graphs-with-gmt/gmt_6.pdf ← just great GMT examples
https://www.aviso.altimetry.fr/fileadmin/documents/OSTST/2009/oral/Scharroo_GMT_GE.pdf another good example

GMT OPenMP

FYI, OpenMP on GMT can be enabled by compiling from scratch, but it isn't in the nearneighbor module, so it doesn't speed up what I am using

Rename Files

#!/bin/bash
 
for file in *.qlv2; do
    mv "$file" "$(basename -s .qlv2 "$file").ql"
done

From HDF5 Directly

Just a scratch pad of gridding from a quicklook file

gmt info lvis_20210806_155400_00000316.geolocresid.h5?/BOUNCEPOINT_LONGITUDE_FR3RX3_A,/BOUNCEPOINT_LATITUDE_FR3RX3_A,/BOUNCEPOINT_ALTITUDE_FR3RX3_A -I0.01

region=$(gmt info lvis_20210806_155400_00000316.geolocresid.h5?/BOUNCEPOINT_LONGITUDE_FR3RX3_A,/BOUNCEPOINT_LATITUDE_FR3RX3_A,/BOUNCEPOINT_ALTITUDE_FR3RX3_A -I0.01)

gmt convert lvis_20210806_155400_00000316.geolocresid.h5?/BOUNCEPOINT_LONGITUDE_FR3RX3_A,/BOUNCEPOINT_LATITUDE_FR3RX3_A,/BOUNCEPOINT_ALTITUDE_FR3RX3_A > out.bin
gmt nearneighbor out.bin -Gout.grd -I0.10c=/0.10c= -S5.0c/5.0c -E-9999 -V -R283.01/283.12/41.38/41.45
gmt nearneighbor out.bin -Gout.grd -I0.10c=/0.10c= -S5.0c/5.0c -ENaN -V -R283.01/283.12/41.38/41.45
gmt nearneighbor out.bin -Gout.grd -I0.10c=/0.10c= -S5.0c/5.0c -ENaN -V $region
gmt makecpt -T-0/600 -Z > rainbow.cpt
gmt grdimage -JI15c out.grd > out.ps
gmt grdimage -Crainbow.cpt -JI15c out.grd > out.ps
gmt grdimage -Crainbow.cpt -JM6i out.grd > out.ps
gmt grd2cpt out.grd -E255 -V > color.cpt
gmt grdimage -Ccolor.cpt -Q -JM6i out.grd > out.ps
gs out.ps

gmt grdimage -Ccolor.cpt -Q -JM6i out.grd -Aout.png

GMT HDF5 Directly

Looks like once you grid it, the grd file looses the geospatial reference? Maybe need to add a projection or geotiff kind of headers to the file to retain WHERE these pixels are, because gmt info now doesn't show a reasonable 'region' for the grid. So how could you combine grids? Or how are we combining these data sets?

Extract the region information from the raw data

region=$(gmt info lvis_20210806_155400_00000316.geoloc.h5?/BOUNCEPOINT_LONGITUDE_FR3RX3_A,/BOUNCEPOINT_LATITUDE_FR3RX3_A,/BOUNCEPOINT_ALTITUDE_FR3RX3_A -I0.01)

You can grid directly from the H5 file:

gmt nearneighbor lvis_20210806_155400_00000316.geoloc.h5?/BOUNCEPOINT_LONGITUDE_FR3RX3_A,/BOUNCEPOINT_LATITUDE_FR3RX3_A,/BOUNCEPOINT_ALTITUDE_FR3RX3_A -Gout.grd -I0.10c=/0.10c= -S5.0c/5.0c -ENaN -V $region

OR convert the raw HDF5 data into GMT binary

gmt convert lvis_20210806_155400_00000316.geoloc.h5?/BOUNCEPOINT_LONGITUDE_FR3RX3_A,/BOUNCEPOINT_LATITUDE_FR3RX3_A,/BOUNCEPOINT_ALTITUDE_FR3RX3_A > out.bin

And then grid the data using the region variable

gmt nearneighbor out.bin -Gout.grd -I0.10c=/0.10c= -S5.0c/5.0c -ENaN -V $region

Create a color table using the grid data

gmt grd2cpt out.grd -E255 > color.cpt

Or grab and use a standard one from online

wget http://soliton.vm.bytemark.co.uk/pub/cpt-city/mby/mby.cpt
gmt grdimage -Cmby.cpt -Q -JM6i out.grd -Aout.png

Create image from the grid

gmt grdimage -Ccolor.cpt -Q -JM6i out.grd -Aout.png

Note: running this on one file is about 2 minutes, running on 4 files took about 13 minutes. Since the grid is getting larger, searching gets exponentially larger.

GMT Combining Grids

So, you want to take a grid from each file and combine it all into one big file! Well, follow this recipe. (this is pretty fast)

Right now, this has a line at the start/end of each grid. May need to

blend things?
make the regions bigger?
google some more (-8

Found this here: https://forum.generic-mapping-tools.org/t/grid-blending-of-gridded-files/119

gmt cut
gmt sample
gmt blend
gmt grd2xyz
5 gmt surface
gmt sample.
these 6 steps helped me to successfully deal with the issue. am grateful for everyone contribution.

Create an ASCII file that lists each grid, region and weight for example, I made one called combined.txt for 4 example grids

Edit this file
```
nano combined.txt
```

Here are the values I put in (but getting the region for each file)

lvis_20210806_155400_00000316.grd -R283.02/283.11/41.38/41.45 1
lvis_20210806_155430_00000317.grd -R283.08/283.18/41.42/41.48 1
lvis_20210806_155500_00000318.grd -R283.15/283.25/41.46/41.52 1
lvis_20210806_155530_00000319.grd -R283.22/283.31/41.49/41.56 1

Combine these grids with grdblend

gmt grdblend combined.txt -Gblend.grd -R283.00/283.32/41.38/41.56 -I0.10c=/0.10c= -V

Generate an image from this new super grid

gmt grdimage -Cmby.cpt -Q -JM6i blend.grd -Aout.png

Gridding with GDAL

(shouldn't this be able to read the HDF5 file directly?)

Convert the data to an ascii file first

h5dump -d /BOUNCEPOINT_LONGITUDE_FR3RX3_A \
   -o h5lon.asci -y -w 1 --format=%.8f lvis_20210806_155400_00000316.geoloc.h5

h5dump -d /BOUNCEPOINT_LATITUDE_FR3RX3_A \
   -o h5lat.asci -y -w 1 --format=%.8f lvis_20210806_155400_00000316.geoloc.h5

h5dump -d /BOUNCEPOINT_ALTITUDE_FR3RX3_A \
   -o h5alt.asci -y -w 1 --format=%.3f lvis_20210806_155400_00000316.geoloc.h5

# make these files into columns ascii
paste h5lon.asci h5lat.asci h5alt.asci > dem.csv

Make a header file so gdal can read the csv file

Edit a file
```
nano dem.vrt
```

It should contain this

<OGRVRTDataSource>
    <OGRVRTLayer name="dem">
        <SrcDataSource>dem.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <GeometryField encoding="PointFromColumns" x="field_1" y="field_2" z="field_3"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

Grid with GDAL 0.0001 is about a 100m grid (at equator anyway, longitude shrinks as you go north, maybe should change for latitude)

time gdal_grid -a invdist:power=2.0:smoothing=1.0:radius1=0.0001:radius2=0.0001:nodata=NaN -outsize 420 420 -of GTiff -ot Float64 -l dem dem.vrt dem.tiff --config GDAL_NUM_THREADS ALL_CPUS

Future

Could convert each grid back to xyz with https://gdal.org/programs/gdal2xyz.html to make a big image with all the data?
I think we want to make an image for each file and just use them as tiles to build up a big KML which I think is what Sarah's python is doing.

Scratch Pad

Divide the X size by sin(latitude) to try and get equal spacing in grid.

This is pretty slow

time gdal_grid -a invdist:power=2.0:smoothing=1.0:radius1=0.0001:radius2=0.0001:nodata=NaN \
   -outsize 4200 2700 \
   -of GTiff \
   -ot Float64 \
   -l dem \
   dem.vrt \
   dem.tiff \
   -a_srs EPSG:4326 \
   --config GDAL_NUM_THREADS ALL_CPUS

This is faster (still see diamonds in the resulting grid, something with grid code?)

time gdal_grid -a invdistnn:power=2.0:smoothing=0.0:radius=0.0005:nodata=NaN \
   -outsize 4200 2700 \
   -of GTiff \
   -ot Float64 \
   -l dem \
   dem.vrt \
   dem.tiff \
   -a_srs EPSG:4326 \
   --config GDAL_NUM_THREADS ALL_CPUS

should be able to fix our grid spacing and let it figure out the number of elements in x.y

-tr <xres> <yres>
# but you need to set
-txe <xmin> <xmax>
-tye <xmin> <xmax>

Build up a tiles from a bunch of images

https://gis.stackexchange.com/questions/140877/create-web-map-tiles-from-large-images

gdalbuildvrt -o merged.vrt file1.jp2 file2.jp2 ....
gdal2tiles.py merged.vrt output_folder/

This could work! kind of does (too many options with gdal, there has to be a righter way)

gdal_merge.py -n nan -o all.tif -of gtiff lvis_*.tif

david.rabine

Table of Contents