Skip to content

Conversation

@htonchia
Copy link
Contributor

Commented out the part plotting the overlapping cells (despite a smaller zorder).
Now the overlapping cells larger than half the projection span stay masked.

Rationale

There is already a patch in pcolormesh that masks overlapping cells which are the cells crossing the longitude boundary and therefore overlapping the plot.
In that patch, the overlapping cells, those who size is more than half the projection span in x, are masked and plot then again with a zorder supposedly small enough so that they don't appear.
Unfortunately they still appeared on the plot.
The idea in plotting the bad cells is that someone may want to use them, by splitting them at the boundary and plot the split part. But they shouldn't be plotted as there are.

So the fix is to comment out the part plotting the bad cells.

This patch works only for pcolormesh, there should be one of the same kind to mask the data generating overlapping cells or contours in contourf, pcolor and contour.

This patch removes the cell greater than half the projection span in x, which works well for Plate Carree but not for the narrow part of a projection like Mollweide. This part is not fixed.

The image below displays ice concentration provided in a polar stereographic grid.
It had been plotted after the proposed fix.
The left plots are pcolormesh which includes the patch for overlapping cells, the right plots are contourf and displays many overlapping cells.
The pcolormesh plot in Mollweide projection still displays some small overlapping cells (less than half the projection span) on the narrow part of the plot.

There are no overlapping cells in the stereographic projection.
pcolormesh_contourf

The data can be found here.

import matplotlib.pyplot as plt
from netCDF4 import Dataset
import cartopy.crs as ccrs

# file available here:
# https://github.com/htonchia/cartopy-plotting-examples/blob/master\
#/ice_conc_nh_polstere-100_multi_201912091200.nc

FILE = "ice_conc_nh_polstere-100_multi_201912091200.nc"
VARNAME = 'ice_conc'
STEP = 1  

def projection4test():
    """
    define the list of projection for the tests
    """
#    proj_data = ccrs.Stereographic(
#        central_longitude=-45, central_latitude=90,
#        true_scale_latitude=70)

    proj2s = {}
    proj2s['PlateCarre_90'] = ccrs.PlateCarree(central_longitude=90)
    proj2s['Robinson_0'] = ccrs.Robinson(central_longitude=0)
    proj2s['Mollweide_90'] = ccrs.Mollweide(central_longitude=90)
    proj2s['polar 45'] = ccrs.Stereographic(
        central_longitude=45, central_latitude=90,
        true_scale_latitude=70)
    
    return proj2s

def main(data, varname, step=1):
    """
    Run plotting tests on different set of projections
    """

    lats = data.variables['lat'][::step, ::step]
    lons = data.variables['lon'][::step, ::step]

    projs = projection4test()

    dplotkwargs = {}
    dplotkwargs['pcolor'] = {'cmap' : 'bwr'}
    dplotkwargs['pcolormesh'] = {'cmap' : 'bwr'}
    dplotkwargs['contourf'] = {'cmap' : 'bwr'}
    dplotkwargs['contour'] = {'colors' : 'black', 'linewidths' : 1}

    opes = ['pcolormesh', 'contourf']

    plt.figure(figsize=(12, 6.4), dpi=100)
    plt.subplots_adjust(wspace=0.01, left=0.01, right=0.99, top=0.95,
                        hspace=0.01, bottom=0.01)

    iii = 0
    
    plt.suptitle("{} {}".format(*opes))
    for proj in projs:
        for ope in opes:
            iii += 1
            mdata = data.variables[varname][0, ::step, ::step]
            axe = plt.subplot(4, 2, iii,
                              projection=projs[proj])
            axe.set_title("{} {}".format(proj, ope))
            if (hasattr(projs[proj], 'proj4_params') and
                    projs[proj].proj4_params['proj'] in
                    ['geos', 'nsper', 'gnom']
                    and ope == 'pcolormesh'):
                # will crash if we try to plot
                axe.coastlines()
                continue
            xy2 = axe.projection.transform_points(
                    ccrs.Geodetic(),
                    lons, lats)
            cx2, cy2 = xy2[..., 0], xy2[..., 1]
            getattr(axe,ope)(cx2, cy2, mdata, **dplotkwargs[ope])
            axe.coastlines()

    plt.savefig("img/{}_{}.png".format(*opes),
                    dpi=200, pad_inches=0)


