I am trying to create a contour of a misfit function on a global map where the minimum misfit is close to the pole. As you can see in the plot, the contour is plotted correctly, but it draws a line connecting the two points on the left and right of the figure. Does anyone know of a way to fix this?
This image was created with pygmt, but I suspect that this is also a problem with gmt itself.
hi, please add your code
These are hard bugs to find. We know of a few circumstances where it happens. Several of these had been fixed in the past but some cases insist in resisting.
Can you please dump that contour and attach it to this post?
It turns out that this in fact works correctly with gmt, and that the problem is with pygmt.
In pygmt, the grdcontour method works fine if you pass the name of a netcdf file. However, when passing an xarray dataarray, this is when you get the odd behaviour shown above.
I’ll try to come up with a reproducible example and then open an issue with them.
In the end, GMT is still the one to blame. My guess is that a pixel vs grid registration difference is what is what triggers the bug via pygmt. Again, can you dump that contour with both pyGMT and plain GMT and see if how they differ?
Here is the code used in ipython to (1) create an xarray datarray using sphinterpolate, (2) export this as a netcdf file, and (3) export the contour using a gmt shell command:
da = pygmt.sphinterpolate(np.array([lon, lat, misfit]).T, spacing=0.2, region='g')
da.to_netcdf('test.nc')
!gmt grdcontour test.nc -C0.40494052, -D > contour.dat
This shows that the dataarray is by default using pixel registration and geographic coordinates
In [92]: da.gmt.registration, da.gmt.gtype
Out[92]: (0, 1)
And here is the contour file
contour.dat (64.8 KB)
And for completeness, here are the two ways of creating the pygmt image, the first doesn’t work, the second does:
fig = pygmt.Figure()
fig.grdimage(da, projection='W0/10c')
fig.grdcontour(da, levels='+0.40494052')
fig.show()
fig2 = pygmt.Figure()
fig2.grdimage('test.nc', projection='W0/10c')
fig2.grdcontour('test.nc', levels='+0.40494052')
fig2.show()
And here is the pygmt issue: problems with grdcontour when using xarray DataArrays · Issue #3479 · GenericMappingTools/pygmt · GitHub
Nope. 0 means grid registration.
You are right. That was a typo. In any case, it has no influence on the underlying problem. You get the same result setting it to 0 or 1.
But you don’t set it to 0 or 1. What determines its value is if the grid is pixel or gridline registered, and on the PyGMT issue I see these important warning messages that clearly indicate that there is an error in creating that grid
grdimage [WARNING]: Guessing of registration in conflict between x and y, using gridline
grdcontour [WARNING]: Guessing of registration in conflict between x and y, using gridline
Plotting that contour alone works just fine.
plot contour.dat -Rd -JW0/10c -Baf -BWSen -png lixo