Salem¶

Subsetting and selecting data¶

Let’s open a WRF model output file:

In [1]: import salem

In [2]: from salem.utils import get_demo_file

In [3]: ds = salem.open_xr_dataset(get_demo_file('wrfout_d01.nc'))

Let’s take a time slice of the variable T2 for a start:

In [4]: t2 = ds.T2.isel(Time=2)

In [5]: t2.salem.quick_map()
Out[5]: <salem.graphics.Map at 0x7efe27aa3410>

Although we are on a Lambert Conformal projection, it’s possible to subset the file using longitudes and latitudes:

In [6]: t2_sub = t2.salem.subset(corners=((77., 20.), (97., 35.)), crs=salem.wgs84)

In [7]: t2_sub.salem.quick_map()
Out[7]: <salem.graphics.Map at 0x7efe2715f710>

It’s also possible to use geometries or shapefiles to subset your data:

In [8]: shdf = salem.read_shapefile(get_demo_file('world_borders.shp'))

In [9]: shdf = shdf.loc[shdf['CNTRY_NAME'].isin(['Nepal', 'Bhutan'])]  # GeoPandas' GeoDataFrame

In [10]: t2_sub = t2_sub.salem.subset(shape=shdf, margin=2)  # add 2 grid points

In [11]: t2_sub.salem.quick_map()
Out[11]: <salem.graphics.Map at 0x7efe286415d0>

Based on the same principle, one can mask out the useless grid points:

In [12]: t2_roi = t2_sub.salem.roi(shape=shdf)

In [13]: t2_roi.salem.quick_map()
Out[13]: <salem.graphics.Map at 0x7efe26f990d0>

Plotting¶

Maps can be pimped with topographical shading, points of interest, and more:

In [14]: smap = t2_roi.salem.get_map(data=t2_roi-273.15, cmap='RdYlBu_r', vmin=-14, vmax=18)

In [15]: _ = smap.set_topography(get_demo_file('himalaya.tif'))

In [16]: smap.set_shapefile(shape=shdf, color='grey', linewidth=3)

In [17]: smap.set_points(91.1, 29.6)

In [18]: smap.set_text(91.2, 29.7, 'Lhasa', fontsize=17)

In [19]: smap.visualize()

Maps are persistent, which is useful when you have many plots to do. Plotting further data on them is possible, as long as the geolocalisation information is shipped with the data (in that case, the DataArray’s attributes are lost in the conversion from Kelvins to degrees Celsius so we have to set it explicitly):

In [20]: smap.set_data(ds.T2.isel(Time=1)-273.15, crs=ds.salem.grid)

In [21]: smap.visualize(title='2m temp - large domain', cbar_title='C')

Reprojecting data¶

Salem can also transform data from one grid to another:

In [22]: dse = salem.open_xr_dataset(get_demo_file('era_interim_tibet.nc'))

In [23]: t2_era_reproj = ds.salem.transform(dse.t2m)

In [24]: assert t2_era_reproj.salem.grid == ds.salem.grid

In [25]: t2_era_reproj.isel(time=0).salem.quick_map()
Out[25]: <salem.graphics.Map at 0x7efe280ea850>

In [26]: t2_era_reproj = ds.salem.transform(dse.t2m, interp='spline')

In [27]: t2_era_reproj.isel(time=0).salem.quick_map()
Out[27]: <salem.graphics.Map at 0x7efe27d8da50>

Is your library mature for production code?¶

No. The API is not always as clever as I wish it would, and it will probably change in the future. It is quite well tested though, at least for the cases I encounter in my daily work.

What others tools should I know about?¶

If you want to plot on maps, cartopy is probably one of the best tools you could pick. For reprojection workflows on large or numerous gridded files you probably want to use rasterio.

The python atmopsheric sciences community is a bit spread between iris and xarray for N-Dimensional data handling (I picked the later for it’s strong interaction with pandas). Several great libraries are available to meteorologists and climatologists, for example MetPy, windspharm, xgcm, and all the ones I forgot to mention.

