Search code examples
pythonmatplotlibcartopy

How to draw only the outlines of a bunch of circles

I want to draw the boundary of the intersection of a couple of circles that I created. Here is what I have so far:

EDITED

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.patches import PathPatch
from matplotlib.path import Path
import matplotlib.ticker as mticker
import numpy as np

fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Mercator())

ax.set_extent([11, 38, 56, 70], crs=ccrs.PlateCarree())

ax.add_feature(cfeature.BORDERS, linestyle='-', linewidth=1)
ax.coastlines(linewidth=1)

# Adding the circles
radius = 120 # km
n_samples = 80 # Number of points on the circle
lons = (21.64, 22.5, 23.82, 24.5, 25.44, 26.32, 26.9, 27.11, 27.38, 29.45, 29.8)
lats = (60.13, 61.81, 63.1, 60.56, 62.3, 64.77, 67.14, 60.9, 62.86, 63.84, 61.91)

for lon, lat in zip(lons,lats):
    circle = Geodesic().circle(lon, lat, radius*1000., n_samples=n_samples)
    feature = cartopy.feature.ShapelyFeature([Polygon(circle)], ccrs.Geodetic(), fc='orange', alpha=0.4, ec="black", lw=0.5, linestyle="-")
    ax.add_feature(feature)

I have searched the web, but could not find a solution for this so far. If someone could help me out, that'd be great :)

Solution

  • You don't actually mention what the issue is, so it's difficult the answer.

    But something that you probably unintentionally do at the moment is mixing different units to create your circles. You define your map as being Orthographic, the center center in lat/lon and the radius in meters. What do you expect the result to be?

    Given that you assign the circles having the PlateCarree() projection, the radius should also be defined as such. In your current code everything works fine, but you set the radius to be 250000 degrees!!

    If you for example set the radius to be 2.5 degrees, it looks like the image shown below:

    enter image description here

    You'll have to decide on whether you want a circle to be a circle on the earth (geodetic), on the map projection or on the image (figure coordinates) etc. And then choose the units accordingly in a consistent manner.

    edit:

    The Geodetic class is certainly a good choice, the pyproj package also has functionality for this type of operation.

    The union of several circles can be drawn as shown in the example below, it uses Shapely to perform the union.

    import matplotlib.pyplot as plt

from cartopy.geodesic import Geodesic
import cartopy.crs as ccrs
import cartopy.feature as cfeature

from shapely.geometry import Polygon
from shapely.ops import unary_union

fig, ax = plt.subplots(figsize=(10, 6), subplot_kw=dict(projection=ccrs.Mercator()))

ax.set_extent([11, 38, 56, 70], crs=ccrs.PlateCarree())
ax.add_feature(cfeature.BORDERS, linestyle='-', linewidth=.3)
ax.coastlines(linewidth=.3, resolution="50m")

# Adding the circles
radius = 120 # km
n_samples = 80 # Number of points on the circle
lons = (21.64, 22.5, 23.82, 24.5, 25.44, 26.32, 26.9, 27.11, 27.38, 29.45, 29.8)
lats = (60.13, 61.81, 63.1, 60.56, 62.3, 64.77, 67.14, 60.9, 62.86, 63.84, 61.91)

geo = Geodesic()
circles = [Polygon(geo.circle(lon, lat, radius*1000., n_samples=n_samples)) for lon, lat in zip(lons, lats)]

feature = cfeature.ShapelyFeature(unary_union(circles), ccrs.PlateCarree(), fc='orange', alpha=0.4, ec="k", lw=2, linestyle="-")
# feature = cfeature.ShapelyFeature(circles, ccrs.PlateCarree(), fc='orange', alpha=0.4, ec="k", lw=2, linestyle="-")
ax.add_feature(feature)

    enter image description here

    Without the union, just plotting circles it looks like:

    enter image description here