10: Using Rasterio and Numpy to examine ice loss on Mt. Rainier

This exercise focuses on working with GIS raster data. It is based on a GeoHackWeek tutorial. We will use the rasterio and rasterstats libraries, along with numpy and matplotlib to look changes in the surface elevations of mapped glaciers since 1970.

Key takeaways:

  • raster formats such as GeoTiff map array data on a regular grid using a simple transform (offset, rotation and pixel size).

  • working with raster data in python is inherently a little difficult, because the raster data model is designed to work efficiently with large datasets on disk (without loading whole thing into memory). Other available tools, such as the gdal commandline tools or ArcPy macros, tend to stay in this paradigm by taking files as input to functions that write files as output.

  • conversely, when we work with arrays in-memory with numpy, we typically think it terms of rows and columns

The main idea is to

  1. use rasterio (or another tool) to handle the files and any transformations needed to align the data. As we’ll see below, most of these operations aren’t especially intuitive or easy to remember. Thankfully, there are recipes available online that are easy to copy for performing common tasks. The rasterio manual is a great place to start.

  2. once the data are in numpy, use numpy and other tools in the scientific python stack (e.g. scipy) to do the computations

Datasets:

  • 3 Digital Elevation Models (DEMs) of Mt. Rainier elevations in 1970, 2008 and 2015; each on a different grid.

  • shapefile of glacier extents

Operations:

  • read and write rasters to/from numpy arrays

  • plotting raster data

  • warping (resampling) raster data to different grids

  • reprojecting raster data to different coordinate reference systems

  • getting raster data at discrete points

  • zonal statistics

  • mapping vector features onto rasters (“rasterizing”)

References:

[1]:

import numpy as np
import geopandas as gp
import rasterio
from rasterio.plot import show
from rasterio import features
import matplotlib.pyplot as plt
from matplotlib.colors import LightSource
import pathlib as pl

# for elevation hillshade if we want to use it
ls = LightSource(azdeg=315, altdeg=45)

Define out input data

We’ll be working with three rasters and a single shapefile for this exercise

[2]:

data_path = pl.Path("data/rasterio")

dem_file_1970 = data_path / "19700901_ned1_2003_adj_warp.tif"
dem_file_2008 = data_path / "20080901_rainierlidar_30m-adj.tif"
dem_file_2015 = data_path / "20150818_rainier_summer-tile-30.tif"

glaciers_shapefile = data_path / "rgi60_glacierpoly_rainier.shp"

input_rasters = {
    1970: dem_file_1970,
    2008: dem_file_2008,
    2015: dem_file_2015
}

Quickly view a raster with rasterio

[3]:

with rasterio.open(input_rasters[1970]) as src:
    show(src)
../../_images/notebooks_part0_python_intro_10_Rasterio_5_0.png

Get raster values at point locations

if we have a list of (x, y) tuples, we can unpack it into two lists of x and y values: zip(*points)

zip(*points) returns a tuple of the two lists

First let’s define visualize our point locations on our raster

[4]:

points = [
    (594508.70, 5189550.08), # summit, in utm zone 10N
    (596702, 5187675) # camp muir
]
x, y = zip(*points)
fig, ax = plt.subplots(figsize=(8, 8))
with rasterio.open(input_rasters[2015]) as src:
     ax = show(src, ax=ax, zorder=1, cmap="plasma", vmin=500, vmax=4500)

ax.scatter(x, y, c="k", s=50, zorder=5)
ax.set_facecolor("lightgrey")
plt.colorbar(ax.images[0], shrink=0.7);
../../_images/notebooks_part0_python_intro_10_Rasterio_8_0.png 
[5]:

xpts, ypts = list(zip(*points))
# points

Now let’s sample elevations at these points

[6]:

with rasterio.open(input_rasters[2015]) as src:
    meta = src.meta
    data = src.read(1)
    i, j = [], []
    for ix in range(len(xpts)):
        ti, tj = src.index(xpts[ix], ypts[ix])
        i.append(ti)
        j.append(tj)