But then, why developing Salem?¶

As an atmospheric scientist, I hate to have to take care about projections and maps. Salem was created to hide all these concerns. By the time I started, it seemed a good idea to provide map transformation tools without depending on GDAL (thanks to conda-forge GDAL is now much easier to install). It is still possible to do reprojection work in Salem using scipy and pyproj alone. It is al done in python and on memory, so don’t expect miracles on that side.

Furthermore, I use the atmospheric model WRF quite often in my work. Its output files are absolutely not compliant with the CF conventions. To my knowledge, there is no tool to plot and manipulate WRF data with Python, and Salem will be further developed with this model in mind.

Why aren’t you using Xarray/Cartopy for your maps?¶

This is actually a good question, and I think that there are no obstacle for Salem to use it’s geolocation information to plot gridded data on cartopy’s maps. I never really got to understand how cartopy really works, but I’m sure that salem’s Grid object would be able to provide the right info to cartopy with a little bit of coding. Contributions welcome!

On the other hand I kind of like how Salem’s Map output looks, and I find it nice and easy to use. But that, of course, might not be your opinion. Other advantages of Salem’s Maps are their persistence (usefull if you want to plot data on the same map several times), their internal optimisations based on disk cache, and the ability to add lon-lat orientation grids to virtually any map projection.

What’s this ”.salem_cache” directory in my home folder?¶

At the first import, Salem will create a hidden directory called .salem_cache in your home folder. It will be used to download Salem’s demo files and standard shapefiles. This directory is also used by joblib to store the result of slow operations such as reading and transforming shapefiles, or downloading google maps from the internet. The cache should not become too large, but if it does: simply delete it.

v0.1.0 (22 October 2016)¶

Big refactoring (PR15), partly backwards incompatible (mostly renaming). Improved xarray accessors, WRF tools, merged Cleo into the codebase, added a first draft of documentation.

v0.0.9 (22 October 2016)¶

Initial release.

Required dependencies¶

• Python 2.7 or 3+ (Py3 recommended)
• numpy (of course)
• scipy: for its interpolation tools, among other things
• pyproj: for map projections
• netCDF4: to read most geoscientific files
• joblib: for it’s Memory class
• six: for Py2 compatibility

Optional dependencies¶

Instructions¶

The very best (unique?) way to install Salem without too much hassle is to install its dependencies with conda and conda-forge:

conda config --add channels conda-forge
conda install <package-name>

Currently, Salem can only be installed via pip:

pip install salem

If you want to install the latest master:

pip install git+https://github.com/fmaussion/salem.git

Grids¶

A Grid is defined by a projection, a reference point in this projection, a grid spacing and a number of grid points:

In [1]: import numpy as np

In [2]: import salem

In [3]: from salem import wgs84

In [4]: grid = salem.Grid(nxny=(3, 2), dxdy=(1, 1), ll_corner=(0.5, 0.5), proj=wgs84)

In [5]: x, y = grid.xy_coordinates

In [6]: x
Out[6]: 
array([[ 0.5,  1.5,  2.5],
       [ 0.5,  1.5,  2.5]])

In [7]: y
Out[7]: 
array([[ 0.5,  0.5,  0.5],
       [ 1.5,  1.5,  1.5]])

The default in Salem is to define the grids according to the pixels center point:

In [8]: smap = salem.Map(grid)

In [9]: smap.set_data(np.arange(6).reshape((2, 3)))

In [10]: lon, lat = grid.ll_coordinates

In [11]: smap.set_points(lon.flatten(), lat.flatten())

In [12]: smap.visualize(addcbar=False)

But with the pixel_ref keyword you can use another convention. For Salem, the two conventions are identical:

In [13]: grid_c = salem.Grid(nxny=(3, 2), dxdy=(1, 1), ll_corner=(0, 0), proj=wgs84, pixel_ref='corner')

In [14]: assert grid_c == grid

