Skip to content

Instantly share code, notes, and snippets.

@hevgyrt
Created January 21, 2020 13:32
Show Gist options
  • Save hevgyrt/c453163a160926a860bc6b796c42aa30 to your computer and use it in GitHub Desktop.
Save hevgyrt/c453163a160926a860bc6b796c42aa30 to your computer and use it in GitHub Desktop.
generic script for plotting AOI in py3 using cartopy
# inspired by https://scitools.org.uk/cartopy/docs/v0.16/gallery/effects_of_the_ellipse.html
import matplotlib.pyplot as plt
import matplotlib
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.patheffects import Stroke
import numpy as np
import shapely.geometry as sgeom
import matplotlib.ticker as mticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
# Creating plot
#fig, ax = plt.subplots(frameon=False,figsize=(20,15),subplot_kw={'projection': ccrs.Mercator()})
#fig, ax = plt.subplots(frameon=False,figsize=(20,15),subplot_kw={'projection': ccrs.AzimuthalEquidistant()})
fig, ax = plt.subplots(frameon=False,figsize=(15,15),subplot_kw={'projection': ccrs.Mercator()})
geodetic = ccrs.Geodetic(globe=ccrs.Globe(datum='WGS84'))
# Deciding extent
#lonmin, lonmax, latmin, latmax = 1,21,57,73
lonmin, lonmax, latmin, latmax = 9.6,15.1,67.1,69
ax.set_extent([lonmin, lonmax, latmin, latmax],geodetic) #x0, x1, y0, y1) of the map in the given coordinate system.
# adding coastline
lscale = 'high' # Scale for coasline
# Add the Stamen aerial imagery at zoom level 7.
#tiler = Stamen('terrain-background')
#ax.add_image(tiler, 7)
#fig = plt.figure()
#ax = fig.add_subplot(1, 1, 1, projection=tiler.crs)
f = cfeature.GSHHSFeature(scale=lscale, levels=[1],
facecolor=cfeature.COLORS['land'])
ax.add_geometries(
f.intersecting_geometries([lonmin, lonmax, latmin, latmax]),
ccrs.PlateCarree(),
facecolor=cfeature.COLORS['land_alt1'],
edgecolor='black')
# Create an inset GeoAxes showing the location of the AOI
sub_ax = fig.add_axes([0.15, 0.64, 0.2, 0.2],
projection=ccrs.AzimuthalEquidistant())
sub_ax.set_extent([-15,35,40,79], geodetic)
# Make a nice border around the inset axes.
effect = Stroke(linewidth=4, foreground='wheat', alpha=0.5)
sub_ax.outline_patch.set_path_effects([effect])
# Add the land, coastlines and the extent of the AOI.
sub_ax.add_feature(cfeature.LAND,facecolor=cfeature.COLORS['land_alt1'])
sub_ax.coastlines(resolution='50m')
extent_box = sgeom.box(lonmin, latmin, lonmax, latmax)
sub_ax.add_geometries([extent_box], ccrs.PlateCarree(), facecolor='none',
edgecolor='blue', linewidth=2)
# Adding references
# in-situ station
ax.plot(12.91313,67.81295,'mo', markersize=5,transform=ccrs.PlateCarree())
mosk_lon, mosk_lat = 12.752111, 67.747667 #island of mosken
vest_lon, vest_lat = 13.8, 67.95 # vestfjorden
ns_lon, ns_lat = 12.0, 68.2 # norwegian sea
#ax.plot(mosk_lon, mosk_lat,'ro', markersize=5,transform=ccrs.PlateCarree())
ax.text(mosk_lon + .03, mosk_lat - .01, 'Mosken',
horizontalalignment='left',
transform=ccrs.Geodetic(),fontsize=12)
ax.text(vest_lon, vest_lat, 'Vestfjorden',
horizontalalignment='left',
transform=ccrs.Geodetic(),fontsize=15,rotation=45)
ax.text(ns_lon, ns_lat, 'Norwegian Sea',
horizontalalignment='center',
transform=ccrs.Geodetic(),fontsize=15,rotation=0)
# Adding ticks
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='gray', alpha=0.3, linestyle='--')
gl.xlabels_top = False
gl.ylabels_left = False
#gl.xlines = False
#gl.xlocator = mticker.FixedLocator([-180, -45, 0, 45, 180])
#gl.xlocator = mticker.FixedLocator(list(np.linspace(lonmin, lonmax,8)))
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.ylabel_style = {'size': 15, 'color': 'black'}
gl.xlabel_style = {'size': 15, 'color': 'black'}
#gl.xlabel_style = {'color': 'red', 'weight': 'bold'}
#plt.show()
#ax.set_xticks(list(np.linspace(lonmin, lonmax,8)), crs=ccrs.PlateCarree())
#ax.set_yticks(list(np.linspace(latmin, latmax)), crs=ccrs.PlateCarree())
# transparent background and boarder
ax.outline_patch.set_visible(False)
ax.background_patch.set_visible(False)
fig.savefig(fname='AOI.png', dpi=150,transparent=True)
#plt.show()
plt.close(fig)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment