Plot Grid cells in actual size

Added by Andreas 9 months ago

Hello,

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 xyz2grd, grdview and 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?

Cheers,
Andi


Replies (2)

RE: Plot Grid cells in actual size - Added by Peter 9 months ago

Hi Andi

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 mapproject
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 mapproject
4. Plot and fill each cell as a polygon using psxy

(perhaps there are simpler ways...)

Best \p

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

(1-2/2)