Search code examples
pythonshapefilematplotlib-basemap

Fill oceans in high resolution to hide low resolution contours in basemap


When plotting low-resolution contours over a high-resolution coastline I get the following result

enter image description here

I would like to fill the area outside of the coastlines (caused by the low resolution of the underlining filled contour plot) with the ocean color at high resolution.

I tried to use the land-sea mask option without coloring the land

m.drawlsmask(land_color=(0, 0, 0, 0), ocean_color='#2081C3',
resolution='h', lakes=True, zorder=2, grid=1.25)

but the 1.25 resolution is not enough for this level of detail (see second image)

enter image description here

Unfortunately there is no builtin method that fills the ocean (and lakes) with the same resolution used for the coastlines ('h' in my case). As a workaround is there any way to fill the area "outside" of the coastline using the original resolution?

I could use a high resolution land-sea mask in drawlsmask but that's a waste of resource since basemap already has indirectly that information with the polygons given by the coastlines.

General notes:

  1. It looks like other questions on Stack Overflow suggest to use the builtin land sea mask of basemap. I can't because it is too low resolution at this zoom level.
  2. Unfortunately I cannot use Cartopy. I already built my entire pipeline on Cartopy but it is way too slow for what I have to do.

Solution

  • I ended up using the solution posted in Fill oceans in basemap adapted to my needs. Note that, in order to retain the lakes, I had to do multiple passes of fillcontinents, so that's how I did

    # extents contain the projection extents as [lon1, lon2, lat1, lat2]
    m = Basemap(projection='merc',
                    llcrnrlat=extents[2],
                    urcrnrlat=extents[3],
                    llcrnrlon=extents[0],
                    urcrnrlon=extents[1],
                    lat_ts=20,
                    resolution='h')
    
    m.fillcontinents(color='#c5c5c5', lake_color='#acddfe', zorder=1)
    # Fill again the lakes over the contour plot
    m.fillcontinents(color=(0, 0, 0, 0), lake_color='#acddfe', zorder=3)
    
    ax = plt.gca()
    
    # Workaround to add high resolution oceans 
    x0,x1 = ax.get_xlim()
    y0,y1 = ax.get_ylim()
    map_edges = np.array([[x0,y0],[x1,y0],[x1,y1],[x0,y1]])
    
    # getting all polygons used to draw the coastlines of the map
    polys = [p.boundary for p in m.landpolygons]
    polys = [map_edges]+polys[:]
    
    codes = [
       [Path.MOVETO] + [Path.LINETO for p in p[1:]]
         for p in polys
      ]
    
    polys_lin = [v for p in polys for v in p] 
    codes_lin = [c for cs in codes for c in cs]
    path = Path(polys_lin, codes_lin)
    patch = PathPatch(path, facecolor='#acddfe', lw=0, zorder=2)
    
    ax.add_patch(patch)
    
    m.drawcountries(linewidth=0.6)
    m.readshapefile(f'{SHAPEFILES_DIR}/ITA_adm_shp/ITA_adm2',
                            'ITA_adm2', linewidth=0.1, color='gray', zorder=5)
    

    which gives something like this

    enter image description here

    Not perfect (because the shapefile which defines the coastline has a different resolution), but definitely better than before.