if __name__ == '__main__':
    main(Dataset(FILE), varname=VARNAME, step=STEP)


Implications

I don't see any implication.

Commented out the part plotting the overlapping cells (despite a smaller zorder).
Now the overlapping cells larger than half the projection span stay masked.
# this method
collection._wrapped_collection_fix = pcolor_col

# the part below in comment retrieves the masked

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

W291 trailing whitespace

collection._wrapped_collection_fix = pcolor_col

# the part below in comment retrieves the masked
# overlapping cells and plot them with pcolor and a

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

W291 trailing whitespace

# the part below in comment retrieves the masked
# overlapping cells and plot them with pcolor and a
# smaller zorder to hide them
# but there are not always hidden

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

W291 trailing whitespace

# overlapping cells and plot them with pcolor and a
# smaller zorder to hide them
# but there are not always hidden
# if one wants to retrieve those cells, split them at the

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

W291 trailing whitespace

# smaller zorder to hide them
# but there are not always hidden
# if one wants to retrieve those cells, split them at the
# boundary and plot them, the code below could be

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

W291 trailing whitespace

# if one wants to retrieve those cells, split them at the
# boundary and plot them, the code below could be
# useful.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

W293 blank line contains whitespace

# pcolor_data = np.ma.array(C, mask=~mask)
#
# pts = coords.reshape((Ny, Nx, 2))
#

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

W291 trailing whitespace

fix trailing white spaces
@dopplershift
Copy link
Contributor

We shouldn't be commenting out code, but just remove it if it's wrong. Before doing so, I'd really want to see some before and after comparisons, because that original code was added for a reason (that I haven't dug into yet to figure out).

@htonchia
Copy link
Contributor Author

The part, I have commented out, plots with pcolor the overlapping cells that have been masked earlier. It plots them with a lower zorder. So if there was no land_mask or it the data were atmospheric data, they would not show. With oceanographic data, like ice concentration, those cells still appear on the plot despite their lower zorder.
It seems to me that its is a very common issue when one designs plot functions for atmospheric data and uses them for oceanographic data.

My understanding is that by keeping those cells, one can later process them to split them at the boundary and plot them correctly, that's why I kept the code commented out. Removing it would be cleaner and if one wants to process those cells, then one would have to change the code.

Copy link
Contributor

@greglucas greglucas left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems like this should be investigated a little more as to why the code was there in the first place.

Is the problem that these masked cells get passed through everything and then later on they get split into two cells (one on each side of the boundary) and only one of those cells keeps the mask, or it somehow changes the masking of the cell?

It seems like there is a more fundamental underlying problem here and that the masks should work in this situation and that would be worthwhile to look into. I'm thumbs down on this brute force comment out the block of code as it currently stands.

# useful.

# #############
#
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What about making this a kwarg option, rather than commenting it out?

if kwargs['keep_bad_cells']:
    ...

# if np.any(~pcolor_data.mask):
# # plot with slightly lower zorder to avoid odd issue
# # where the main plot is obscured
# zorder = collection.zorder - .1
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did you try making this a much lower zorder like -1, so it is below all other artists?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I tried a much lower order like -1.
In that case, it doesn't always work so I tried an ever larger number, and it worked.
Now I think there is no guarantee that it will always hid the cell, not plotting them is a better guarantee.
I think the issue comes with masked data, like oceanographic data with a land mask. I would say that when this was written, pushing the "hidden cell" behind the good one, was good enough, but when there is no data on the land mask, then the cells behind still appear on the masked area, not obscuring the data but the mask, and still making a mess of the plot.

Making it a kwargs is a possible option.

The idea was that those cell had to be hidden, so let's hide them efficiently. I admit that I spend quite some time thinking about it and I am sure of one thing is that cartopy needs to get rid of the cells...

@greglucas
Copy link
Contributor

I just took a little closer look at this, and I think it is important to keep those lines in there for boundary wrapping purposes. For your cases could you use mesh._wrapped_collection_fix.set_visible(False)?

In general, though, I think those cells need to be replotted to appear on the boundaries as demonstrated below.
Figure_1

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

fig, (ax1, ax2) = plt.subplots(nrows=2, subplot_kw={'projection': ccrs.Mercator()})

x = np.linspace(0, 360)
y = np.linspace(-45, 45)
X, Y = np.meshgrid(x, y)
Z = np.ma.masked_greater(X**2 + Y**2, 200**2)

