Can grdimage plot gridded data where the latitudinal bands are unequally spaced? I am thinking particularly of the case where the grid was defined to facilitate integrations using Gauss-Legendre quadrature.
Now, if I have a netcdf grid with unevenly spaced latitude bands, gridimagewill plot it. However, if you use grdinfo on the file, it provides a single latitudinal spacing. So either grdinfo is providing misleading information, or grdimage is ignoring the values of the latitudes and assuming that they are equally spaced.
@mrak this is a much requested feature in GMT but I think itās not supported (@pwessel can correct me here if Iām wrong). It would be helpful if you had an example grid and plot for testing. Out of curiosity, what output do you get from grdimage?
y_min: -75.5028834059 y_max: 75.5028834059 y_inc: 18.8757208515 name: latitude n_rows: 9
The increment appears to be just (max-min)/(n-1). I personally donāt care if gmt supports this or not, but if they donāt, a warning or error should be issued (pygmt catches this sometimes, but not always).
No, GMT (any module) cannot deal with non-regular grids.
I remember that in the past we had issues with grids that had coordinate vectors not exactly equi-spaced. Itās probably because of that an average increment is calculated ā¦ without looking into the std.
Best work-around is to dump your grid to a table using grd2xyz and regrid onto an equidistant lattice, then plot that grid. Depending on the size of the table (number of points) I would recommend greenspline (especially for global grids) provided < 10,000 points.
I have two questions here:
(1) How to adjust the background color to gray or other color to make it looks well.
(2) Actually, the data range I use is: latitude (39 Ā° -40 Ā°), longitude (111 Ā° -113 Ā°), but it can be seen from Figure 2 that these points are completely offset from the correct position. How do I modify the code? Thank you.
Since you change your region in the gmt plot call then nothing will line up properly anymore. Remove that -R - you want to keep it the same as the background region.
If you want a grayshade for the topo then change the cpt you are using for the FRM to gray or any other - that is your choice to make.
I have this issue with some data sets that i frequently work with (the NOAA global tsunami propagation database), but dumping the data using grd2xyz just gives me an xyz table with the xās and yās locations at the average increments, not at the actual values in the original data file. So then when you go to re-grid it, you just get the same thing back again.
weāve had to set up a MATLAB based work around to interpolate it on to a regular lattice, then bring it back to GMT for plotting.
Iād love to have a way to do it all internally in GMT, but havenāt been successful in getting it to work.
I have the examples worked out on another computer so i can show you what I mean at some point.
grdinfo [ERROR]: (x_max-x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001.
grdinfo (gmtlib_read_grd_info): Use grdedit -A on your grid file to make region and increments compatible [../data/NPac_Tohokuh.nc?max_height]
../data/NPac_Tohokuh.nc: Title: MOST tsunami propagation model output ki24b1
../data/NPac_Tohokuh.nc: Command: 2020-04-04 21:18:23 NZDT: propcomb -v -t -m -n -h -3 -x 121 200 0 60 -p NPac_Tohoku -l ../propdb/compressed -d 4.66*ki24b+12.23*ki25b+26.31*ki26a+21.27*ki26b+22.75*ki27a+4.98*ki27b
../data/NPac_Tohokuh.nc: Remark:
../data/NPac_Tohokuh.nc: Gridline node registration used [Geographic grid]
../data/NPac_Tohokuh.nc: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
../data/NPac_Tohokuh.nc: x_min: 120.55 x_max: 200.2833 x_inc: 0.266666555184 (16 min) name: longitude n_columns: 300
../data/NPac_Tohokuh.nc: y_min: -0.2834 y_max: 60.2625 y_inc: 0.211698951049 name: latitude n_rows: 287
../data/NPac_Tohokuh.nc: z_min: 0 z_max: 793.649353027 name: Maximum Wave Amplitude [cm]
../data/NPac_Tohokuh.nc: scale_factor: 1 add_offset: 0
../data/NPac_Tohokuh.nc: format: classic
Then using ncdump -h: (note there are two sets of lat/long dimensions ālatā and ālonā are coarser spaced for the entire model output domain and āgrid-latā and āgrid-lonā which are at higher resolution and only cover the earthquake deformation area, these are for plotting the initial deformation field. We are concerned with ālatā and ālonā but both have the same issue.)
if i use grdsample i get the same thing. You can see how the data gets progressively shifted from reality the further from the pole you are.
Writing the data out to an xyz using grd2xyz puts the lats at the average spacing of 0.2117 deg, so, if i regrid it, it comes out the same and plots up as above.
When we use our MATLAB script, we write out the array of lats and longs and then that goes in to an irregularly spaced xyz which then gets regridded to a spacing which we can control. In this case i used 0.1 deg (artificially interpolated, i knowā¦)
So yeah, the problem as I see it is that i canāt write the actual x,y, positions from an irregularly spaced netCDF matrix to an xyz for regridding using xyz2grd since it always uses the average spacing (max_value - min_value)/n_spaces.
But it is entirely possible (and quite likely!) that I am doing something stupidā¦
Anyway, for your enjoyment, here is an animation I made with GMT from this data set for the 150th anniversary of the 1868 Arica tsunami
I am forever grateful for the hack/tip you gave me a few years ago to make the counter for the hours!
I know what you mean but we donāt have a way to handle that in GMT. We would need to have the equivalent of the meshgrid function, which we donāt.
We did a similar thing in our GMTMEX paper. You can find the Julia translation (very similar to the matlab code) of that example here
BTW, I think we met in a AGU some years ago and talked about tsunamis.
I am not sure exactly what meshgrid does, but that is what we use (as well as interp2) to produce the interpolated grid in MATLAB. then we write it back in to an interpolated netCDF and run that through the plotting script.
seems like it should be easy to read the specific array (lon,lat) values in an irregularly spaced netCDF grid and write them directly to an (x,y,z), but by āeasyā i mean in my imagination and not in my practical ability!
AGUā¦ man itās been a long time! I hope to get back one of these years when weāre allowed to travel and see people again!
@pwessel@Joaquim just got around to downloading and installing 6.1 and tried out the new feature in grd2xyz and it works great! eliminates the need for me to go to MATLAB for that one workaround step for these particular data sets.