print(data[i, j])

[4361.6406 3055.5547]

Viewing the 3 rasters side by side as subplots

The data are all on different grids with varying extents and projections. Let’s plot them up to visualize this.

[7]:

fig, axes = plt.subplots(1, 3, figsize=(21, 7))

meta = {}
for i, (year, f) in enumerate(input_rasters.items()):
    with rasterio.open(f) as src:

        # get the metadata
        meta[year] = src.meta

        # plot the data
        if i == 0:
            data = src.read(1)
            vmin, vmax = data.min(), data.max()
        ax=axes.flat[i]
        show(src, ax=ax, vmin=vmin, vmax=vmax)
        ax.set_title(year)
        ax.set_facecolor("lightgrey")
        plt.colorbar(ax.images[0], ax=ax,
                     orientation='horizontal', label='elevation, m')
../../_images/notebooks_part0_python_intro_10_Rasterio_13_0.png

We see what these are by fetching the metadata dictionary (src.meta attribute) for each one. Each grid is described in the metadata:

[8]:

meta

[8]:

{1970: {'driver': 'GTiff',
  'dtype': 'float32',
  'nodata': -9999.0,
  'width': 884,
  'height': 824,
  'count': 1,
  'crs': CRS.from_wkt('PROJCS["NAD83 / UTM zone 11N",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","26911"]]'),
  'transform': Affine(30.0, 0.0, 125349.01533479267,
         0.0, -30.0, 5213641.54503958)},
 2008: {'driver': 'GTiff',
  'dtype': 'float32',
  'nodata': -9999.0,
  'width': 1249,
  'height': 1093,
  'count': 1,
  'crs': CRS.from_wkt('PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]'),
  'transform': Affine(30.0, 0.0, 581379.9999992945,
         0.0, -30.0, 5206284.00012134)},
 2015: {'driver': 'GTiff',
  'dtype': 'float32',
  'nodata': 0.0,
  'width': 827,
  'height': 761,
  'count': 1,
  'crs': CRS.from_wkt('PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]'),
  'transform': Affine(30.0, 0.0, 583464.253406,
         0.0, -30.0, 5202549.0989)}}

Warp the 1970 and 2008 rasters onto the 2015 grid

The code below was copied from the rasterio manual. The result is a new set of files aligned to the destination grid.

[9]:

from rasterio.enums import Resampling
from rasterio.vrt import WarpedVRT
from rasterio import shutil as rio_shutil


vrt_options = {
    'resampling': Resampling.cubic,
    'crs': meta[2015]['crs'],
    'transform': meta[2015]['transform'],
    'height': meta[2015]['height'],
    'width': meta[2015]['width'],
}

aligned_rasters = {}
for year, path in input_rasters.items():

    with rasterio.open(path) as src:

        with WarpedVRT(src, **vrt_options) as vrt:

            # At this point 'vrt' is a full dataset with dimensions,
            # CRS, and spatial extent matching 'vrt_options'.

            # Read all data into memory.
            data = vrt.read()

            # Process the dataset in chunks.  Likely not very efficient.
            for _, window in vrt.block_windows():
                data = vrt.read(window=window)

            # Dump the aligned data into a new file.  A VRT representing
            # this transformation can also be produced by switching
            # to the VRT driver.
            name = path.name
            outfile = data_path / f'aligned-{name}'
            rio_shutil.copy(vrt, outfile, driver='GTiff')
            print('wrote {}'.format(outfile))
            aligned_rasters[year] = outfile

wrote data/rasterio/aligned-19700901_ned1_2003_adj_warp.tif
wrote data/rasterio/aligned-20080901_rainierlidar_30m-adj.tif
wrote data/rasterio/aligned-20150818_rainier_summer-tile-30.tif

Now let’s make sure all three are on the same grid

[10]:

