Grdedit -A seems not to work

I have a GMT grid that I created with grdlandmask to match a grid that I created with gdal_translate. When I try to use grdmath to multiply the mask and the data, grdmath complains:

grdmath [ERROR]: (x_max-x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001.
grdmath (gmtlib_read_grd_info): Use grdedit -A on your grid file to make region and increments compatible [landmask.grd]

I then tried to run “gmt grdedit -A landmask.grd -Glandmask-A.grd” on my landmask.grd, but the output still gives the same error about the increment mismatch.

gmt grdinfo landmask-A.grd
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 [landmask-A.grd]
landmask-A.grd: Title: z
landmask-A.grd: Command: grdlandmask -Glandmask.grd -R/Users/fielding/Google_Drive_JPL/Mexico-Drive/Oaxaca-EQ-2020/S1AB/A107/ARIA_proc/results/S1_A107_20200625_20200619_disp_cm.grd -Df
landmask-A.grd: Remark: Derived from the full resolution shorelines
landmask-A.grd: Pixel node registration used [Geographic grid]
landmask-A.grd: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
landmask-A.grd: x_min: -96.5451065365 x_max: -93.9876553945 x_inc: 0.00083331741351 (3 sec) name: longitude n_columns: 3069
landmask-A.grd: y_min: 15.046045981 y_max: 16.7804535652 y_inc: 0.000833449103383 name: latitude n_rows: 2081
landmask-A.grd: z_min: 0 z_max: 1 name: z
landmask-A.grd: scale_factor: 1 add_offset: 0
landmask-A.grd: format: netCDF-4 chunk_size: 134,131 shuffle: on deflation_level: 3

Am I missing something about how to use “grdedit -A”?
I am using GMT 6.0.0 installed with MacPorts.

Hm, if you can make that large grid available then I can have a look. Does the ARIA grid really have slightly different dx/dy or is this all sloppy rounding on that end?

The original ARIA InSAR product is available online (https://aria-share.jpl.nasa.gov/20200623_Oaxaca_Mexico_EQ/Displacements/S1-GUNW_COSEISMIC-A-R-107-tops-20200625_20200619-003043-16794N_15045N-RR-35f2-v2_0_3.nc), and it should be exactly 3-arcsecond spacing. It is stored in a complicated multilayer NetCDF HDF5 format, so I use some Python tools to extract the displacement grid into a GeoTIFF, then I used gdal_translate to convert the GeoTIFF to a GMT-format grid. I am not sure why the dx and dy are slightly different but it seems that sloppy rounding must have happened somewhere. I was not able to read the GeoTIFF directly with GMT6.0.

Long, but simple.

gmt grdconvert S1-GUNW_COSEISMIC-A-R-107-tops-20200625_20200619-003043-16794N_15045N-RR-35f2-v2_0_3.nc=gd?NETCDF:"S1-GUNW_COSEISMIC-A-R-107-tops-20200625_20200619-003043-16794N_15045N-RR-35f2-v2_0_3.nc:/science/grids/data/amplitude" -Gamp.grd

increments are exactly equal

grdinfo amp.grd
amp.grd: Title: Produced by grdconvert
amp.grd: Command: grdconvert S1-GUNW_COSEISMIC-A-R-107-tops-20200625_20200619-003043-16794N_15045N-RR-35f2-v2_0_3.nc=gd?NETCDF:S1-GUNW_COSEISMIC-A-R-107-tops-20200625_20200619-003043-16794N_15045N-RR-35f2-v2_0_3.nc:/science/grids/data/amplitude -Gamp.grd
amp.grd: Remark:
amp.grd: Gridline node registration used [Geographic grid]
amp.grd: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
amp.grd: x_min: -96.5533333333 x_max: -93.91 x_inc: 0.000833333333333 (3 sec) name: longitude n_columns: 3173
amp.grd: y_min: 15.045 y_max: 16.795 y_inc: 0.000833333333333 (3 sec) name: latitude n_rows: 2101
amp.grd: z_min: 0 z_max: 754838.0625 name: z
amp.grd: scale_factor: 1 add_offset: 0
amp.grd: format: netCDF-4 chunk_size: 133,132 shuffle: on deflation_level: 3
GEOGCS["unknown",
    DATUM["Unknown based on WGS84 ellipsoid",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Longitude",EAST],
    AXIS["Latitude",NORTH]]

Thanks @Joaquim, so the gdal_translate route messes up node registration and resets increments etc. I think that is the root of the problems, @EJFielding. The number of nodes are very different from the two files - sure it is the same file?

Yes, the data came from the same original file, but obviously there was some resampling or cropping at some step. I did not try to use grdconvert directly on the original file. That looks like the best option for the simple layers. There are some layers that are stored at different resolutions that have to be resampled to be used, but the amplitude and unwrapped phase of the interferogram are simple layers. Thanks, @Joaquim !

I heard that GDAL really only knows about pixel registration for grids, so it must be resampling the data to force the pixel registration. The original ARIA interferogram product is geocoded to the SRTM DEM, which is grid registration, so I am glad to see that reading the original file directly sees that and gets the correct information.

The Python tools for reading the ARIA standard products are at:

I will open an issue there to discuss the problem of the dx and dy changing.

Yes, it’s true that for GDAL everything is pixel reg, but I’m not sure that it will resample. It’s just much simpler to to adjust the limits by half a cell width. That is exactly what we do when we ask GDAL to read a file/layer and return the data to GMT. Like the grdconvert example above.

It is probably the ariaExtract.py program from ARIA-tools that is resampling. Maybe you saw there are some layers in the ARIA interferogram product that are 3-D layers at coarse resolution. Those have to be resampled to match the unwrapped phase and amplitude images. I guess that means the code also resamples the full-resolution layers. Bottom line is that I think the ultimate bug is in the ARIA-tools ariaExtract writing out files with inconsistent limits and spacing.

It is still strange that GMT grdedit -A can’t adjust the spacing or region so that other GMT programs can use the grid files.

Yes, I saw it but I’m almost sure you can do the resampling with GMT alone by giving -Ramp.grd (to use my example).

But regarding grdedit -A, I suspect you are right. A couple of months ago there was another forum thread where it apparently also failed to do the promised.

Offer to look at your file that fails in grdedit -A still stands!

I have uploaded the GMT grd file that I used as the region for the grdlandmask command that made the landmask.grd file that grdedit -A fails to correct. It is at:
https://aria-share.jpl.nasa.gov/20200623_Oaxaca_Mexico_EQ/Displacements/S1_A107_20200625_20200619_disp_cm.grd

Note that I will likely delete this file later as it is obviously not a correct GMT grid.

I am afraid that GMT can’t do the resampling of the 3D grid layers directly from the ARIA product. They require a 3D interpolation using an additional DEM, not a basic 2D interpolation. I could use the “hammer” approach of resampling the files to a correct GMT grid with -Ramp.grd after extracting them from the ARIA product with ariaExtract. The extra layers are quite smooth so extra resampling should not hurt much.

OK, but keep an eye on #2966

I got your grid Eric so you can remote it any time.

1 Like

I can reproduce the eps messages. Interestingly, they don’t come if I run grdinfo your grid (which has hardly any netcdf metadata in it) but once GMT writes out that landmask it complains if I run grdinfo on that one, despite having exactly same -R -I -r as your start file…

Thanks. Good to see you folks are working on the 3D grids.

The ARIA interferogram format is a prototype for the NISAR mission interferogram products, so it would be great to eventually add direct reading of this kind of file in GMT. NISAR launch is not until 2022, so there is still plenty of time ;-).

I had to use

gmt grdconvert 'S1-GUNW_COSEISMIC-A-R-107-tops-20200625_20200619-003043-16794N_15045N-RR-35f2-v2_0_3.nc=gd?NETCDF:"S1-GUNW_COSEISMIC-A-R-107-tops-20200625_20200619-003043-16794N_15045N-RR-35f2-v2_0_3.nc":/science/grids/data/amplitude' -Gamp.grd

on my Mac.

@pwessel I am glad you can reproduce the eps messages. It seems there is something going subtly wrong somewhere.

Note that the ARIA-tools software has been revised to provide output with original 3-arcsecond spacing and consistent region values.