Can grdimage plot grids where the latitude bands are unequally spaced?

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, gridimage will 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?

The image looks “about” right, but that is because the latitude nodes are approximately equally spaced.

Here are the “real” latitudes in the netcdf file

[ 75.50288341, 56.72335728, 37.83368077, 18.92035001, 0. , -18.92035001, -37.83368077, -56.72335728, -75.50288341]

And here is what grdinfo says

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.

Hi everyone, I encountered a problem here. I used the shell command in Figure 1 to get the result in Figure 2. It looks very bad.
Figure 1:


Figure 2:

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.

@SUN-Ravenclaw please do not reuse old topics. Just start a new one. Ah, and also post your commands in plain text, not screen captures.

Hi Paul,

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.

Thanks!

-jose

I agree that grd2xyz is not going to do what we want here. Maybe you can run ncdump -h on that grid and share the output here?

Hey Paul, here’s some info on the grid.

Firstly, using grdinfo -L0 i get:

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

netcdf NPac_Tohokuh {
dimensions:
	lon = 300 ;
	lat = 287 ;
	grid_lon = 1200 ;
	grid_lat = 1148 ;
	time = 1441 ;
variables:
	double lon(lon) ;
		lon:long_name = "longitude" ;
		lon:units = "degrees_east" ;
		lon:point_spacing = "even" ;
		lon:resolution = "16.00 arcmin" ;
	double lat(lat) ;
		lat:long_name = "latitude" ;
		lat:units = "degrees_north" ;
		lat:point_spacing = "uneven" ;
		lat:resolution_min = "7.95 arcmin" ;
		lat:resolution_max = "16.00 arcmin" ;
	int crs ;
		crs:grid_mapping_name = "latitude_longitude" ;
		crs:longitude_of_prime_meridian = 0. ;
		crs:semi_major_axis = 6378137. ;
		crs:inverse_flattening = 298.257223563 ;
		crs:crs_wkt = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]" ;
		crs:spatial_ref = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]" ;
	double grid_lon(grid_lon) ;
		grid_lon:long_name = "Grid Longitude" ;
		grid_lon:units = "degrees_east" ;
		grid_lon:point_spacing = "even" ;
		grid_lon:resolution = "4.00 arcmin" ;
	double grid_lat(grid_lat) ;
		grid_lat:long_name = "Grid Latitude" ;
		grid_lat:units = "degrees_north" ;
		grid_lat:point_spacing = "uneven" ;
		grid_lat:resolution_min = "1.99 arcmin" ;
		grid_lat:resolution_max = "4.00 arcmin" ;
	float deformation(grid_lat, grid_lon) ;
		deformation:long_name = "Grid Deformation" ;
		deformation:units = "meters" ;
		deformation:grid_mapping = "crs" ;
		deformation:missing_value = -1.e+34f ;
		deformation:_FillValue = -1.e+34f ;
	float max_height(lat, lon) ;
		max_height:long_name = "Maximum Wave Amplitude" ;
		max_height:units = "cm" ;
		max_height:grid_mapping = "crs" ;
		max_height:_FillValue = -1.e+34f ;
		max_height:missing_value = -1.e+34f ;
	float travel_time(lat, lon) ;
		travel_time:long_name = "Travel Time" ;
		travel_time:units = "hours" ;
		travel_time:grid_mapping = "crs" ;
		travel_time:_FillValue = -1.e+34f ;
		travel_time:missing_value = -1.e+34f ;

and here’s a sample of the ‘ncdump -v lat’

data:
 lat = -0.2834, -0.0167, 0.25, 0.5166, 0.7833, 1.0499, 1.3165, 1.5831, 
    1.8496, 2.1162, 2.3826, 2.649, 2.9154, 3.1817, 3.4479, 3.714, 3.9801, 
    4.2461, 4.512, 4.7778, 5.0434, 5.309, 5.5745, 5.8398, 6.105, 6.3701, 
    6.6351, 6.8999, 7.1646, 7.4291, 7.6934, 7.9576, 8.2216, 8.4854, 8.7491, 
    9.0126, 9.2759, 9.5389, 9.8018, 10.0645, 10.3269, 10.5892, 10.8512, 
    11.113, 11.3745, 11.6358, 11.8969, 12.1577, 12.4183, 12.6785, 12.9386, 
    13.1983, 13.4578, 13.717, 13.9759, 14.2346, 14.4929, 14.7509, 15.0087, 
    15.2661, 15.5232, 15.78, 16.0364, 16.2925, 16.5483, 16.8038, 17.0589, 
    17.3136, 17.568, 17.8221, 18.0758, 18.3291, 18.5821, 18.8346, 19.0868, 
    19.3386, 19.5901, 19.8411, 20.0917, 20.342, 20.5918, 20.8412, 21.0902, 
    21.3388, 21.587, 21.8348, 22.0821, 22.329, 22.5754, 22.8214, 23.067, 
    23.3121, 23.5568, 23.801, 24.0448, 24.2881, 24.5309, 24.7733, 25.0151, 
    25.2566, 25.4975, 25.738, 25.9779, 26.2174, 26.4564, 26.6949, 26.9329, 
    27.1704, 27.4073, 27.6438, 27.8798, 28.1153, 28.3502, 28.5846, 28.8185, 
    29.0519, 29.2848, 29.5171, 29.7489, 29.9801, 30.2108, 30.441, 30.6707, 
    30.8997, 31.1283, 31.3563, 31.5837, 31.8106, 32.0369, 32.2627, 32.4879, 
    32.7126, 32.9367, 33.1602, 33.3832, 33.6055, 33.8273, 34.0486, 34.2692, 
    34.4893, 34.7088, 34.9278, 35.1461, 35.3639, 35.581, 35.7976, 36.0136, 
    36.229, 36.4438, 36.658, 36.8717, 37.0847, 37.2971, 37.509, 37.7202, 
    37.9308, 38.1409, 38.3503, 38.5591, 38.7674, 38.975, 39.182, 39.3884, 
    39.5942, 39.7994, 40.0039, 40.2079, 40.4112, 40.614, 40.8161, 41.0176, 
    41.2185, 41.4188, 41.6184, 41.8175, 42.0159, 42.2137, 42.4109, 42.6075, 
    42.8035, 42.9988, 43.1935, 43.3876, 43.5811, 43.774, 43.9662, 44.1578, 
    44.3488, 44.5392, 44.729, 44.9181, 45.1066, 45.2945, 45.4818, 45.6685, 
    45.8545, 46.04, 46.2248, 46.4089, 46.5925, 46.7754, 46.9577, 47.1395, 
    47.3205, 47.501, 47.6808, 47.8601, 48.0387, 48.2167, 48.394, 48.5708, 
    48.7469, 48.9225, 49.0974, 49.2717, 49.4454, 49.6184, 49.7909, 49.9627, 
    50.134, 50.3046, 50.4746, 50.644, 50.8128, 50.981, 51.1486, 51.3156, 
    51.482, 51.6477, 51.8129, 51.9775, 52.1414, 52.3048, 52.4675, 52.6297, 
    52.7912, 52.9522, 53.1126, 53.2723, 53.4315, 53.5901, 53.7481, 53.9054, 
    54.0623, 54.2185, 54.3741, 54.5291, 54.6836, 54.8374, 54.9907, 55.1434, 
    55.2955, 55.4471, 55.598, 55.7484, 55.8982, 56.0474, 56.1961, 56.3442, 
    56.4916, 56.6386, 56.7849, 56.9307, 57.0759, 57.2206, 57.3647, 57.5082, 
    57.6512, 57.7936, 57.9354, 58.0767, 58.2175, 58.3577, 58.4973, 58.6363, 
    58.7748, 58.9128, 59.0502, 59.1871, 59.3234, 59.4592, 59.5944, 59.7291, 
    59.8633, 59.9969, 60.13, 60.2625 ;
}

the lat spacing gets continuously smaller to the north.

If I just plot it straight using grdimage, i get:

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

anyhoo, when i plot it out it comes out as below:

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!

Also, you might like this one too:

Thanks!

-jose

Hi José

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.

1 Like

Cheers @Joaquim… that’s what I thought.

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!

-jose

meshgrid does this

>> [X,Y] = meshgrid([1,2,3],[11,12,13,14])

X =

     1     2     3
     1     2     3
     1     2     3
     1     2     3


Y =

    11    11    11
    12    12    12
    13    13    13
    14    14    14

that is, replicates X length(Y) times and vice-versa for Y

It shouldn’t be that difficult to have grd2xyz output the correct x,y,z triplets for those grids but it was never implemented.

AGU days will come back … but money for them will take longer.

can it be a feature request?

Sure. One can always ask.

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

So yeah,

Thanks!
Tusen takk!
Muito Obrigado!

-jose