Malloc error with grdview using a draped TIFF image

Q&A
1

Hi, using GMT master (updated today) I’m getting a malloc() error when using grdview with a draped 3-band GeoTIFF image. I don’t get an error when I use -C to color the DEM using a CPT file, and I also don’t get an error if I split my GeoTiff image into red, green, blue TIFFs first. I imagine this error could be related to the malloc() error in grdimage that was recently fixed?(Getting malloc errors with grdimage of GeoTIFF in oblique Mercator projection(-JOa))

gmt --version
6.4.0_d8429f0_2022.02.20

gmt grdview profiles/P_IR_grid3_newgrid.tif -Gprofiles/P_IR_grid3_colored_hillshade.tif -p -Qi300 -R -J -JZ -O >> profiles/P_IR_profile.ps

gmt(31096,0x10f9a4dc0) malloc: can’t allocate region
:*** mach_vm_map(size=18446744073709514752, flags: 100) failed (error code=3)
gmt(31096,0x10f9a4dc0) malloc: *** set a breakpoint in malloc_error_break to debug
grdview [ERROR]: gdalread: failure to allocate enough memory
grdview [ERROR]: ERROR reading image with gdalread.
grdview (gmtapi_import_image): Could not read from file [profiles/P_IR_grid3_colored_hillshade.tif]
[Session gmt (0)]: Error returned from GMT API: GMT_IMAGE_READ_ERROR (22)
[Session gmt (0)]: Error returned from GMT API: GMT_IMAGE_READ_ERROR (22)

These commands work without a problem:
gmt grdview profiles/P_IR_grid3_newgrid.tif -Cgeo -p -Qi300 -R -J -JZ -O >> profiles/P_IR_profile.ps
gmt grdview profiles/P_IR_grid3_newgrid.tif -Gprofiles/P_IR_grid3_red.tif -Gprofiles/P_IR_grid3_green.tif -Gprofiles/P_IR_grid3_blue.tif -p -Qi300 -R -J -JZ -O >> profiles/P_IR_profile.ps

Cheers,
Kyle

2

Happy to debug this if you can make data available

3

Here’s a ZIP file with the various TIFFs and a test script that will run the two commands that work and also the command that produces the malloc error.

https://www.dropbox.com/s/4xe7ycw9sqx2i9t/malloc_error_data.zip?dl=1

I have also noticed that gmt grdinfo doesn’t report correct Z values for drape.tif, while gdalinfo can read the TIFF just fine:

gmt grdinfo drape.tif

drape.tif: Title: Grid imported via GDAL
drape.tif: Command:
drape.tif: Remark:
drape.tif: Pixel node registration used [Cartesian grid]
drape.tif: Grid file format: gd = Import/export through GDAL
drape.tif: x_min: 0 x_max: 390.2 x_inc: 1.00051282051 name: x n_columns: 390
drape.tif: y_min: -49.825549 y_max: 49.824651 y_inc: 1.00656767677 name: y n_rows: 99
drape.tif: v_min: 1.79769313486e+308 v_max: -1.79769313486e+308 name: z
drape.tif: scale_factor: 1 add_offset: 0
drape.tif: Default CPT:
+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs

gdalinfo drape.tif

Driver: GTiff/GeoTIFF
Files: drape.tif
Size is 390, 99
Coordinate System is:
PROJCRS[“WGS 84 / UTM zone 12N”,
BASEGEOGCRS[“WGS 84”,
DATUM[“World Geodetic System 1984”,
ELLIPSOID[“WGS 84”,6378137,298.257223563,
LENGTHUNIT[“metre”,1]]],
PRIMEM[“Greenwich”,0,
ANGLEUNIT[“degree”,0.0174532925199433]],
ID[“EPSG”,4326]],
CONVERSION[“UTM zone 12N”,
METHOD[“Transverse Mercator”,
ID[“EPSG”,9807]],
PARAMETER[“Latitude of natural origin”,0,
ANGLEUNIT[“degree”,0.0174532925199433],
ID[“EPSG”,8801]],
PARAMETER[“Longitude of natural origin”,-111,
ANGLEUNIT[“degree”,0.0174532925199433],
ID[“EPSG”,8802]],
PARAMETER[“Scale factor at natural origin”,0.9996,
SCALEUNIT[“unity”,1],
ID[“EPSG”,8805]],
PARAMETER[“False easting”,500000,
LENGTHUNIT[“metre”,1],
ID[“EPSG”,8806]],
PARAMETER[“False northing”,0,
LENGTHUNIT[“metre”,1],
ID[“EPSG”,8807]]],
CS[Cartesian,2],
AXIS[“(E)”,east,
ORDER[1],
LENGTHUNIT[“metre”,1]],
AXIS[“(N)”,north,
ORDER[2],
LENGTHUNIT[“metre”,1]],
USAGE[
SCOPE[“Engineering survey, topographic mapping.”],
AREA[“Between 114°W and 108°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Alberta; Northwest Territories (NWT); Nunavut; Saskatchewan. Mexico. United States (USA).”],
BBOX[0,-114,84,-108]],
ID[“EPSG”,32612]]
Data axis to CRS axis mapping: 1,2
Origin = (0.000000000000000,-49.825549000000002)
Pixel Size = (1.000512820512820,1.006567676767677)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 0.0000000, -49.8255490) (115d29’19.48"W, 0d 0’ 1.62"S)
Lower Left ( 0.0000000, 49.8246510) (115d29’19.48"W, 0d 0’ 1.62"N)
Upper Right ( 390.200, -49.826) (115d29’ 6.89"W, 0d 0’ 1.62"S)
Lower Right ( 390.200, 49.825) (115d29’ 6.89"W, 0d 0’ 1.62"N)
Center ( 195.100, -0.000) (115d29’13.19"W, 0d 0’ 0.00"S)
Band 1 Block=390x7 Type=Byte, ColorInterp=Red
Band 2 Block=390x7 Type=Byte, ColorInterp=Green
Band 3 Block=390x7 Type=Byte, ColorInterp=Blue

