Plot Grid cells in actual size
I'm pretty new to GMT, but so far I really like the possibilities. I have written a small tool (in Python), which generates ground tracks which are stored in a csv file. This file is then processed by GMT, in order to see how often the same spot on earth (in a 30km x 30km grid) is visited. However, I can not find the right option to fill a 30kmx30km cell completly with a color (the color is dependend on how many times the cell is visited).
Here is what I do:
# Draw basemap with coastlines gmt pscoast -Rd -JQ6i -B60f30g30 -N1 -Sblue -Ggray -P -K > points.ps # Resample input data to get continous track (since the actual sample rate of the points is 30s) gmt sample1d -I30k -Afr -gd50k -h1 -i7,6,0 "%%F" > "%%F_resampled.csv" # Select all the points over water surfaces gmt select "%%F_resampled.csv" -A -Nk/s/k/s/k > "%%F_ocean.csv" # Create grid and count the number of occurances gmt blockmean -C -Sn -I30k "%%F_ocean.csv" -h1 -i0,1 -Rd > "%%F_blockmean_ocean.csv" gmt psxy "%%F_blockmean_ocean.csv" -Rd -JQ6i -h0 -i0,1,2s0.025 -Scc -Ggreen -O -K >> points.ps
These files (*_blockmean_ocean.csv) are then processed further and in the end used by:
gmt pscoast -R-180/180/70/90 -JS0/90/11c -B20g:."Polar (70 to 90 deg)": -N1 -Sblue -Ggray -P -K > "%%F_pol.ps" gmt psscale -Dx5i/0i+w3i/0.25i -Csnr.cpt -B1:."Occurance": -O -K >> "%%F_pol.ps" # Plot the points ... Here is the command I'm missing gmt psxy "%%F" -R-180/180/70/90 -JS0/90/11c -h1 -i1,2,3 -Sc0.025c -Csnr.cpt -O >> "%%F_pol.ps"
The result has points plotted at the center of the grid, in the correct color (i.e. what I expected). But as mentioned above, the grid cell is not filled. I have played around with
surface, but the plots do not look what I expext (especially around the poles). Is this the correct way to do? If so, what would be the right command to get the complete cells filled?
RE: Plot Grid cells in actual size - Added by Peter 9 months ago
Afaik only the
-SE option can handle geographical distances currently (and it seems you want a rectangle rather than an ellipse). Would be nice though is the size parameter also accepted the units described under Units in the man page so perhaps something for the developers to think about but might require that you post a feature request for it)
A work around for now would be to plot and fill a polygon instead. This would require some pre-processing but the basic steps would be to:
1. transform the mean position of the cells from the output of
blockmean that you want to plot to projected units using
2. for each cell find the position of the corners by subtracting/adding 15 km.
3. transform the position of the cell corners back to geographical coordinates using
4. Plot and fill each cell as a polygon using
(perhaps there are simpler ways...)
RE: Plot Grid cells in actual size - Added by Paul 9 months ago
You cannot make a 30x30 km grid since it would not be equidistant or have the same number of columns per row. The simplest option would be to use a geographic grid (fixed dx, dy in degree) and repeat your exercise, then adjust for the shrinking area with latitude. Also note that your resampled line may have several points inside any given cell and thus be counted more than once. I assume you only want to increment the count for each line going through a cell, which is slightly different. If you really want equal area cells then you will need to convert your lines to equal area projection coordinates and do the counting on a projected grid (which would be equal area) and then plot that grid using a linear projection (since the grid is already projected). I recommend using cylindrical equal area (-JY).