While it’s good to know how grids work, most of the time grids should be inferred directly from the data files (if not, we have to implement it!):

In [15]: ds = salem.open_xr_dataset(salem.get_demo_file('himalaya.tif'))

In [16]: grid = ds.salem.grid

In [17]: grid.proj.srs
Out[17]: '+units=m +init=epsg:4326 '

In [18]: grid.extent
Out[18]: 
[78.312513473116013,
 93.937514288103202,
 24.237496572347762,
 31.545830286884001]

Grids come with several convenience functions, for example for transforming a point onto the grid coordinates:

In [19]: grid.transform(85, 27, crs=salem.wgs84)
Out[19]: (801.99834136842378, 544.9996059728029)

Or for reprojecting structured data (the xarray accessors transform() method calls map_gridded_data() internally):

In [20]: dse = salem.open_xr_dataset(salem.get_demo_file('era_interim_tibet.nc'))

In [21]: t2_era_reproj = ds.salem.transform(dse.t2m.isel(time=0))

In [22]: t2_era_reproj.salem.quick_map()
Out[22]: <salem.graphics.Map at 0x7efe27d8e890>

Color handling with DataLevels¶

DataLevels is the base class for handling colors. It is there to ensure that there will never be a mismatch between data and assigned colors.

In [1]: from salem import DataLevels

In [2]: a = [-1., 0., 1.1, 1.9, 9.]

In [3]: dl = DataLevels(a)

In [4]: dl.visualize(orientation='horizontal', add_values=True)

Discrete levels¶

In [5]: dl.set_plot_params(nlevels=11)

In [6]: dl.visualize(orientation='horizontal', add_values=True)

vmin, vmax¶

In [7]: dl.set_plot_params(nlevels=9, vmax=3)

In [8]: dl.visualize(orientation='horizontal', add_values=True)

Out-of-bounds data¶

In [9]: dl.set_plot_params(levels=[0, 1, 2, 3])

In [10]: dl.visualize(orientation='horizontal', add_values=True)

Note that if the bounds are not exceeded, the colorbar extensions are gone:

In [11]: dl.set_data([0., 0.2, 0.4, 0.6, 0.8])

In [12]: dl.visualize(orientation='horizontal', add_values=True)

This might be undesirable, so you can set a keyword to force out-of-bounds levels:

In [13]: dl.set_plot_params(levels=[0, 1, 2, 3], extend='both')

In [14]: dl.visualize(orientation='horizontal', add_values=True)

Using DataLevels with matplotlib¶

Here with the example of a scatterplot:

In [15]: x, y = np.random.randn(1000), np.random.randn(1000)

In [16]: z = x**2 + y**2

In [17]: dl = DataLevels(z, cmap='RdYlBu_r', levels=np.arange(6))

In [18]: fig, ax = plt.subplots(1, figsize=(6, 4))

In [19]: ax.scatter(x, y, color=dl.to_rgb(), s=64);

In [20]: cbar = dl.append_colorbar(ax, "right")  # DataLevel draws the colorbar

In [21]: plt.show()

Maps¶

Map is a sublass of DataLevels, but adds the georeferencing aspects to the plot. A Map is initalised with a Grid:

In [22]: from salem import mercator_grid, Map, get_demo_file, open_xr_dataset

In [23]: grid = mercator_grid(center_ll=(10.76, 46.79), extent=(9e5, 4e5))

In [24]: emap = Map(grid)

In [25]: emap.visualize(addcbar=False)

The map image has it’s own pixel resolution (set with the keywords nx or ny), and the cartographic information is simply overlayed on it. When asked to plot data, the map will automatically transform it to the map projection:

In [26]: ds = open_xr_dataset(get_demo_file('histalp_avg_1961-1990.nc'))

In [27]: emap.set_data(ds.prcp)

In [28]: emap.visualize()

Add topographical shading to a map¶

You can add topographical shading to a map with DEM files:

In [29]: grid = mercator_grid(center_ll=(10.76, 46.79), extent=(18000, 14000))

