====== Gridding with GMT ======
===== RTP High Level =====
==== Scripts Running ====
- rtp_raw2ql.sh - every 20 seconds, look for a new md5sum
- rtp_grid.sh - one for each thread so multiple grids can be happening at a time
- rtp_dem.sh - every 15 minutes build up a summary
###### VEGAS config settings ######
geolocation_only
obs_bias -137.807 0.0
instrument_roll_bias 1.9265100000000000 0.0
instrument_pitch_bias -1.4537700000000000 0.0
instrument_yaw_bias 0.359650000000000 0.0
scan_angle_roll_scale -1.0018271200000000 0.0
scan_angle_pitch_scale 0.9970223500000000 0.0
roll_pitch_alignment_z 90.0000000000000000 0.0
time_tag_bias_pitch_to_obs -0.0015568100000000 0.0
time_tag_bias_roll_to_obs -0.0006150700000000 0.0
fiber_speed_factor 1.0
use_channel_a
use_channel_b
use_channel_c
use_range_frrx
measurements_per_batch 30000
disable_gphyscorr_all
ATT_FILE "/Volumes/XSan4/david/rtp/traj59432_rtp_192.168.0.63_quicklook.att"
TRAJ_FILE "/Volumes/XSan4/david/rtp/traj59432_rtp_192.168.0.63_quicklook.traj"
MEAS_FILE "/Users/drabine/rtp/src/lvis_20210806_155400_00000316.meas.h5"
GEOLOC_RESID_FILE "/Users/drabine/rtp/src/lvis_20210806_155400_00000316.geoloc.h5"
===== 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
===== 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 =====
* 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
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
{{: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]]
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
{{: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
h5dump -d /BOUNCEPOINT_LONGITUDE_FR3RX3_A \
-o h5lon.asci -y -w 1 -m %.8f lvis_20210806_155400_00000316.geoloc.h5
h5dump -d /BOUNCEPOINT_LATITUDE_FR3RX3_A \
-o h5lat.asci -y -w 1 -m %.8f lvis_20210806_155400_00000316.geoloc.h5
h5dump -d /BOUNCEPOINT_ALTITUDE_FR3RX3_A \
-o h5alt.asci -y -w 1 -m %.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
dem.csv
wkbPoint
- 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
{{: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
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
FIXME should be able to fix our grid spacing and let it figure out the number of elements in x.y
-tr
# but you need to set
-txe
-tye
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