How to use the salem.get_demo_file function in salem

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

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

# -*- 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)

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

"""

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)

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')

# 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()