In [30]: smap = Map(grid, countries=False)

In [31]: smap.set_topography(get_demo_file('hef_srtm.tif'));

In [32]: smap.visualize(addcbar=False, title='Topographical shading')

Note that you can also use the topography data to make a colourful plot:

In [33]: z = smap.set_topography(get_demo_file('hef_srtm.tif'))

In [34]: smap.set_data(z)

In [35]: smap.set_cmap('topo')

In [36]: smap.visualize(title='Topography', cbar_title='m a.s.l.')

Diagnostic variables¶

Salem adds a layer between xarray and the underlying NetCDF file. This layer computes new variables “on the fly” or, in the case of staggered variables, “unstaggers” them:

In [8]: ds.U
Out[8]: 
<xarray.DataArray 'U' (time: 3, bottom_top: 27, south_north: 150, west_east: 150)>
[1822500 values with dtype=float32]
Coordinates:
    lon          (time, south_north, west_east) float32 70.7231 70.9749 ...
    lat          (time, south_north, west_east) float32 7.78907 7.82944 ...
  * time         (time) datetime64[ns] 2008-10-26T12:00:00 ...
  * south_north  (south_north) float64 -2.235e+06 -2.205e+06 -2.175e+06 ...
  * west_east    (west_east) float64 -2.235e+06 -2.205e+06 -2.175e+06 ...
  * bottom_top   (bottom_top) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ...
Attributes:
    FieldType: 104
    MemoryOrder: XYZ
    description: x-wind component
    units: m s-1
    stagger: X
    pyproj_srs: +units=m +proj=lcc +lat_1=29.0400009155 +lat_2=29.0400009155 +lat_0=29.0399971008 +lon_0=89.8000030518 +x_0=0 +y_0=0 +a=6370000 +b=6370000

This computation is done only on demand (just like a normal NetCDF variable), this new layer is therefore relatively cheap.

In addition to unstaggering, Salem adds a number of “diagnostic” variables to the dataset. Some are convenient (like T2C, temperature in Celsius instead of Kelvins), but others are more complex (e.g. SLP for sea-level pressure, or PRCP which computes step-wize total precipitation out of the accumulated fields). For a list of diagnostic variables (and TODOs!), refer to GH18.

In [9]: ds.PRCP.isel(time=-1).salem.quick_map(cmap='Purples', vmax=5)
Out[9]: <salem.graphics.Map at 0x7efe26c3afd0>

Geogrid simulator¶

Salem provides a small tool which comes handy when defining new WRF domains. It parses the WPS configuration file (namelist.wps) and generates the grids and maps corresponding to each domain.

geogrid_simulator() will search for the &geogrid section of the file and parse it:

In [10]: fpath = get_demo_file('namelist_mercator.wps')

In [11]: with open(fpath, 'r') as f:  # this is just to show the file
   ....:     print(f.read())
   ....: 
&geogrid
 parent_id         =   1, 1, 2, 3,
 parent_grid_ratio =   1, 3, 5, 5,
 i_parent_start    =   1, 28, 54, 40,
 j_parent_start    =   1, 24, 97, 50,
 e_we              =  211, 178, 101, 136,
 e_sn              =  131, 220, 141, 136,
 dx = 30000,
 dy = 30000,
 map_proj = 'mercator',
 ref_lat   =  -8.0,
 ref_lon   =  -62.0,
 truelat1  =  -8.0,
/



In [12]: from salem import geogrid_simulator

In [13]: g, maps = geogrid_simulator(fpath)

In [14]: maps[0].set_rgb(natural_earth='lr')  # add a background image

In [15]: maps[0].visualize(title='Domains 1 to 4')

Grid¶

Georeferencing utils¶

Graphics¶

Input/output¶

Xarray¶

Old-style datasets¶

Old-style Datasets (prior to xarray), kept for backwards compatibility reasons and because they are quite lightweight. They might be replaced by xarray’s datasets one day.

Other¶