mesh1 = ax1.pcolormesh(X, Y, Z, transform=ccrs.PlateCarree())
mesh2 = ax2.pcolormesh(X, Y, Z, transform=ccrs.PlateCarree())
mesh2._wrapped_collection_fix.set_visible(False)

ax1.set_title('Before')
ax2.set_title('After')

plt.show()

@htonchia
Copy link
Contributor Author

htonchia commented Feb 18, 2020 via email

@greglucas
Copy link
Contributor

I apologize, I wasn't very clear in the previous message. I think the "Before" in the image above is the correct way of doing things. In the second "After" image, the cells that cross the boundary are masked and removed from the plot (same as commenting out the lines of code).

Note that I haven't looked into the full transformation path, but I think what is happening is:

  • Make a quadmesh and mask out any cells spanning half the plot.
  • Replot the masked quadmesh cells from above with pcolor, which returns a PolyCollection
  • Now, here is the speculation on my part, I think that the Polygons in the PolyCollection go through the _project_geometry code path which detects boundary intersections and makes two separate Polygons, one on each side of the boundary.

That is why I think this code path is important to leave in, we want those two extra Polygons to go right up to the border on each side of the map, not to just ignore them.

Perhaps this means that the boundary needs to be defined at a higher resolution for the actual clip paths being used? I'm not sure if #1422 will help here with the clipping or not.
One of the worst offenders of your plot demonstrations is the Mollweide, which it looks like has some issues at the poles. See #1333 and the fix of #1334

I think that your other PR #1420 will help with identifying which cells to mask/wrap better. But, I think there is likely a better way to deal with the masked/wrapped cells than just dropping them altogether.

@htonchia
Copy link
Contributor Author

Thanks for your comment.

I have been studying the issue and will issue soon (I hope) a detailed comment.
I have learned a little bit more on the subject.
I didn't understand why your example didn't work as the data used that highlighted the issue to fix, and now I think I have part of the answer.

The issue of overlapping cells not managed properly by pcolor comes from data provided in lat lon for which we don't know the projection used to create the grid.
So I am working on a fix that would check this case and manage it.

Your comment and your example have been very useful.

@htonchia
Copy link
Contributor Author

This PR is closed because taken into account by #1467

After further analysis and using both my example (ice data in stereographic projection) and the example above, I came to the conclusion that pcolor was failing to plot the overlapping cells when the data coordinates were given in the axis projection.

Here is the code to used to show that point.

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs


PROJ_DATA = ccrs.Stereographic(
        central_longitude=-45, central_latitude=90,
        true_scale_latitude=70)