Thanks again,
Kyle

4

Thanks, will have a look, but just a few comments first:

  • GMT does not really support geographic grids > 360 degrees longitude, since there is no such thing (well, maybe I can image some sort of time series where the 360-390 is another time, but there are better ways to do that). So we may get confused with such a file. Seems like just wasting ~10% space to repeat the data. But, since you are treating it as a Cartesian data set (-JX) the region can be whatever of course.
  • Since drape.tif is an image, grdinfo will not return useful z-ranges. But region and dimensions are right.
  • The -Rxmin/ymin/xmax/ymax+r syntax is used for oblique geographic regions, but you have a Cartesian system (well, see above). So there may be issue related to your 390 degree range here as well.
  • Finally, looks like the image is upside down with y positive downwards. @Joaquim, that gives us a negative nYSize in gmt_gdalread.c and things go downhill from there.
5

Yes my grids are unusual because they are in profile coordinates rather than geographic coordinates. I plot the grids above a swath profile, so the X coordinate is the along-profile distance in km and the Y coordinate is the across-profile distance in km, where Y=0 is the profile line. I extract the R/G/B values from an existing GeoTIFF using grdtrack, and I build a new raster in ‘profile space’ from those. That’s why the raster has a Cartesian coordinate system and the image might appear upside down when inspected.

The profile plots I make end up looking like this:

image
image1676×684 184 KB

6

Also, for some but not all profiles, I am encountering the following memory allocation error from grdview; here is the -Vd output for one example and the grid is attached. I haven’t yet discovered any relationship between grids that work and grids that fail at this step. This grdview call fails with the -Cgeo option so it’s not exactly the same problem as the -G errors above.

P_AE_grid3_newgrid.tif (104.7 KB)

grdview [DEBUG]: History: Process -p
grdview [DEBUG]: History: Process -R
grdview [DEBUG]: History: Process -J
grdview [DEBUG]: History: Process -JZ
grdview [DEBUG]: Look for file 0/-10.001698/664.312/10.001698/-7.26092/1.83124r in /Users/kylebradley/.gmt
grdview [DEBUG]: Look for file 0/-10.001698/664.312/10.001698/-7.26092/1.83124r in /Users/kylebradley/.gmt/cache
grdview [DEBUG]: Look for file 0/-10.001698/664.312/10.001698/-7.26092/1.83124r in /Users/kylebradley/.gmt/server
grdview [DEBUG]: Got regular w/e/s/n for region (0/-10.001698/664.312/10.001698/-7.26092/1.83124r)
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Replace file profiles/P_AE_grid3_newgrid.tif with path profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Replace file profiles/P_AE_grid3_newgrid.tif with profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_red.tif
grdview [DEBUG]: Replace file profiles/P_AE_grid3_red.tif with path profiles/P_AE_grid3_red.tif
grdview [DEBUG]: Replace file profiles/P_AE_grid3_red.tif with profiles/P_AE_grid3_red.tif
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_green.tif
grdview [DEBUG]: Replace file profiles/P_AE_grid3_green.tif with path profiles/P_AE_grid3_green.tif
grdview [DEBUG]: Replace file profiles/P_AE_grid3_green.tif with profiles/P_AE_grid3_green.tif
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_blue.tif
grdview [DEBUG]: Replace file profiles/P_AE_grid3_blue.tif with path profiles/P_AE_grid3_blue.tif
grdview [DEBUG]: Replace file profiles/P_AE_grid3_blue.tif with profiles/P_AE_grid3_blue.tif
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Replace file profiles/P_AE_grid3_newgrid.tif with path profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Replace file profiles/P_AE_grid3_newgrid.tif with path profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Object ID 0 : Registered Grid File profiles/P_AE_grid3_newgrid.tif as an Input resource with geometry Surface [n_objects = 1]
grdview [DEBUG]: gmtapi_begin_io: Input resource access is now enabled [container]
grdview [DEBUG]: gmtapi_import_grid: Passed ID = 0 and mode = 131073
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Calling nc_open on profiles/P_AE_grid3_newgrid.tif, ncid = 0, err = -51
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Calling nc_open on profiles/P_AE_grid3_newgrid.tif, ncid = 0, err = -51
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Calling nc_open on profiles/P_AE_grid3_newgrid.tif, ncid = 0, err = -51
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Calling nc_open on profiles/P_AE_grid3_newgrid.tif, ncid = 0, err = -51
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Calling nc_open on profiles/P_AE_grid3_newgrid.tif, ncid = 0, err = -51
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Calling nc_open on profiles/P_AE_grid3_newgrid.tif, ncid = 0, err = -51
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Calling nc_open on profiles/P_AE_grid3_newgrid.tif, ncid = 0, err = -51
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: Look for file profiles/P_AE_grid3_newgrid.hdr in /Users/kylebradley/.gmt
grdview [DEBUG]: Look for file profiles/P_AE_grid3_newgrid.hdr in /Users/kylebradley/.gmt/cache
grdview [DEBUG]: Look for file profiles/P_AE_grid3_newgrid.hdr in /Users/kylebradley/.gmt/server
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [INFORMATION]: File profiles/P_AE_grid3_newgrid.tif reads with GDAL driver GTiff
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [INFORMATION]: Cartesian input grid
grdview [INFORMATION]: Cartesian input grid
grdview [DEBUG]: GMT_End_IO: Input resource access is now disabled
grdview [DEBUG]: Projected values in meters: 0 664.312 -10.0017 10.0017
grdview [INFORMATION]: Linear projection implies x-axis distance exaggeration relative to the y-axis by a factor of 0.999999
grdview [INFORMATION]: Map scale is 0.0373629 km per cm or 1:3736.29.
grdview [INFORMATION]: Processing shape grid
grdview [DEBUG]: gmtapi_begin_io: Input resource access is now enabled [container]
grdview [DEBUG]: gmtapi_import_grid: Passed ID = 0 and mode = 131074
grdview [INFORMATION]: Reading grid from file profiles/P_AE_grid3_newgrid.tif
grdview [INFORMATION]: gmt_gdalread: Read band(s): 1
grdview [DEBUG]: Found readable file profiles/P_AE_grid3_newgrid.tif
grdview [ERROR]: gdalread: failure to allocate enough memory
grdview [ERROR]: ERROR reading file with gdalread.
grdview (gmtapi_import_grid): Could not open file [profiles/P_AE_grid3_newgrid.tif]
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_READ_ERROR (18)
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_READ_ERROR (18)
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_READ_ERROR (18)
grdview [DEBUG]: gmtlib_garbage_collection: Destroying object: C=0 A=0 ID=0 W=Input F=Grid M=File S=Unused P=7fb19bb04200 N=profiles/P_AE_grid3_newgrid.tif
grdview [DEBUG]: GMTAPI_Garbage_Collection freed 1 memory objects
grdview [DEBUG]: gmtlib_unregister_io: Unregistering object no 0 [n_objects = 0]
gmt [DEBUG]: Entering GMT_Destroy_Session

7

Well, I would assume anything with max > min will give negative y-dimension and hence when computing required memory it will flip around and be seen as a giant number; hence the message about unable to allocate memory.

We will see what @Joaquim says when he gets up since he wrote gmt_gdalread.c to see what we would need to do to support a left-handed coordinate system. As it is it just cannot work since the assumption is right-handed.

8

Alternatively, if you can flip your calculations so that y is up and x is to the right then I think it should work.

9

Hi Paul, I can confirm that when I change my coordinate system so that Y is always positive, I don’t seem to be running into either of these memory allocation errors. It would be more convenient for me to keep the -Y to +Y coordinates as that preserves the profile line at Y=0 and I don’t have to find out the minimum Y value ahead of time and then add it back to the Y coordinates, so if the memory allocation issue can be ironed out I would prefer to keep my code as-is. I can’t always expect X to be positive either, as I sometimes align profiles at a crossing line where X=0, so X goes negative. Thanks again for the insight, the code is at least working now!

10

Update: I still get allocation errors with a grid that has only positive X and Y coordinates, if the maximum X coordinate is too large (>800).

grdview [ERROR]: gdalread: failure to allocate enough memory
grdview [ERROR]: ERROR reading file with gdalread.
grdview (gmtapi_import_grid): Could not open file [profiles/P_AX_grid3_newgrid.tif]
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_READ_ERROR (18)
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_READ_ERROR (18)
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_READ_ERROR (18)

gmt grdinfo profiles/P_AX_grid3_newgrid.tif
profiles/P_AX_grid3_newgrid.tif: Title: Grid imported via GDAL
profiles/P_AX_grid3_newgrid.tif: Command:
profiles/P_AX_grid3_newgrid.tif: Remark:
profiles/P_AX_grid3_newgrid.tif: Pixel node registration used [Cartesian grid]
profiles/P_AX_grid3_newgrid.tif: Grid file format: gd = Import/export through GDAL
profiles/P_AX_grid3_newgrid.tif: x_min: 0 x_max: 863 x_inc: 1 name: x n_columns: 863
profiles/P_AX_grid3_newgrid.tif: y_min: 0 y_max: 99.7279 y_inc: 1.00735252525 name: y n_rows: 99
profiles/P_AX_grid3_newgrid.nc: v_min: -4730.88214153 v_max: 964.602611661 name: z
profiles/P_AX_grid3_newgrid.nc: scale_factor: 1 add_offset: 0
profiles/P_AX_grid3_newgrid.nc: Default CPT:
+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs

11

The origin in GDAL is the upper left corner and Y is positive down but we take care of that. NetCDF can even be written top-down or bottom-up. Don’t know if the same applies to GeoTIFF. So, if that is the cause of the crash we may be missing some check on that.

However, these grids have strange coordinates. Also GDAL reports bad corner coordinates.

gdalinfo grid.tif
...
Corner Coordinates:
Upper Left  (   0.0000000, -49.8255490) (115d29'19.48"W,  0d 0' 1.62"S)
Lower Left  (   0.0000000,  49.8246510) (115d29'19.48"W,  0d 0' 1.62"N)
Upper Right (     390.200,     -49.826) (115d29' 6.89"W,  0d 0' 1.62"S)
Lower Right (     390.200,      49.825) (115d29' 6.89"W,  0d 0' 1.62"N)
Center      (     195.100,      -0.000) (115d29'13.19"W,  0d 0' 0.00"S)

gdalinfo drape.tif
...
Corner Coordinates:
Upper Left  (   0.0000000, -49.8255490) (115d29'19.48"W,  0d 0' 1.62"S)
Lower Left  (   0.0000000,  49.8246510) (115d29'19.48"W,  0d 0' 1.62"N)
Upper Right (     390.200,     -49.826) (115d29' 6.89"W,  0d 0' 1.62"S)
Lower Right (     390.200,      49.825) (115d29' 6.89"W,  0d 0' 1.62"N)
Center      (     195.100,      -0.000) (115d29'13.19"W,  0d 0' 0.00"S)

Those UTM meters are actually km aren’t they?

12

It is not the sign of the coordinates that is the problem, it is the positive direction. The previous grid you shared has y as positive down rows while normal grids have y positive up the rows so that we have a right-handed coordinate system. If gdalinfo says the upper left has y = 0 and lower left has y = 99.72 then it is still going downwards, and gmt_gdalread.c des not check for this case.

13

OK I believe I understand now. The issue arose from creating the rasters using a gdal_grid command of the form:

gdal_grid -txe $dem_minx $dem_maxx -tye $dem_miny $dem_maxy -outsize $numx $numy …

when it needed to be (less intuitively I think):

gdal_grid -txe $dem_minx $dem_maxx -tye $dem_maxy $dem_miny -outsize $numx $numy …

Searches reveal that this behavior of gdal_grid is pretty well known. I apparently made grdview work with the incorrect raster orientation by flipping the sign of the Y data in the input file to gdal_grid, producing some kind of double-flipped raster. This didn’t come back to bite me until I upgraded gdal/GMT recently. Not sure whether you want to be able to read these kinds of rasters in, or just output a warning when detecting the flipped coordinates?

Thanks again for the help and the chance to learn - Kyle

14

Thanks Kyle, we will probably start with putting in an error message and exit rather than crashing. If we are able to handle the reverse y then we will add that when we find the time and wisdom. Cheers, Paul

15

Thankfully, @Joaquim was able to update the code to handle the upside down y-cases. With that PR now merged into the repo you should be able to run your scripts - your test.sh works for me, for instance.