How to plot woven with potatoes and matplotlib?
For building skymaps, I just switched from Basemap to map, I like it a lot more.
(The main reason was segfaulting Basemap on some computers, which I couldn't fix).
The only thing I am struggling with is getting the cloth circle (used to display the cone of our telescope).
This is a sample code showing random stars (I'm using a directory for the real thing):
import matplotlib.pyplot as plt
from cartopy import crs
import numpy as np
# create some random stars:
n_stars = 100
azimuth = np.random.uniform(0, 360, n_stars)
altitude = np.random.uniform(75, 90, n_stars)
brightness = np.random.normal(8, 2, n_stars)
fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection=crs.NorthPolarStereo())
ax.background_patch.set_facecolor('black')
ax.set_extent([-180, 180, 75, 90], crs.PlateCarree())
plot = ax.scatter(
azimuth,
altitude,
c=brightness,
s=0.5*(-brightness + brightness.max())**2,
transform=crs.PlateCarree(),
cmap='gray_r',
)
plt.show()
How would I add a woven circle with a specific radius in degrees to this image? https://en.wikipedia.org/wiki/Tissot%27s_indicatrix
I mean to go back and add two functions from GeographicLib that provide forward and backward geodesic computations, while this is just a matter of calculating the geodesic circle by sampling at the appropriate azimuths for a given lat / lon / radius value.Alas, I haven't done that yet. but there is a pretty primitive (but efficient) wrapper in pyproj function.
To implement the cloth indicatrix, then the code might look something like this:
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
from pyproj import Geod
import shapely.geometry as sgeom
def circle(geod, lon, lat, radius, n_samples=360):
"""
Return the coordinates of a geodetic circle of a given
radius about a lon/lat point.
Radius is in meters in the geodetic coordinate system.
"""
lons, lats, back_azim = geod.fwd(np.repeat(lon, n_samples),
np.repeat(lat, n_samples),
np.linspace(360, 0, n_samples),
np.repeat(radius, n_samples),
radians=False,
)
return lons, lats
def main():
ax = plt.axes(projection=ccrs.Robinson())
ax.coastlines()
geod = Geod(ellps='WGS84')
radius_km = 500
n_samples = 80
geoms = []
for lat in np.linspace(-80, 80, 10):
for lon in np.linspace(-180, 180, 7, endpoint=False):
lons, lats = circle(geod, lon, lat, radius_km * 1e3, n_samples)
geoms.append(sgeom.Polygon(zip(lons, lats)))
ax.add_geometries(geoms, ccrs.Geodetic(), facecolor='blue', alpha=0.7)
plt.show()
if __name__ == '__main__':
main()
