fig = plt.figure()
# Define the Cartopy Projection
maxlat = 0.
minlat = -50
maxlon = 190
minlon = 100
lat0 = (maxlat + minlat)/2.0
lon0 = (maxlon + minlon)/2.0
proj_kwargs={}
proj_kwargs['central_latitude'] = lat0
proj_kwargs['central_longitude'] = lon0
proj_kwargs['standard_parallels'] = [lat0,lat0]
proj = ccrs.AlbersEqualArea(**proj_kwargs)
map_ax = fig.add_axes([0.1,0.1,0.8,0.8],projection=proj)
x0, y0 = proj.transform_point(lon0, lat0, proj.as_geodetic())
deg2m_lat = 2 * np.pi * 6371 * 1000 / 360
deg2m_lon = deg2m_lat * np.cos(lat0 / 180 * np.pi)
height = (maxlat - minlat) * deg2m_lat
width = (maxlon - minlon) * deg2m_lon
map_ax.set_xlim(x0 - width / 2, x0 + width / 2)
map_ax.set_ylim(y0 - height / 2, y0 + height / 2)
# Plot the Coastlines
map_ax.coastlines()
# Plot the Parallels and Meridians
map_ax.gridlines()
st.merge()
# Plot the source with a dot
map_ax.scatter(lon, lat, marker="o", s=120, zorder=10,
color="k", edgecolor="w", transform = proj.as_geodetic())
# Plot the station with data with a triangle, and the path source-reciever
for station in inv[0]:
if len(st.select(station=station.code))!=0:
map_ax.scatter(station.longitude, station.latitude, marker="v", s=120, zorder=10,
color="k", edgecolor="w", transform = proj.as_geodetic())
plt.plot([station.longitude, lon], [station.latitude, lat], color='0.8', transform=ccrs.Geodetic())
st.select(station=station.code)[0].stats.distance=gps2dist_azimuth(station.latitude, station.longitude, lat, lon, a=6378137.0, f=0.0033528106647474805)[0]