david.rabine

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== 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 =====

<code bash>
#!/bin/bash

for file in *.qlv2; do
    mv "$file" "$(basename -s .qlv2 "$file").ql"
done
</code>

===== From HDF5 Directly =====

Just a scratch pad of gridding from a quicklook file

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

</code>

===== GMT HDF5 Directly =====

  * FIXME 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 <code>
region=$(gmt info lvis_20210806_155400_00000316.geoloc.h5?/BOUNCEPOINT_LONGITUDE_FR3RX3_A,/BOUNCEPOINT_LATITUDE_FR3RX3_A,/BOUNCEPOINT_ALTITUDE_FR3RX3_A -I0.01)
</code>
  - You can grid directly from the H5 file: <code>
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
</code>
    - OR convert the raw HDF5 data into GMT binary <code>
gmt convert lvis_20210806_155400_00000316.geoloc.h5?/BOUNCEPOINT_LONGITUDE_FR3RX3_A,/BOUNCEPOINT_LATITUDE_FR3RX3_A,/BOUNCEPOINT_ALTITUDE_FR3RX3_A > out.bin
</code>
    - And then grid the data using the region variable <code>
gmt nearneighbor out.bin -Gout.grd -I0.10c=/0.10c= -S5.0c/5.0c -ENaN -V $region
</code>
  - Create a color table using the grid data <code>
gmt grd2cpt out.grd -E255 > color.cpt
</code>
    - Or grab and use a standard one from online <code>
wget http://soliton.vm.bytemark.co.uk/pub/cpt-city/mby/mby.cpt
gmt grdimage -Cmby.cpt -Q -JM6i out.grd -Aout.png
</code>
  - Create image from the grid <code>
gmt grdimage -Ccolor.cpt -Q -JM6i out.grd -Aout.png
</code>

{{:instrument:development:lvis:2020:rtproducts:gridding:out.png?400|}}

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

FIXME 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]]
<code>
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.
</code>

  - 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 <code>
nano combined.txt
</code>
    - Here are the values I put in (but getting the region for each file)<code>
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
</code>
  - Combine these grids with **grdblend** <code>
gmt grdblend combined.txt -Gblend.grd -R283.00/283.32/41.38/41.56 -I0.10c=/0.10c= -V
</code>
  - Generate an image from this new super grid <code>
gmt grdimage -Cmby.cpt -Q -JM6i blend.grd -Aout.png
</code>

{{:instrument:development:lvis:2020:rtproducts:gridding:blend4_try1.png?400|}}
===== Gridding with GDAL =====

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

  - Convert the data to an ascii file first <code>
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
</code>
  - Make a header file so **gdal** can read the csv file
    - Edit a file <code>
nano dem.vrt
</code>
    - It should contain this <code>
<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>
</code>
  - Grid with GDAL 0.0001 is about a 100m grid (at equator anyway, longitude shrinks as you go north, maybe should change for latitude) <code>
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
</code>

{{:instrument:development:lvis:2020:rtproducts:gridding:qgis_dem_load.png?400|}}

==== 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
<code>
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
</code>

This is faster (still see diamonds in the resulting grid, something with grid code?)
<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
</code>

FIXME should be able to fix our grid spacing and let it figure out the number of elements in x.y
<code>
-tr <xres> <yres>
# but you need to set
-txe <xmin> <xmax>
-tye <xmin> <xmax>
</code>

Build up a tiles from a bunch of images
  * [[https://gis.stackexchange.com/questions/140877/create-web-map-tiles-from-large-images]]
<code>
gdalbuildvrt -o merged.vrt file1.jp2 file2.jp2 ....
gdal2tiles.py merged.vrt output_folder/
</code>

This could work!  kind of does  (too many options with gdal, there has to be a righter way)
<code>
gdal_merge.py -n nan -o all.tif -of gtiff lvis_*.tif
</code>