What is the correct way to plot unstructured triangular grids?

[guidocioni]

I'm trying to slowly replace basemap with cartopy but finding myself struggling in even the most simple applications. For example I often have to deal with data on an unstructured grid, that is 1-dimensional data where the information on every cell position is described by two 1-D non-ordered vectors of longitude and latitude values.

With basemap I would usually take the lat/lon vectors from the file, convert them from radians to degrees and then plot them with tricontourf or contourf with Tri=True. This works in most of the cases (although not always). So with Cartopy I was thinking of doing something similar but it just does not work.

fig = plt.figure(figsize=(16,9))
ax = plt.axes(projection=ccrs.NearsidePerspective(
                            central_latitude=50.,
                            central_longitude=-25.,
                            satellite_height=4e6))

ax.add_feature(cfeature.BORDERS), linestyle='-', alpha=.5,
           edgecolor='white', linewidth=2.)

ax.tricontourf(lon, lat, dset.t[0,0,:], transform=ccrs.PlateCarree())

A blank plots filled with NaN is produced

If I omit the transform keyword then the data gets correctly plotted but without having the correct representation on the map, as tricontourf would usually do also without cartopy.

Am I missing something really basic here? I have not seen ANY documentation on how to use cartopy with unstructured grids besides the side project psy-maps, but I want something that is more flexible and works with matplotlib objects.

[dopplershift]

So I'm guessing it doesn't work because Cartopy doesn't wrap tricontourf and there is usually some massaging that needs to be done. You should still be able to take the basemap approach:

fig = plt.figure(figsize=(16,9))
proj = ccrs.NearsidePerspective(
                            central_latitude=50.,
                            central_longitude=-25.,
                            satellite_height=4e6)

x, y, _ = proj.transform_points(ccrs.PlateCarree(), lon, lat).T
ax = plt.axes(projection=proj)

ax.add_feature(cfeature.BORDERS), linestyle='-', alpha=.5,
           edgecolor='white', linewidth=2.)

ax.tricontourf(x, y, dset.t[0,0,:])

I expect this won't work perfectly, but I'd be interested to see how/if it fails.

[guidocioni]

It fails because of the same bug that was in the basemap implementation (see matplotlib/matplotlib#4424) The problem is that points outside of the globe in the ortographic projections attain really large values (in basemap it was 1e30, in cartopy apparently they are inf). This can be easily avoided by masking the array beforehand.

The following code runs in about 25 seconds.

fig = plt.figure(figsize=(16,9))
proj = projection=ccrs.NearsidePerspective(
                            central_latitude=50.,
                            central_longitude=-25.,
                            satellite_height=4e6)
x, y, _ = proj.transform_points(ccrs.PlateCarree(), lon, lat).T

mask = np.invert(np.logical_or(np.isinf(x), np.isinf(y)))
x = np.compress(mask, x)
y = np.compress(mask, y)

ax = plt.axes(projection=proj)

ax.tricontourf(x, y, dset.t[0,0, mask] )

ax.coastlines()

and produces the following picture. So I guess it works somehow.

This masking should be somehow implemented in cartopy as the bug has been well known since basemap and would avoid any confusion in the future. Also by private conversation with Philip Sommer I understand matplotlib way of doing the triangulation is not the fastest one (see https://github.com/Chilipp/psyplot/issues/6)...It would be nice to implement also his method in cartopy somehow, but I don't know whether this is possible or not.

[guidocioni]

Any updates on this issue 4 years later? :)

This is the current status of cartopy as of 2024. The only way to plot unstructured data is to use PlateCarree because all other projections somehow fail with different problems.

This is the code that I used today

_ = plt.figure(1, figsize=(15, 15))

# proj = ccrs.NearsidePerspective(central_latitude=40,
#                                 central_longitude=-100,
#                                 satellite_height=2e6)

# proj = ccrs.PlateCarree()
# proj = ccrs.Orthographic(central_latitude=40, central_longitude=-100)
# proj = ccrs.Stereographic(central_latitude=40, central_longitude=-100)
# proj = ccrs.Geostationary(
#                                 central_longitude=-100,
#                                 satellite_height=4e6)
proj = ccrs.LambertAzimuthalEqualArea(central_latitude=40,
                                      central_longitude=-100)

ax = plt.axes(projection=proj)
ax.set_extent([-134, -70, 17, 52],
              crs=ccrs.PlateCarree())     

ax.add_feature(cfeature.OCEAN.with_scale('110m'),
                facecolor='#2081C3', zorder=2)
ax.add_feature(cfeature.COASTLINE.with_scale('110m'),
                zorder=2)
ax.add_feature(cfeature.BORDERS.with_scale('110m'),
            zorder=2)
ax.add_feature(cfeature.STATES.with_scale('110m'),
               zorder=2, linewidth=.2)

ax.add_feature(ShapelyFeature(
            Reader('/Users/guidocioni/Downloads/2024eclipse_shapefiles/upath_lo.shp').geometries(),
            ccrs.PlateCarree(),
            edgecolor='yellow',
            facecolor='yellow',
            alpha=0.5),
            zorder=5)

cs = ax.tricontourf(
    ds['clon'],
    ds['clat'],
    ds['prob_cloudy'],
    transform=ccrs.PlateCarree(),
    levels=levels,
    extend='both',
    cmap=cmap)

plt.colorbar(cs, orientation='vertical',
             label='Probability cloud cover > 50% [%]',
                         fraction=0.02, pad=0.03)

and all the different results I get with the projections.

PlateCarree (correct)

Others projections

To me it seems really hard to believe that no solutions was found in all these years. I'm still forced to use basemap in many of my projects because of these small issues that somehow cannot be ironed out.

[rcomer]

It's actually 5 years later ;-)

It looks like you had a working solution at #1293 (comment). Does that also work for your current case? It seems to me that it could be incorporated into Cartopy via a transform_first keyword, as we already have that for contour and contourf (link).

Transformation of filled paths is a much harder problem. We have several open issues about contourf not showing all the levels properly, e.g. #1076. I would be interested to know if your case looks sensible with unfilled contours (i.e. tricontour).

If you would be interested in contributing a pull request based on the "transform first" approach, we would welcome that. For anyone else to work on it, the first thing they would need to do is reproduce the problem. So it would be very helpful if you can supply a self-contained code example that we can run. Note that Cartopy is almost entirely developed and maintained by volunteers, so we cannot make any promises about if/when someone might get to it.