[Python] Converting user-drawn polygons on a Cartopy map to a binary mask in Python?

I'm trying to make a program using Tkinter, Cartopy and Matplotlib that would allow users to draw polygons on a map. Once the polygon is closed (right click), I want to convert it to a mask where pixels inside the polygon are set to 1 and pixels outside are set to 0, while keeping the longitude and latitude information. I need those masks to later apply on data simulated over Greenland.

Here's a minimal version of my code:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import tkinter as tk

def on_click(event):
global polygon_closed, pointsX, pointsY

if event.button == 1 and not polygon_closed:
if len(pointsX) > 0:
ax.plot([pointsX[-1], event.xdata], [pointsY[-1], event.ydata], 'r-')

pointsX.append(event.xdata)
pointsY.append(event.ydata)
ax.plot(event.xdata, event.ydata, 'r+')
canvas.draw()

elif event.button == 3 and not polygon_closed:
ax.plot([pointsX[0], pointsX[-1]], [pointsY[0], pointsY[-1]], 'r-')
maskedZone = list(zip(pointsX, pointsY))
polygon = Polygon(maskedZone, closed=True, facecolor='red', edgecolor='black', alpha=0.5)
ax.add_patch(polygon)
canvas.draw()
polygon_closed = True
print(maskedZone)

root = tk.Tk()
root.title("Mask drawing")
root.geometry("500x700")

projection = ccrs.Stereographic(central_latitude=70, central_longitude=-42)
fig, ax = plt.subplots(facecolor="w", subplot_kw=dict(projection=projection))

ax.set_extent([-65, -25, 59, 85], crs=ccrs.PlateCarree())
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.gridlines(draw_labels=True)

canvas = FigureCanvasTkAgg(fig, master=root)
canvas.draw()
canvas.get_tk_widget().pack(fill=tk.BOTH, expand=True)
canvas.mpl_connect("button_press_event", on_click)

pointsX = []
pointsY = []
polygon_closed = False

root.mainloop()


And here is an exemple of the output after closing a polygon:

[(np.float64(38555.16314861295), np.float64(645247.4786037477)),
(np.float64(-126380.37478850619), np.float64(221127.52390829823)),
(np.float64(374316.793949177), np.float64(256470.85346625233))]


I tried creating a shapely.Polygon and then using rasterio.features.rasterize to convert into a mask but the raster ends up being full of just 0 even though there is a shape to rasterize. (Perhaps an output of 10x10 is too small so it is not "viewing" the shape)

import shapely
import rasterio.features

maskedZone = list(zip(pointsX, pointsY))
maskPolygon = shapely.Polygon(maskedZone)
raster = rasterio.features.rasterize([(maskPolygon, 1)], out_shape=(10, 10), all_touched=True)


Question: How can I convert the drawn polygon to a 2D binary mask? I would like to define a regular grid over the same extent as the map, and for each grid point, assign:

• 1 if the point is inside the polygon
• 0 if it is outside

I'm unsure how to handle things properly with Cartopy's projections and coordinate transforms. I read people use rasterio but I'm not sure how to integrate projections correctly.

Any help or direction would be appreciated.

Continue reading...