Filling ocean with Natural Earth shapefile does not produce the right result

[guidocioni]

As per https://stackoverflow.com/questions/74433797/fill-oceans-in-high-resolution-to-hide-low-resolution-contours-in-basemap I'm trying to fill the ocean with high resolution polygons to mask out some low resolution filled contours. Unfortunately the built-in land sea mask of basemap is not enough for zoomed in regions.

For this reason I'm trying to fill the ocean shapefiles of natural earth https://www.naturalearthdata.com/downloads/10m-physical-vectors/ using the following code

m = Basemap(projection='merc',
                llcrnrlat=36,
                urcrnrlat=47,
                llcrnrlon=6,
                urcrnrlon=19,
                lat_ts=20,
                resolution='h')

shp = m.readshapefile('../input/shapefiles/ne_10m_ocean/ne_10m_ocean',
                'ne_10m_ocean', drawbounds = True)
patches   = []
for info, shape in zip(m.ne_10m_ocean_info, m.ne_10m_ocean):
    patches.append( Polygon(np.array(shape), True) )

plt.gca().add_collection(PatchCollection(patches, facecolor= 'red', edgecolor='k', linewidths=1., zorder=2))

This, however, does also fill some of the major islands with the same color of the ocean, as the following figure shows.

Am I doing something wrong? I checked the original shapefile and the polygons seem to be well defined. Here is the same plot done with Cartopy

[molinav]

Hi, @guidocioni! First of all, sorry for replying so late. I have downloaded the ne_10m_ocean.zip file and for the moment I can confirm you that I get into the same problem as you.

Since you have tried with cartopy and it is plotted correctly, I would think that something might be wrong in Basemap.readshapefile that makes some polygons remain open and their inside is therefore considered part of the background. I will need to check a bit deeper.

[molinav]

By the way, this error seems to be around for some time already, see this other very similar issue: https://github.com/matplotlib/basemap/issues/485

[molinav]

I was diving a bit more into the problem. Can you try to reproduce the following snippet? I am basically doing the same as in your snippet but, instead of centering the map in Italy, I move the map boundaries so that they cover the 180-degree meridian:

from matplotlib.patches import Polygon
from mpl_toolkits.basemap import Basemap

# Create Basemap instance.
plt.clf()
ax = plt.gca()
bmap = Basemap(projection="merc",
               llcrnrlat=40,
               urcrnrlat=80,
               llcrnrlon=120,
               urcrnrlon=200,
               lat_ts=0,
               resolution="c")
bmap.drawmapboundary()

# Collect polygons within the map boundary.
shppath = "../input/shapefiles/ne_10m_ocean/ne_10m_ocean"
bmap.readshapefile(shppath, "ne_10m_ocean", drawbounds=False)
patches = []
for info, verts_i in zip(bmap.ne_10m_ocean_info, bmap.ne_10m_ocean):
    verts_i = np.array(verts_i)
    if ((verts_i[:, 0] >= bmap.llcrnrx) & (verts_i[:, 0] <= bmap.urcrnrx)).any():
        if ((verts_i[:, 1] >= bmap.llcrnry) & (verts_i[:, 1] <= bmap.urcrnry)).any():
            patch = Polygon(verts_i, closed=True, facecolor="yellow",
                            edgecolor="black", linewidth=1, zorder=2)
            patches.append(patch)

# Draw and fill only the first polygon.
for patch in patches[:1]:
    ax.add_patch(patch)

And this is what I get:

I am plotting only the first polygon, with fill color set to yellow. The first polygon here is the main polygon for Eurafrasia, but if you look at the 180-degree meridian, the polygon is actually being closed following the oceans and not the land bodies, so the "inside" of the polygon is actually the ocean, and therefore it is the ocean what it gets filled with the facecolor attribute.

[molinav]

The same snippet if I use "ne_10m_land" instead of "ne_10m_ocean" (from the same website):

from matplotlib.patches import Polygon
from mpl_toolkits.basemap import Basemap

# Create Basemap instance.
plt.clf()
ax = plt.gca()
bmap = Basemap(projection="merc",
               llcrnrlat=40,
               urcrnrlat=80,
               llcrnrlon=120,
               urcrnrlon=200,
               lat_ts=0,
               resolution="c")
bmap.drawmapboundary()

# Collect polygons within the map boundary.
shppath = "../input/shapefiles/ne_10m_land/ne_10m_land"
bmap.readshapefile(shppath, "ne_10m_land", drawbounds=False)
patches = []
for info, verts_i in zip(bmap.ne_10m_land_info, bmap.ne_10m_land):
    verts_i = np.array(verts_i)
    if ((verts_i[:, 0] >= bmap.llcrnrx) & (verts_i[:, 0] <= bmap.urcrnrx)).any():
        if ((verts_i[:, 1] >= bmap.llcrnry) & (verts_i[:, 1] <= bmap.urcrnry)).any():
            patch = Polygon(verts_i, closed=True, facecolor="yellow",
                            edgecolor="black", linewidth=1, zorder=2)
            patches.append(patch)

# Draw and fill only the first polygon.
for patch in patches[:1]:
    ax.add_patch(patch)

And this is what I get:

[molinav]

Finally, to reproduce the example in your original message, I create the Basemap instance, I draw the map boundary with fill color set to red, and then I add the polygons from "ne_10m_land" with facecolor set to white:

from matplotlib.collections import PatchCollection
from matplotlib.patches import Polygon
from mpl_toolkits.basemap import Basemap

# Create Basemap instance.
plt.clf()
ax = plt.gca()
bmap = Basemap(projection="merc",
               llcrnrlat=36,
               urcrnrlat=47,
               llcrnrlon=6,
               urcrnrlon=19,
               lat_ts=20,
               resolution="c")
bmap.drawmapboundary(fill_color="red")

# Collect all polygons.
shppath = "../input/shapefiles/ne_10m_land/ne_10m_land"
bmap.readshapefile(shppath, "ne_10m_land", drawbounds=False)
patches = []
patches_kwds = dict(facecolor="white", edgecolor="black", linewidths=1, zorder=2)
for info, verts_i in zip(bmap.ne_10m_land_info, bmap.ne_10m_land):
    patches.append(Polygon(np.array(verts_i), closed=True))

# Draw and fill all polygons.
ax.add_collection(PatchCollection(patches, **patches_kwds))

And then I get the same result as in your example with cartopy:

I would be curious to know about the snippet that you used for cartopy, because if you also used the "ne_10m_ocean" shapefile with cartopy, then they must have some additional magic on how to close the polygons so that they close the land instead of the ocean.

[guidocioni]

Sorry, honestly this was a long time ago and I don't remember anymore if (and how) I fixed it. Judging from the stack overflow issue I opened (https://stackoverflow.com/questions/74433797/fill-oceans-in-high-resolution-to-hide-low-resolution-contours-in-basemap), it was working at the end...so I'm not sure how I fixed it. Or maybe I just ignored it.

What is the difference between the last example you provided and what I was trying at the beginning? Because if it's working with basemap that we can just take this last example as the way to go and that's fine.

[molinav]

The only relevant change in my snippet is the shapefile downloaded from the Natural Earth website.

The original message was using "ne_10m_ocean". Inspecting this shapefile, it defines ocean polygons (i.e. the polygons are defined in a way that the water bodies are the inside of the polygons). So when using facecolor="red" for this shapefile, one is actually requesting to fill the inside of the polygons (i.e. the water bodies) in red.

My last snippet is using "ne_10m_land" from the same website. In this shapefile, the polygons are defined for the land bodies (i.e. the inside of the polygons is the land). So when using facecolor with this shapefile, facecolor now is the color for the lands.

To reproduce your goal image, I combined a fill_color="red" for Basemap.drawmapboundary (the fill color for the water) and a facecolor="white" for the land polygons (the fill color for the land).

PS: it is totally understandable that you do not remember anymore. I try to address the issues when I find some spare time, unfortunately this sometimes comes with a lot of delay. Sorry!

[guidocioni]

The only relevant change in my snippet is the shapefile downloaded from the Natural Earth website.

The original message was using "ne_10m_ocean". Inspecting this shapefile, it defines ocean polygons (i.e. the polygons are defined in a way that the water bodies are the inside of the polygons). So when using facecolor="red" for this shapefile, one is actually requesting to fill the inside of the polygons (i.e. the water bodies) in red.

My last snippet is using "ne_10m_land" from the same website. In this shapefile, the polygons are defined for the land bodies (i.e. the inside of the polygons is the land). So when using facecolor with this shapefile, facecolor now is the color for the lands.

To reproduce your goal image, I combined a fill_color="red" for Basemap.drawmapboundary (the fill color for the water) and a facecolor="white" for the land polygons (the fill color for the land).

PS: it is totally understandable that you do not remember anymore. I try to address the issues when I find some spare time, unfortunately this sometimes comes with a lot of delay. Sorry!

Thanks for summing up. No worries about the delay, it's just good enough that you're finding the time to maintain this repo.

Then I'd say we can close this issue because you're clearly proposing a solution. What do you think?

I'd have to check again what I'm doing right now in my repo to avoid this issue but clearly it seems I did find a workaround...just don't remember which one :D

[molinav]

Is this repo publicly available? At least I would be curious to know about how you solved your problem (at that time).

I think the issue can be closed, Basemap.readshapefile seems to be doing what is supposed to do; it is more a matter of selecting the shapefile that best fits the plotting goals.

[guidocioni]

Is this repo publicly available? At least I would be curious to know about how you solved your problem (at that time).

I think the issue can be closed, Basemap.readshapefile seems to be doing what is supposed to do; it is more a matter of selecting the shapefile that best fits the plotting goals.

Unfortunately not, but it seems what I used in the end was the solution I described in the Stackoverflow post. Going to close for now