PK 8KX) plot_hydrosheds.py# -*- coding: utf-8 -*-
"""
===================
Plotting shapefiles
===================
Put some colors and labels on shapefiles
In this script, we use data from the `HydroSHEDS `_
database to illustrate some functionalities of salem Maps. The data shows the
sub-basins of the Nam Co Lake catchment in Tibet. We navigate between the
various tributary catchments of the lake.
"""
import salem
import matplotlib.pyplot as plt
# read the shapefile
shpf = salem.get_demo_file('Lev_09_MAIN_BAS_4099000881.shp')
gdf = salem.read_shapefile(shpf)
# Get the google map which encompasses all geometries
g = salem.GoogleVisibleMap(x=[gdf.min_x.min(), gdf.max_x.max()],
y=[gdf.min_y.min(), gdf.max_y.max()],
maptype='satellite', scale=2,
size_x=400, size_y=400)
ggl_img = g.get_vardata()
# Get each level draining into the lake, then into the last level, and so on
gds = []
prev_id = [gdf.iloc[0].MAIN_BAS]
while True:
gd = gdf.loc[gdf.NEXT_DOWN.isin(prev_id)]
if len(gd) == 0:
break
gds.append(gd)
prev_id = gd.HYBAS_ID.unique()
# make a map of the same size as the image
sm = salem.Map(g.grid, factor=1)
sm.set_rgb(ggl_img) # add the background rgb image
# add all the draining basins
cmap = plt.get_cmap('Blues')
for i, gd in enumerate(gds):
# here we use a trick. set_shapefile uses PatchCollections internally,
# which is fast but does not support legend labels.
# so we use set_geometry instead:
for g, geo in enumerate(gd.geometry):
# we don't want more than one label per level
label = 'Level {:02d}'.format(i+1) if g == 0 else None
sm.set_geometry(geo, facecolor=cmap(i/(len(gds)-1)),
alpha=0.8, label=label)
# Get the polygon of the last sink (i.e. the lake) and plot it
gds_0 = gdf.loc[gdf.HYBAS_ID == gdf.iloc[0].MAIN_BAS]
sm.set_shapefile(gds_0, linewidth=2)
# Compute the outline of the entire basin and plot it
gds_1 = gdf.geometry.unary_union
sm.set_geometry(gds_1, linewidth=4)
# plot!
f, ax = plt.subplots(figsize=(6, 4))
ax.set_position([0.05, 0.06, 0.7, 0.9])
sm.visualize(ax=ax)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.show()
PK !KX # plot_polarstereo.py# -*- coding: utf-8 -*-
"""
===============================
Polar stereographic projections
===============================
This example illustrates that salem graphics have a small issue at the pole
and proposes an alternative based on cartopy.
The only problem at the pole is the plotting of the longitude contours.
There is probably no easy way so solve this problem in salem,
hence the suggestion to use cartopy in this case. Note that if the pole isn't
in the domain, salem works just fine.
"""
from salem import open_xr_dataset, get_demo_file
import matplotlib.pyplot as plt
# prepare the figure
f = plt.figure(figsize=(8, 7))
# WRF polar file(s)
d1 = open_xr_dataset(get_demo_file('geo_em_d01_polarstereo.nc'))
d2 = open_xr_dataset(get_demo_file('geo_em_d02_polarstereo.nc'))
# Plot with salem
ax = plt.subplot(2, 2, 1)
d1.HGT_M.salem.quick_map(ax=ax, cmap='Oranges')
ax = plt.subplot(2, 2, 3)
d2.HGT_M.salem.quick_map(ax=ax, cmap='Oranges')
# Now with cartopy
proj = d1.salem.cartopy()
ax = plt.subplot(2, 2, 2, projection=proj)
ax.coastlines()
ax.gridlines()
d1.HGT_M.plot(ax=ax, transform=proj, cmap='Oranges')
ax.set_extent(d1.salem.grid.extent, crs=proj)
# D2 can use a higher resolution coastline
proj = d2.salem.cartopy()
ax = plt.subplot(2, 2, 4, projection=proj)
ax.coastlines(resolution='50m')
ax.gridlines()
d2.HGT_M.plot(ax=ax, transform=proj, cmap='Oranges')
ax.set_extent(d2.salem.grid.extent, crs=proj)
plt.show()
PK .KXoܝ6 6 plot_googlestatic.py# -*- coding: utf-8 -*-
"""
===============================
Plot on a google map background
===============================
Google static maps API
In this script, we use `motionless `_
to download an image from the `google static map API `_
and plot it on a :py:class:`~salem.Map`. We then add information to the map
such as a glacier outline (from the `RGI `_),
and ground penetrating radar measurements (GPR, from the
`GlaThiDa `_ database).
The GPR measurements were realized in 1997, the glacier outlines are from
2003, and the map background is from 2016. This illustrates the retreat of
the Kesselwandferner glacier.
"""
import numpy as np
import pandas as pd
import salem
from salem import get_demo_file, DataLevels, GoogleVisibleMap, Map
import matplotlib.pyplot as plt
# prepare the figure
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
# read the shapefile and use its extent to define a ideally sized map
shp = salem.read_shapefile(get_demo_file('rgi_kesselwand.shp'))
# I you need to do a lot of maps you might want
# to use an API key and set it here with key='YOUR_API_KEY'
g = GoogleVisibleMap(x=[shp.min_x, shp.max_x], y=[shp.min_y, shp.max_y],
scale=2, # scale is for more details
maptype='satellite') # try out also: 'terrain'
# the google static image is a standard rgb image
ggl_img = g.get_vardata()
ax1.imshow(ggl_img)
ax1.set_title('Google static map')
# make a map of the same size as the image (no country borders)
sm = Map(g.grid, factor=1, countries=False)
sm.set_shapefile(shp) # add the glacier outlines
sm.set_rgb(ggl_img) # add the background rgb image
sm.set_scale_bar(location=(0.88, 0.94)) # add scale
sm.visualize(ax=ax2) # plot it
ax2.set_title('GPR measurements')
# read the point GPR data and add them to the plot
df = pd.read_csv(get_demo_file('gtd_ttt_kesselwand.csv'))
dl = DataLevels(df.THICKNESS, levels=np.arange(10, 201, 10), extend='both')
x, y = sm.grid.transform(df.POINT_LON.values, df.POINT_LAT.values)
ax2.scatter(x, y, color=dl.to_rgb(), s=50, edgecolors='k', linewidths=1)
dl.append_colorbar(ax2, label='Ice thickness (m)')
# make it nice
plt.tight_layout()
plt.show()
PK /KXH H plot_landsea_mask.py# -*- coding: utf-8 -*-
"""
===============
Shape to raster
===============
Compute a land/sea mask from a shapefile.
This example illustrates the two different methods available to compute a
raster mask from shapefile polygons.
"""
import salem
import matplotlib.pyplot as plt
# make a local grid from which we will compute the mask
# we make it coarse so that we see the grid points
grid = salem.Grid(proj=salem.wgs84, x0y0=(-18, 3), nxny=(25, 15), dxdy=(1, 1))
# read the ocean shapefile (data from http://www.naturalearthdata.com)
oceans = salem.read_shapefile(salem.get_demo_file('ne_50m_ocean.shp'),
cached=True)
# read the lake shapefile (data from http://www.naturalearthdata.com)
lakes = salem.read_shapefile(salem.get_demo_file('ne_50m_lakes.shp'),
cached=True)
# The default is to keep only the pixels which center is within the polygon:
mask_default = grid.region_of_interest(shape=oceans)
mask_default = grid.region_of_interest(shape=lakes, roi=mask_default)
# But we can also compute a mask from all touched pixels
mask_all_touched = grid.region_of_interest(shape=oceans, all_touched=True)
mask_all_touched = grid.region_of_interest(shape=lakes, all_touched=True,
roi=mask_all_touched)
# Make a map to check our results
sm = salem.Map(grid, countries=False)
sm.set_shapefile(oceans, edgecolor='k', facecolor='none', linewidth=2)
sm.set_shapefile(lakes, edgecolor='k', facecolor='none', linewidth=2)
sm.set_plot_params(cmap='Blues', vmax=2)
# prepare the figure
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3))
# plot 1
sm.set_data(mask_default)
sm.visualize(ax=ax1, addcbar=False, title='Default')
# plot 2
sm.set_data(mask_all_touched)
sm.visualize(ax=ax2, addcbar=False, title='All touched')
# plot!
plt.tight_layout()
plt.show()
PK +KXK. . plot_subgrid_mask.py# -*- coding: utf-8 -*-
"""
==========================
Shape to raster -- subgrid
==========================
Compute a percentage raster from a shapefile.
This example extends the land/sea shape to raster example with a subgrid land
cover mask.
"""
import salem
import matplotlib.pyplot as plt
# make a local grid from which we will compute the mask
# we make it coarse so that we see the grid points
grid = salem.Grid(proj=salem.wgs84, x0y0=(-18, 3), nxny=(25, 15), dxdy=(1, 1))
# make a high-res subgrid grid to compute the mask
hr_grid = grid.regrid(factor=10)
# read the ocean shapefile (data from http://www.naturalearthdata.com)
oceans = salem.read_shapefile(salem.get_demo_file('ne_50m_ocean.shp'),
cached=True)
# read the lake shapefile (data from http://www.naturalearthdata.com)
lakes = salem.read_shapefile(salem.get_demo_file('ne_50m_lakes.shp'),
cached=True)
# Now compute the hr mask
hr_mask = hr_grid.region_of_interest(shape=oceans)
hr_mask = hr_grid.region_of_interest(shape=lakes, roi=hr_mask)
# And reduce it to the original grid
perc_mask = salem.reduce(hr_mask, factor=10)
# Make a map to check our results
sm = salem.Map(grid, countries=False)
sm.set_plot_params(cmap='Blues', vmax=1)
# prepare the figure
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3))
# plot 1
sm.set_data(hr_mask, hr_grid)
sm.visualize(ax=ax1, addcbar=False, title='High-res grid')
# plot 2
sm.set_data(perc_mask)
sm.visualize(ax=ax2, addcbar=True, title='Percentage mask')
# plot!
plt.tight_layout()
plt.show()
PK (KXNq plot_wind_overlay.py# -*- coding: utf-8 -*-
"""
=============
Plot overlays
=============
Add contours and wind arrows to a salem plot
"""
import salem
import numpy as np
import matplotlib.pyplot as plt
# get the data at the latest time step
ds = salem.open_wrf_dataset(salem.get_demo_file('wrfout_d01.nc')).isel(time=-1)
# get the wind data at 10000 m a.s.l.
u = ds.salem.wrf_zlevel('U', 10000.)
v = ds.salem.wrf_zlevel('V', 10000.)
ws = ds.salem.wrf_zlevel('WS', 10000.)
# get the axes ready
f, ax = plt.subplots()
# plot the salem map background, make countries in grey
smap = ds.salem.get_map(countries=False)
smap.set_shapefile(countries=True, color='grey')
smap.plot(ax=ax)
# transform the coordinates to the map reference system and contour the data
xx, yy = smap.grid.transform(ws.west_east.values, ws.south_north.values,
crs=ws.salem.grid.proj)
cs = ax.contour(xx, yy, ws, cmap='viridis', levels=np.arange(0, 81, 10),
linewidths=2)
# Quiver only every 7th grid point
u = u[4::7, 4::7]
v = v[4::7, 4::7]
# transform their coordinates to the map reference system and plot the arrows
xx, yy = smap.grid.transform(u.west_east.values, u.south_north.values,
crs=u.salem.grid.proj)
xx, yy = np.meshgrid(xx, yy)
qu = ax.quiver(xx, yy, u.values, v.values)
qk = plt.quiverkey(qu, 0.7, 0.95, 50, '50 m s$^{-1}$',
labelpos='E', coordinates='figure')
# done!
plt.show()
PK KXrjY Y plot_change_style.py# -*- coding: utf-8 -*-
"""
====================
Customize Salem maps
====================
How to change the look of a Map?
"""
import salem
import matplotlib.pyplot as plt
# get the map from a WRF file
ds = salem.open_wrf_dataset(salem.get_demo_file('wrfout_d01.nc'))
smap = ds.salem.get_map(countries=False)
# Change the country borders
smap.set_shapefile(countries=True, color='C3', linewidths=2)
# Add oceans and lakes
smap.set_shapefile(oceans=True)
smap.set_shapefile(rivers=True)
smap.set_shapefile(lakes=True, facecolor='blue', edgecolor='blue')
# Change the lon-lat countour setting
smap.set_lonlat_contours(add_ytick_labels=False, interval=5, linewidths=1.5,
linestyles='-', colors='C1')
# Add a scalebar (experimental)
smap.set_scale_bar(location=(0.87, 0.04), add_bbox=True)
# done!
smap.visualize()
plt.show()
PK KX{' ' plot_topo_shading.py# -*- coding: utf-8 -*-
"""
===================
Topographic shading
===================
Add topographic shading to a plot
"""
from salem import mercator_grid, Map, get_demo_file
import matplotlib.pyplot as plt
# prepare the figure
f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(8, 7))
# map extent
grid = mercator_grid(center_ll=(10.76, 46.79), extent=(18000, 14000))
sm = Map(grid, countries=False)
sm.set_lonlat_contours(interval=0)
sm.set_scale_bar()
# add topography
fpath = get_demo_file('hef_srtm.tif')
sm.set_topography(fpath)
sm.visualize(ax=ax1, addcbar=False, title='relief_factor=0.7 (default)')
# stronger shading
sm.set_topography(fpath, relief_factor=1.4)
sm.visualize(ax=ax2, addcbar=False, title='relief_factor=1.4')
# add color shading
z = sm.set_topography(fpath)
sm.set_data(z)
sm.set_cmap('topo')
sm.visualize(ax=ax3, title='Color with shading', addcbar=False)
# color without topo shading
sm.set_topography()
sm.visualize(ax=ax4, title='Color without shading', addcbar=False)
# make it nice
plt.tight_layout()
plt.show()
PK 7KX} plot_topo_comparison.py# -*- coding: utf-8 -*-
"""
========================================
Compare datasets of different resolution
========================================
An example of use for salem's lookup_transform
.. currentmodule:: salem
In this example, we compare a model topography defined at 1km resolution with
the topography from the SRTM v4.1 dataset (resolution 3 minutes of arc,
so ~ 90m). For this we use the :py:meth:`DataArrayAccessor.lookup_transform`
method.
From the plot below, we see that the model topography is smoother than the
aggregated SRTM (this is a good thing, as atmospheric models do not like
steep gradients too much). The highest peaks or lowest valley aren't resolved
by the 1km topography.
"""
import numpy as np
from salem import get_demo_file, open_xr_dataset
import matplotlib.pyplot as plt
# get the topography data
srtm = open_xr_dataset(get_demo_file('riosan_srtm_hgt.nc')).srtm
wrf = open_xr_dataset(get_demo_file('riosan_wrf_hgt.nc')).HGT
# transform the high-res topography onto the coarse grid
# we ask for the lookup table to speed up the second transform
srtm_on_wrf, lut = wrf.salem.lookup_transform(srtm, return_lut=True)
srtm_on_wrf_std = wrf.salem.lookup_transform(srtm, method=np.std, lut=lut)
# for fun we compute the max and min for each grid point
srtm_on_wrf_min = wrf.salem.lookup_transform(srtm, method=np.min, lut=lut)
srtm_on_wrf_max = wrf.salem.lookup_transform(srtm, method=np.max, lut=lut)
# then compute the max absolute difference to wrf
absdif = np.abs(np.dstack([srtm_on_wrf_min - wrf, srtm_on_wrf_max - wrf]))
maxabsdif = np.max(absdif, axis=2)
# Get the map defined by the WRF grid
sm = wrf.salem.get_map(cmap='topo')
# remove the lon-lat ticks for clarity
sm.set_lonlat_contours(interval=0)
# prepare the figure and plot
f, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2, 3, figsize=(11, 7))
# absolute values
sm.set_data(srtm)
sm.visualize(ax=ax1, title='SRTM 90m')
sm.set_data(wrf)
sm.visualize(ax=ax2, title='WRF 1km')
sm.set_data(srtm_on_wrf)
sm.visualize(ax=ax3, title='SRTM 1km')
# comparisons
sm.set_data(srtm_on_wrf_std)
sm.set_plot_params(vmin=0, cmap='Purples')
sm.visualize(ax=ax4, title='Std. Dev of SRTM')
sm.set_data(wrf - srtm_on_wrf)
sm.set_plot_params(levels=np.linspace(-250, 250, 11), cmap='RdBu')
sm.visualize(ax=ax5, title='Diff WRF - SRTM')
sm.set_data(maxabsdif)
sm.set_plot_params(vmin=0, cmap='OrRd')
sm.visualize(ax=ax6, title='Max. abs. diff.')
# make it nice
plt.tight_layout()
plt.show()
PK KXXT plot_natural_earth.py# -*- coding: utf-8 -*-
"""
==============================
Add a Natural Earth background
==============================
An alternative to Google Static Maps
"""
import salem
import matplotlib.pyplot as plt
# get the map from a predefined grid
grid = salem.mercator_grid(transverse=False, center_ll=(16., 0.),
extent=(8e6, 9e6))
smap = salem.Map(grid)
# Add the background (other resolutions include: 'mr', 'hr')
smap.set_rgb(natural_earth='lr')
# done!
smap.visualize()
plt.show()
PK 8KX) plot_hydrosheds.pyPK !KX # plot_polarstereo.pyPK .KXoܝ6 6 plot_googlestatic.pyPK /KXH H D plot_landsea_mask.pyPK +KXK. . plot_subgrid_mask.pyPK (KXNq &