I received a message from a colleague (Ola — GitHub link here ) who is working with PyGMT and adapting the “Antarctic Peninsula 3D Maps” repository. The goal is to drape an RGB image over a bathymetric grid of the Weddell Sea using the projection=‘L-60/-65/-72/-60/15c’.
However, the final result shows the image partially clipped, and it’s unclear why.
I ran several tests on my side and initially suspected the issue was related to the workflow for cropping the satellite image and aligning it with the grid. But then I realized something interesting: just by changing the perspective, the clipping changes position. This makes me believe the issue may be related to the rendering process (see image above for an example).
Does anyone have an idea of what might be causing this?
If you’re curious or want to try reproducing the problem, my notebook is available here: [link].
I see this bug was not fixed since the last year. Do you have any recommendation for a workaround? Otherwise, I will just try out different perspectives and projections to see if I find something that looks good at high latitudes and renders correctly.
The workaround in Strange data cuting with grdview is not really an option. Changing the surface type as a surface plot s instead of image i does not make sense, when the whole point of the figure is to drape an image.
Ole,
In the tests I ran, I noticed that this issue occurs when using an image as the drapegrid, but not when using a regular (netcdf) grid.
So, one possible workaround would be to separate the RGB image into three individual grids (R, G, and B) and then recompose them on top of the original grid with a specific cpt. I remember seeing a discussion about a workaround like this in the past. I’ll run some tests and see how it goes. I’ll keep you posted.
Just to be clear…I’m looking deep into grdview.c and I think the most likely origin of the “cut” or cropping observed in the draped image lies in how the gmt_geoz_to_xy() function transforms geographic coordinates (longitude, latitude, elevation) into projected coordinates (x, y) during the generation of the 3D model.
In the code, this transformation occurs at several points, more importantly, when grdview defines the sweep order of the 3D tiles (grid patches) based on the view perspective, it calls gmt_map_setup(…) and grdview_set_loop_order(…).
At that stage, it determines the sweep direction using the visual octant, computing the azimuth as az = d_atan2d(y1 - y0, x1 - x0) and that sets the direction in which tiles are rendered.
When a drape image is loaded, it is converted into three separate R, G, B grids. These grids are forced to match the topographic grid by adopting its extent and resolution, and in this process, grdview resamples the image to fit the topography. However, it assumes that the map projection preserves alignment — which might not hold true in edge cases, especially with -JL (Lambert) projections.
A deeper debug of the projection and rendering pipeline would be necessary to fully resolve this, which adds complexity.
However, the solution I propose in your case is to convert the drape image into a NetCDF grid with three bands (R, G, B), aligned with the topographic grid, and then reconstruct a colormap (CPT) that reflects the image content. This approach should ensure consistent rendering under any projection.
And why do you think an image cannot be draped on top of a -Qs surface? I didn’t test but I don’t remember such a limitation exists (though I may be wrong ofc).
If I change the surface type from image to surface, the rendering jumps from at maximum a few minutes to not being done after 20 minutes. Maybe it is feasible for lower resolution drape grids, but apparently not for my application.
Some of the other projections also trigger the bug and cut off the drapegrid. Mercator works, but distorts my region from 60°S to 80°S more than what I like to accept.
Following up on an issue raised by @opinner, we’ve been investigating a problem when trying to drape a satellite image over antartic bathymetry (high latitude) using grdview. The goal was to produce a 3D perspective of the Weddell Sea using a satellite raster and a projected base grid (as in the Antarctic_Peninsula_3D_Maps example).
Unfortunately, after multiple attempts and a careful inspection of the code, we couldn’t find a definitive fix for the image being partially cropped during the grdview rendering. This seems to occur particularly when:
Using non-cylindrical projections (e.g., polar stereographic),
Combining large raster images with high latitudes,
Or applying drapegrid with reprojected inputs that rely on split the image and adjust the tiles.
At this point, it looks like a complex interplay between projection distortion and image boundaries.
Given that, I’d like to suggest adding a warning or note in the PyGMT documentation (under grdview) indicating that:
“When using drapegrid with high-latitude or strongly projected grids, users may encounter clipping or unexpected behavior, particularly with satellite imagery.”
If the community agrees, I’d be happy to submit a pull request (or open an Issue focused just to correct the docs) with a concise warning in the appropriate section of the docs.
the above is not the same as the below, so the proposed remark appears to me as misleading.
now it’s time for a question: why not subsampling the raster image before the draping? does this help/fixes the rendering, e.g. in combination with -Qs, as a fix for a way too long rendering time? It sounds like the resolution of the draped image may be unnecessarily high.