meta = {}
for year, f in aligned_rasters.items():
    with rasterio.open(f) as src:
        meta[year] = src.meta
meta

[10]:

{1970: {'driver': 'GTiff',
  'dtype': 'float32',
  'nodata': -9999.0,
  'width': 827,
  'height': 761,
  'count': 1,
  'crs': CRS.from_wkt('PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]'),
  'transform': Affine(30.0, 0.0, 583464.253406,
         0.0, -30.0, 5202549.0989)},
 2008: {'driver': 'GTiff',
  'dtype': 'float32',
  'nodata': -9999.0,
  'width': 827,
  'height': 761,
  'count': 1,
  'crs': CRS.from_wkt('PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]'),
  'transform': Affine(30.0, 0.0, 583464.253406,
         0.0, -30.0, 5202549.0989)},
 2015: {'driver': 'GTiff',
  'dtype': 'float32',
  'nodata': 0.0,
  'width': 827,
  'height': 761,
  'count': 1,
  'crs': CRS.from_wkt('PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]'),
  'transform': Affine(30.0, 0.0, 583464.253406,
         0.0, -30.0, 5202549.0989)}}

Success!

Read the aligned grids into a new set of masked numpy arrays

rasterio can be used to construct masked numpy arrays based on the 'nodata' value in the raster metadata

[11]:

aligned_data = {}
for year, f in aligned_rasters.items():
    with rasterio.open(f) as src:
        data = src.read(1, masked=True)
        aligned_data[year] = data

[12]:

aligned_data[2015]

[12]:

masked_array(
  data=[[--, --, --, ..., --, --, --],
        [--, --, --, ..., --, --, --],
        [--, --, --, ..., --, --, --],
        ...,
        [--, 1470.8792724609375, 1477.231201171875, ..., --, --, --],
        [1451.615966796875, 1456.029296875, 1460.4559326171875, ..., --,
         --, --],
        [1429.700439453125, --, --, ..., --, --, --]],
  mask=[[ True,  True,  True, ...,  True,  True,  True],
        [ True,  True,  True, ...,  True,  True,  True],
        [ True,  True,  True, ...,  True,  True,  True],
        ...,
        [ True, False, False, ...,  True,  True,  True],
        [False, False, False, ...,  True,  True,  True],
        [False,  True,  True, ...,  True,  True,  True]],
  fill_value=0.0,
  dtype=float32)

[13]:

fig, ax = plt.subplots(figsize=(8, 8))
im = ax.imshow(aligned_data[2015])
ax.set_facecolor("lightgrey")
plt.colorbar(im, shrink=0.7);
../../_images/notebooks_part0_python_intro_10_Rasterio_24_0.png

Assign information on shared raster grid to variables for later use

[14]:

epsg_32610 = meta[2015]['crs']
transform = meta[2015]['transform']
nrow, ncol = aligned_data[2015].shape
cellsize = transform[0]
extent = src.bounds[0], src.bounds[2], src.bounds[1], src.bounds[3]  # extent for plotting (xmin, xmax, ymin, ymax)

Get the summit elevation

Ss noted in the original GeoHackWeek example, the elevation values in these rasters are relative to the WGS84 datum, not a geoid model that approximate sea level. The geoid offset at Rainier is approximately 18.79 m, meaning that we have to add we have to add this to the values to get the correct summit elevation of 14,411 ft.

[15]:

summit_ij = np.unravel_index(aligned_data[1970].argmax(), aligned_data[1970].shape)
summit_ij

[15]:

(np.int64(433), np.int64(367))

The corrected summit elevations being a little under 14,411 is a little unsatisfying, but one of the costs of downsampling the DEMs to a 30 m resolution to conserve space. At any rate, since we are interested in differences, the geoid offset doesn’t matter for the remainder of our analysis.

[16]:

elevation_offset = 18.79
(aligned_data[1970][summit_ij] + elevation_offset) * 3.2084  # in feet

[16]:

np.float32(14085.793)

Compute the elevation differences

[17]:

elevation_diffs = {}
elevation_diffs['1970-2008'] = aligned_data[2008] - aligned_data[1970]
elevation_diffs['2008-2015'] = aligned_data[2015] - aligned_data[2008]
elevation_diffs['1970-2015'] = aligned_data[2015] - aligned_data[1970]

fig, axes = plt.subplots(1, 3, figsize=(21, 7))

for i, (dates, diffs) in enumerate(elevation_diffs.items()):

    ax=axes.flat[i]
    im = ax.imshow(diffs, cmap='RdBu', vmin=-30, vmax=30)
    ax.set_facecolor("lightgrey")
    ax.set_title(dates)

fig.colorbar(im, ax=axes.ravel(), label='elevation change, m', shrink=0.6, extend='both');
../../_images/notebooks_part0_python_intro_10_Rasterio_32_0.png

Average rates of change

[18]:

elevation_rates = {}
for dates, diffs in elevation_diffs.items():
    yr0, yr1 = map(int, dates.split('-'))
    annual_rate = diffs/(yr1-yr0)
    elevation_rates[dates] = annual_rate

fig, axes = plt.subplots(1, 3, figsize=(21, 7))

for i, (dates, rates) in enumerate(elevation_rates.items()):

    ax=axes.flat[i]
    im = ax.imshow(rates, cmap='RdBu', vmin=-2, vmax=2)
    ax.set_facecolor("lightgrey")
    ax.set_title(dates)

fig.colorbar(im, ax=axes.ravel(), label='average rate of change, m/yr', shrink=0.6, extend='both');
../../_images/notebooks_part0_python_intro_10_Rasterio_34_0.png

Notice the huge changes around the margins of the mountain. These are due to the 1970 and 2008 DEMs representing “bare ground”, while the 2015 DEM includes trees.

Map the shapefile of glacier extents onto the same grid as the rasters

  • first reproject the glacier features to the same crs as the rasters

  • then make a column of unique identifiers (integers) for relating the features to raster pixels

[19]:

glaciers = gp.read_file(glaciers_shapefile)
glaciers.head()

[19]:
RGIId GLIMSId BgnDate EndDate CenLon CenLat O1Region O2Region Area Zmin ... Aspect Lmax Status Connect Form TermType Surging Linkages Name geometry
0 RGI60-02.14393 G238292E46809N 19709999 19949999 -121.70769 46.80930 2 4 0.126 2001 ... 148 397 0 0 0 0 0 9 Williwakas Glacier WA POLYGON ((-121.70642 46.80792, -121.7065 46.80...
1 RGI60-02.14394 G238275E46811N 19709999 19949999 -121.72504 46.81055 2 4 0.010 2207 ... 222 -9 0 0 0 0 0 9 WA POLYGON ((-121.72552 46.80985, -121.72561 46.8...
2 RGI60-02.14395 G238235E46811N 19709999 19949999 -121.76535 46.81088 2 4 0.010 2103 ... 165 -9 0 0 0 0 0 9 WA POLYGON ((-121.76498 46.81006, -121.76511 46.8...
3 RGI60-02.14396 G238243E46810N 19709999 19949999 -121.75702 46.80963 2 4 0.015 2186 ... 153 182 0 0 0 0 0 9 WA POLYGON ((-121.75659 46.80885, -121.75663 46.8...
4 RGI60-02.14397 G238302E46809N 19709999 19949999 -121.69756 46.80850 2 4 0.012 1990 ... 125 160 0 0 0 0 0 9 WA POLYGON ((-121.69777 46.80903, -121.6976 46.80...

5 rows × 23 columns

[20]:

glaciers.to_crs(epsg_32610, inplace=True)
glaciers['id'] = np.arange(1, len(glaciers) + 1)

Use rasterio to rasterize the glacier extents

  • input is a collection of the (feature geometry, unique value) tuples and an Affine transform. Further reading on affine transforms can be found here.

  • output is a numpy array

[21]:

geoms = zip(glaciers.geometry, glaciers.id)
glacier_footprints = features.rasterize(geoms, out_shape=(nrow, ncol), transform=transform)

[22]:

mask = glacier_footprints == 0
masked_footprint = np.ma.masked_array(glacier_footprints, mask=mask)

# plot footprints
fig, ax = plt.subplots(figsize=(8, 8))
im = ax.imshow(masked_footprint)
ax.set_facecolor("lightgrey")
plt.colorbar(im, shrink=0.7);
../../_images/notebooks_part0_python_intro_10_Rasterio_41_0.png

Clip the 1970-2015 raster to get an estimate of total volume change for all glaciers

Using the mask we created before, we can clip the elevation difference arrays to estimate the total volume of change for all of the glaciers

[23]:

masked_diffs = np.ma.masked_array(elevation_diffs["1970-2015"], mask=mask)

[24]:

total_vol_change = masked_diffs.sum() * cellsize ** 2 / 1000**3

fig, ax = plt.subplots(figsize=(8, 8))
im = plt.imshow(masked_diffs, cmap='RdBu', vmin=-80, vmax=80)

ax.text(10, 50, "Total volume change: {:.2f} km$^3$".format(total_vol_change))
ax.set_facecolor("lightgrey")
plt.colorbar(im, label='change in m', shrink=0.6, extend='both');
../../_images/notebooks_part0_python_intro_10_Rasterio_44_0.png

Get volume change by glacier, see if it relates to minimum elevation

[25]:

vol_change = []
for n in glaciers.id.values:
    # compute total volume change
    change_1970_2015 = elevation_diffs['1970-2015'][glacier_footprints == n]
    vc = change_1970_2015.sum() * cellsize ** 2 / 1000**3
    vol_change.append(vc)

glaciers['d_vol_km3'] = vol_change

[26]:

glaciers.sort_values(by='d_vol_km3', inplace=True)
glaciers.head()

[26]:
RGIId GLIMSId BgnDate EndDate CenLon CenLat O1Region O2Region Area Zmin ... Status Connect Form TermType Surging Linkages Name geometry id d_vol_km3
23 RGI60-02.14256 G238227E46901N 19709999 19949999 -121.77327 46.90137 2 4 8.004 1085 ... 0 0 0 0 0 9 Carbon Glacier WA POLYGON ((591944.474 5199551.542, 591934.849 5... 24 -0.112824
111 RGI60-02.14336 G238252E46825N 19709999 19949999 -121.74768 46.82481 2 4 4.444 1425 ... 0 0 0 0 0 1 Nisqually Glacier WA POLYGON ((595661.797 5185442.136, 595671.895 5... 112 -0.071836
26 RGI60-02.14259 G238261E46885N 19709999 19949999 -121.73878 46.88511 2 4 9.195 1437 ... 0 0 0 0 0 9 Winthrop Glacier WA POLYGON ((596567.284 5191536.1, 596500.097 519... 27 -0.069640
64 RGI60-02.14297 G238286E46865N 19709999 19949999 -121.71401 46.86478 2 4 10.594 1605 ... 0 0 0 0 0 1 Emmons Glacier WA POLYGON ((600717.13 5192154.864, 600689.098 51... 65 -0.065486
103 RGI60-02.14329 G238173E46848N 19709999 19949999 -121.82670 46.84820 2 4 4.017 1667 ... 0 0 0 0 0 9 Puyallup Glacier WA POLYGON ((590132.877 5189587.64, 590162.377 51... 104 -0.049973

5 rows × 25 columns

Now we can plot the volume of change vs minimum elevation 'Zmin'

[27]:

fig, ax = plt.subplots(figsize=(11, 8.5))
s = ax.scatter(glaciers.Zmin, glaciers.d_vol_km3, c=glaciers.d_vol_km3, cmap='plasma_r')
for i, r in glaciers.head(10).iterrows():
    ax.text(r.Zmin * 1.03, r.d_vol_km3, r.Name)
ax.set_ylabel('Volume change, 1970-2015, km$^3$')
ax.set_xlabel('Minimum elevation, in meters')
plt.colorbar(s, label='change in km$^3$');
../../_images/notebooks_part0_python_intro_10_Rasterio_49_0.png

We can also plot this relationship spatially using geopandas

[28]:

fig, ax = plt.subplots(figsize=(11, 8.5))

ax.imshow(ls.hillshade(aligned_data[1970], vert_exag=0.5), cmap='gray', extent=extent)
ax = glaciers.plot('d_vol_km3', ax=ax, legend=True, zorder=10, alpha=0.5, cmap="plasma_r")
../../_images/notebooks_part0_python_intro_10_Rasterio_51_0.png

If we regress by area, we see that volume change is coorelated to the original area of the glacier

[29]:

fig, ax = plt.subplots(figsize=(11, 8.5))
s = ax.scatter(glaciers.Area, glaciers.d_vol_km3, c=glaciers.d_vol_km3, cmap='plasma_r')
for i, r in glaciers.head(10).iterrows():
    ax.text(r.Area - 2.03, r.d_vol_km3 - 0.003, r.Name)
ax.set_ylabel('Volume change, 1970-2015, km$^3$')
ax.set_xlabel('Glacial area in km$^2$')
plt.colorbar(s, label='change in km$^3$');
../../_images/notebooks_part0_python_intro_10_Rasterio_53_0.png

Lets plot only the elevation differences in glaciers now

[30]:

just_glaciers = np.ma.masked_array(elevation_diffs['1970-2015'], mask=mask)
fig, ax = plt.subplots(figsize=(11, 8.5))

ax.imshow(
    ls.hillshade(
        aligned_data[1970],
        vert_exag=0.2
    ),
    cmap='gray',
    zorder=-1
)
plt.imshow(just_glaciers, cmap='RdBu', vmin=-30, vmax=30, alpha=0.8)
plt.colorbar(label='Change in glacier elevation, 1970-2015, meters')

[30]:

<matplotlib.colorbar.Colorbar at 0x7fc0663cf850>
../../_images/notebooks_part0_python_intro_10_Rasterio_55_1.png

Writing a raster: save out 1970-2015 elevation changes

To do this, we need to provide rasterio with the same information we got when we read the metadata for the aligned rasters (the elevation changes are on the same grid). * first reset the mask to include all areas outside of the glacier footprints * to simplify things, we can just copy the metadata dict we already have, and feed it to rasterio.open as keyword arguments by unpacking it (** notation) * let’s add compression though to reduce the file size * we can have rasterio write the mask as well * Note that a raster can have more than one band. This code (see ``export_array` function) <https://github.com/aleaf/modflow-export/blob/master/mfexport/array_export.py>`__ illustrates writing a 3D raster with a 3D mask to a single GeoTiff with multiple bands.

[31]:

elevation_diffs['1970-2015'] = np.ma.masked_array(elevation_diffs['1970-2015'], mask=mask)

[32]:

out_meta = meta[2015].copy()
out_meta['compress'] = 'lzw'
out_meta

[32]:

{'driver': 'GTiff',
 'dtype': 'float32',
 'nodata': 0.0,
 'width': 827,
 'height': 761,
 'count': 1,
 'crs': CRS.from_wkt('PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]'),
 'transform': Affine(30.0, 0.0, 583464.253406,
        0.0, -30.0, 5202549.0989),
 'compress': 'lzw'}

[33]:

el_change_outfile = data_path / 'rainier_glaciers_1970-2015_el_change_m.tif'
with rasterio.open(el_change_outfile, 'w', **out_meta) as dest:

    # write the data to band 1
    dest.write(elevation_diffs['1970-2015'], 1)

    # write the mask as well
    dest.write_mask(~elevation_diffs['1970-2015'].mask)

What if we want to reproject to a different CRS?

  • We already have the information we need in memory, but it easiest to simply open a raster with the source grid and then read the metadata from the rasterio file instance

  • We’ll then feed that information to the rasterio.warp.calculate_default_transform function, which will return a new orientation (transform instance) and dimensions for the raster in the destination CRS.

  • We’ll use the new location information to update the metadata from the source grid

  • The code below was taken from the rasterio manual section on reprojection

[34]:

out_meta

[34]:

{'driver': 'GTiff',
 'dtype': 'float32',
 'nodata': 0.0,
 'width': 827,
 'height': 761,
 'count': 1,
 'crs': CRS.from_wkt('PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]'),
 'transform': Affine(30.0, 0.0, 583464.253406,
        0.0, -30.0, 5202549.0989),
 'compress': 'lzw'}

[35]:

from rasterio.crs import CRS
from rasterio.warp import calculate_default_transform, reproject, Resampling

dst_crs = 'epsg:4269'
with rasterio.open(el_change_outfile) as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)

    # copy the source metadata
    kwargs = src.meta.copy()

    # update it with the new location information
    kwargs.update({
    'crs': dst_crs,
    'transform': transform,
    'width': width,
    'height': height
})
name = aligned_rasters[1970].name
name.replace('warp.tif', '4269.tif')
outfile = aligned_rasters[1970].replace(data_path / name)
outfile

[35]:

PosixPath('data/rasterio/aligned-19700901_ned1_2003_adj_warp.tif')

[36]:

with rasterio.open(outfile, 'w', **kwargs) as dst:
    reproject(
        source=aligned_data[1970],
        destination=rasterio.band(dst, 1),
        src_transform=src.transform,
        src_crs=src.crs,
        dst_transform=transform,
        dst_crs=dst_crs,
        resampling=Resampling.bilinear)

Zonal statistics with the rasterstats package

(another way to reduce raster data by polygon features, as we did above in summing the elevation differences by glacier)

[37]:

from rasterstats import zonal_stats

[38]:

help(zonal_stats)

Help on function zonal_stats in module rasterstats.main:

zonal_stats(*args, **kwargs)
    The primary zonal statistics entry point.

    All arguments are passed directly to ``gen_zonal_stats``.
    See its docstring for details.

    The only difference is that ``zonal_stats`` will
    return a list rather than a generator.


[39]:

from rasterstats import gen_zonal_stats
help(gen_zonal_stats)

Help on function gen_zonal_stats in module rasterstats.main:

gen_zonal_stats(vectors, raster, layer=0, band=1, nodata=None, affine=None, stats=None, all_touched=False, categorical=False, category_map=None, add_stats=None, zone_func=None, raster_out=False, prefix=None, geojson_out=False, boundless=True, **kwargs)
    Zonal statistics of raster values aggregated to vector geometries.

    Parameters
    ----------
    vectors: path to an vector source or geo-like python objects

    raster: ndarray or path to a GDAL raster source
        If ndarray is passed, the ``affine`` kwarg is required.

    layer: int or string, optional
        If `vectors` is a path to a fiona source,
        specify the vector layer to use either by name or number.
        defaults to 0

    band: int, optional
        If `raster` is a GDAL source, the band number to use (counting from 1).
        defaults to 1.

    nodata: float, optional
        If `raster` is a GDAL source, this value overrides any NODATA value
        specified in the file's metadata.
        If `None`, the file's metadata's NODATA value (if any) will be used.
        defaults to `None`.

    affine: Affine instance
        required only for ndarrays, otherwise it is read from src

    stats:  list of str, or space-delimited str, optional
        Which statistics to calculate for each zone.
        All possible choices are listed in ``utils.VALID_STATS``.
        defaults to ``DEFAULT_STATS``, a subset of these.

    all_touched: bool, optional
        Whether to include every raster cell touched by a geometry, or only
        those having a center point within the polygon.
        defaults to `False`

    categorical: bool, optional

    category_map: dict
        A dictionary mapping raster values to human-readable categorical names.
        Only applies when categorical is True

    add_stats: dict
        with names and functions of additional stats to compute, optional

    zone_func: callable
        function to apply to zone ndarray prior to computing stats

    raster_out: boolean
        Include the masked numpy array for each feature?, optional

        Each feature dictionary will have the following additional keys:
        mini_raster_array: The clipped and masked numpy array
        mini_raster_affine: Affine transformation
        mini_raster_nodata: NoData Value

    prefix: string
        add a prefix to the keys (default: None)

    geojson_out: boolean
        Return list of GeoJSON-like features (default: False)
        Original feature geometry and properties will be retained
        with zonal stats appended as additional properties.
        Use with `prefix` to ensure unique and meaningful property names.

    boundless: boolean
        Allow features that extend beyond the raster dataset’s extent, default: True
        Cells outside dataset extents are treated as nodata.

    Returns
    -------
    generator of dicts (if geojson_out is False)
        Each item corresponds to a single vector feature and
        contains keys for each of the specified stats.

    generator of geojson features (if geojson_out is True)
        GeoJSON-like Feature as python dict


[40]:

result = zonal_stats(glaciers.geometry.tolist(), el_change_outfile, stats=['mean'], masked=True)

/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/rasterstats/io.py:352: UserWarning: Setting masked to True because dataset mask has been detected
  warnings.warn(

[41]:

result[0:10]

[41]:

[{'mean': -14.126690331304935},
 {'mean': -16.196862951501622},
 {'mean': -7.582302639637432},
 {'mean': -6.186180868474749},
 {'mean': -12.480500358226568},
 {'mean': -6.864627430050872},
 {'mean': -10.58214457460898},
 {'mean': -8.181843650703586},
 {'mean': -9.408220589605067},
 {'mean': -14.84462573412123}]

[42]:

glaciers["mean_el_change_m"] = [d["mean"] for d in result]

[43]:

glaciers.sort_values(by="mean_el_change_m").head()

[43]:
RGIId GLIMSId BgnDate EndDate CenLon CenLat O1Region O2Region Area Zmin ... Connect Form TermType Surging Linkages Name geometry id d_vol_km3 mean_el_change_m
94 RGI60-02.14579 G238389E46854N 19709999 -9999999 -121.61131 46.85395 2 4 0.023 2012 ... 0 0 0 0 9 WA POLYGON ((605930.559 5189945.894, 605935.655 5... 95 -0.001380 -61.327578
111 RGI60-02.14336 G238252E46825N 19709999 19949999 -121.74768 46.82481 2 4 4.444 1425 ... 0 0 0 0 1 Nisqually Glacier WA POLYGON ((595661.797 5185442.136, 595671.895 5... 112 -0.071836 -16.196863
87 RGI60-02.14319 G238159E46860N 19709999 19949999 -121.84114 46.86014 2 4 0.088 1970 ... 0 0 0 0 9 WA POLYGON ((588067.776 5190318.647, 588049.276 5... 88 -0.001280 -14.974627
19 RGI60-02.18818 G238252E46825N 19709999 19949999 -121.74768 46.82481 2 4 1.551 2081 ... 0 0 0 0 9 Wilson Glacier WA POLYGON ((594772.409 5186608.63, 594809.347 51... 20 -0.023033 -14.844626
109 RGI60-02.14591 G238379E46854N 19709999 -9999999 -121.62056 46.85436 2 4 0.405 1867 ... 0 0 0 0 9 Sarvant Glaciers WA POLYGON ((605395.071 5190219.959, 605402.479 5... 110 -0.005892 -14.710493

5 rows × 26 columns

[ ]:



[ ]: