code:work:lvisf:2021:rtp:start
This is an old revision of the document!
Table of Contents
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://www.geologie.ens.fr/~ecalais/teaching/maps-and-graphs-with-gmt/gmt_6.pdf ← just great GMT examples
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
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
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
code/work/lvisf/2021/rtp/start.1640033439.txt.gz · Last modified: 2021/12/20 20:50 by david