def sample_data(proj=PROJ_DATA):
    """
    simulated data in a mesh in projection PROJ_DATA
    with coordinates in latitudes and longitudes
    return (lons, lats) in  projection proj, simulated data, proj
    """
    xcc = np.linspace(-3845, 3745, 760//10) * 1000
    ycc = np.linspace(5845, -5345, 1120//10) * 1000
    xcm, ycm = np.meshgrid(xcc, ycc)
    lls = ccrs.Geodetic().transform_points(proj, xcm, ycm)
    lats = np.deg2rad(lls[..., 1])
    lons = np.deg2rad(lls[..., 0])
    mlats = np.min(lats)
    lats = lats - mlats
    wave = 0.75 * (np.sin(2 * lats) ** 8) * np.cos(4 * lons)
    mean = 0.5 * np.cos(2 * lats) * ((np.sin(2 * lats)) ** 2 + 2)
    data = wave + mean
    circle = (xcm+2000000)**2 + (ycm-200000)**2
    data = np.ma.masked_where(circle < 1000000**2, data)
    return (lls[..., 0], lls[..., 1]), data, proj

def main():
    """
    Test geoaxe pcolormesh with cells with not removed by
    cell border length but by the diagonals length
    """
    coords, data2plot, proj_data = sample_data(proj=PROJ_DATA)
    proj = ccrs.PlateCarree(central_longitude=182)
    plt.figure(figsize=(12, 6.4), dpi=100)

    ### here we consider that we don't know the data projection
    ### and we have the coordinates in lat lon
    xy2 = proj.transform_points(
         ccrs.Geodetic(),
         coords[0], coords[1])
    cx2, cy2 = xy2[..., 0], xy2[..., 1]

    axe0 = plt.subplot(4, 2, 1, projection=proj)
    axe0.pcolormesh(cx2, cy2, data2plot)
    axe0.set_title('pcolormesh _wrapped_collection_fix: auto, ax proj')
    axe1 = plt.subplot(4, 2, 3, projection=proj)
    axe1.set_title('pcolormesh _wrapped_collection_fix visible ax proj')
    mesh1 = axe1.pcolormesh(cx2, cy2, data2plot)
    axe2 = plt.subplot(4, 2, 5, projection=proj)
    mesh2 = axe2.pcolormesh(cx2, cy2, data2plot)
    axe2.set_title('pcolormesh _wrapped_collection_fix not visible ax proj')
    axe3 = plt.subplot(4, 2, 7, projection=proj)
    axe3.pcolor(cx2, cy2, data2plot)
    if hasattr(mesh2, '_wrapped_collection_fix'):
        mesh2._wrapped_collection_fix.set_visible(False)
        mesh1._wrapped_collection_fix.set_visible(True)
    axe3.set_title('pcolor coords ax proj')

    ### here we consider that we know the data projection
    ### and we have the coordinates in data projection
    xy2 = proj_data.transform_points(
         ccrs.Geodetic(),
         coords[0], coords[1])
    cx2, cy2 = xy2[..., 0], xy2[..., 1]

    axe10 = plt.subplot(4, 2, 2, projection=proj)
    axe10.pcolormesh(cx2, cy2, data2plot, transform=proj_data)
    axe10.set_title('pcolormesh _wrapped_collection_fix auto projdata')
    axe4 = plt.subplot(4, 2, 4, projection=proj)
    mesh4 = axe4.pcolormesh(cx2, cy2, data2plot, transform=proj_data)
    axe4.set_title('pcolormesh _wrapped_collection_fix visible projdata')
    axe5 = plt.subplot(4, 2, 6, projection=proj)
    mesh5 = axe5.pcolormesh(cx2, cy2, data2plot, transform=proj_data)
    axe5.set_title('pcolormesh _wrapped_collection_fix not visible  projdata')
    axe6 = plt.subplot(4, 2, 8, projection=proj)
    axe6.pcolor(cx2, cy2, data2plot, transform=proj_data)
    if hasattr(mesh5, '_wrapped_collection_fix'):
        mesh5._wrapped_collection_fix.set_visible(False)
        mesh4._wrapped_collection_fix.set_visible(True)
    axe6.set_title('pcolor coords in projdata')


    plt.show()
    plt.savefig("img/PR1419_test_wrapped_collection.png",
                    dpi=100, pad_inches=0)



if __name__ == '__main__':
    main()

Here is the resulting plot:
PR1419_test_wrapped_collection

Here is a very simple example with big cells showing how it works.
This example uses the geoaxes of #1467.
In axes projection pcolor plot the cells as they come, not managing wrapping, pcolormesh hides the large cells considering it is an overlap (they can be set to visible though as showed in the previous code).

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

axproj = ccrs.Mercator()

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2, 
         figsize=(12, 6.4), dpi=100,
         subplot_kw={'projection': axproj})

lon = np.linspace(120, 240, 4)
lat = np.linspace(0, 75, 4)
lons, lats = np.meshgrid(lon, lat)
Z = lons**2 + lats**2 * 10

xyas = axproj.transform_points(ccrs.Geodetic(), lons, lats)

ax1.pcolormesh(lon, lat, Z, transform=ccrs.PlateCarree())
ax2.pcolormesh(xyas[..., 0], xyas[..., 1],  Z)
ax3.pcolor(lon, lat, Z, transform=ccrs.PlateCarree())
ax4.pcolor(xyas[..., 0], xyas[..., 1], Z)

ax1.set_title('pcolormesh coords in data projection')
ax2.set_title('pcolormesh coords in axis projection')
ax3.set_title('pcolor coords in data projection')
ax4.set_title('pcolor coords in axis projection')
for axe in [ax1, ax2, ax3, ax4]:
    axe.set_global()
    axe.coastlines()

plt.show()
plt.savefig("img/I1416_test_dataproj.png",
                    dpi=100, pad_inches=0)

I1416_test_dataproj

@htonchia htonchia closed this Feb 18, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants