--- title: "Glacier velocity in Python with NASA ITS_LIVE: Jakobshavn Isbræ, 1985–2025" description: "Glacier velocity in Python with NASA ITS_LIVE Zarr datacubes: velocity maps, flowline profiles, filtered speed time series and strain rates." url: https://www.technicalgis.com/tutorials/glacier-velocity-its-live-python/ author: TechnicalGIS Team publisher: TechnicalGIS (https://www.technicalgis.com/about/) published: 2026-10-21 updated: 2026-10-21 language: en-GB --- # Glacier velocity in Python with NASA ITS_LIVE: Jakobshavn Isbræ, 1985–2025 > Glacier velocity in Python with NASA ITS_LIVE Zarr datacubes: velocity maps, flowline profiles, filtered speed time series and strain rates. - Level: Advanced - Tested with: Python 3.12 - Source code: [Download scripts](https://www.technicalgis.com/code/glacier-velocity-its-live-python/) ## Key takeaways 1. ITS\_LIVE publishes glacier velocities for the whole planet as 100 km Zarr datacubes on AWS S3. xarray.open\_dataset(url, engine="zarr") opens one lazily with no login, and a point read of all 80,568 image pairs over Jakobshavn takes about a second. 2. Use the annual composites for maps and profiles, and the image-pair datacube for anything seasonal. The 2014–2024 climatology peaks at 15,195 m/yr, and 598 km² of the cube moves faster than 1 km/yr. 3. A flowline traced through the velocity field runs 89 km up the trunk. 5 km above the ice edge, annual speed rose from 6,208 m/yr in 1985 to 15,276 m/yr in 2007, fell to 9,569 m/yr in 2019 and was 10,489 m/yr in 2024. 4. Raw image pairs need filtering. Keeping pairs at most 48 days apart, with an error below 15% of the speed and within 3 MADs of a 45-day rolling median, kept 72–94% of pairs on the fast trunk; long pairs that fail on the fast ice report about 700 m/yr instead of about 10,000. 5. The lower trunk speeds up by about 50% each summer: from a March–April low to an August–September peak, a range of 49–55% at 5 and 15 km, against 13–14% at 30 and 60 km. ## Before you start - Comfortable with NumPy, pandas and xarray (labelled arrays, .sel and lazy loading) - Python 3.11+ with xarray, zarr 3, dask, fsspec, aiohttp, numpy, pandas, scipy, geopandas, shapely, pyproj, matplotlib and requests - An internet connection: the ITS\_LIVE datacubes are read straight from a public AWS S3 bucket, with no account, API key or download step **Short answer:** to get glacier velocity in Python, use NASA’s **ITS\_LIVE** data. Find the datacube that covers your glacier in the public GeoJSON catalogue, then open it with `xr.open_dataset(zarr_url, engine="zarr")`. No login, API key or download is needed, because the Zarr stores sit in a public AWS S3 bucket. The **annual composites** give clean maps (`v`, `vx`, `vy`) and a 2014–2024 climatology (`v0`). The **image-pair datacube** gives every individual velocity measurement since 1984, which you filter by time separation (`date_dt`) and error (`v_error`) before building time series. Jakobshavn Isbræ (Sermeq Kujalleq in Greenlandic) drains about 6.5% of the Greenland Ice Sheet into Disko Bay through the Ilulissat Icefjord, a UNESCO World Heritage Site. It is one of the fastest glaciers on Earth. After its floating tongue broke up between 1998 and 2003, it roughly doubled in speed (Joughin et al., 2004). It then slowed sharply in 2016–2018, when cooler ocean water reached the fjord (Khazendar et al., 2019). That history makes it an ideal test case: the signal is large, well documented and visible at every timescale from months to decades. You will learn how to: 1. find and open the right **ITS\_LIVE datacube** from the public catalogue, lazily, with xarray; 2. map **surface speed** from the annual composites and the 2014–2024 climatology; 3. trace a **flowline** by integrating the velocity vectors, and extract **speed profiles** along it from 1985 to 2025; 4. pull the **image-pair time series** at points on the flowline and **filter** it by time separation, error and a rolling-median test; 5. measure the **seasonal cycle** and the **multi-year** change, and check them against the composites; 6. compute **strain rates** and plot **flow direction** with a quiver plot. All outputs below are real results from [`glacier_velocity_its_live.py`](https://www.technicalgis.com/code/glacier-velocity-its-live-python/), run on 21 October 2026 with Python 3.12, xarray 2026.9.0, zarr 3.4.1 and scipy 1.18.1. ## How ITS\_LIVE measures glacier velocity ITS\_LIVE measures velocity by **feature tracking**. Take two images of the same glacier a few days or weeks apart. Crevasses, rifts, debris bands and snow dunes move with the ice. For each small window (a **chip**) of the first image, the **autoRIFT** algorithm (Lei et al., 2021) searches the second image for the best match by normalised cross-correlation, and refines it to a fraction of a pixel. The offset in metres, divided by the time between the images, is the surface velocity. Repeating this on a grid gives a velocity field with an x and a y component. ITS\_LIVE runs this on every suitable pair from **Landsat 4–9** and **Sentinel-2** (optical, which needs daylight and no cloud) and **Sentinel-1** (radar amplitude, which works through cloud and in the polar night). Each pair becomes one layer in a datacube, labelled by its **mid\_date** (half way between the two images), with its time separation **date\_dt** in days, error estimates (**v\_error**) and the sensors used. That has three consequences that shape the rest of this tutorial: - **Every measurement is an average** over its time separation. A 12-day Sentinel-1 pair sees the summer speed-up; a 365-day Landsat pair averages over it. - **Matching can fail.** If the ice has changed too much (fast flow, melt, fresh snow) or the offset is larger than the search window, autoRIFT may lock onto stable features and report a speed that is far too low. Long pairs on fast ice are the classic case. - **Errors differ by sensor and by pair.** Filtering on `date_dt` and `v_error` is not optional. ## Step 1: find the datacube that covers the glacier ITS\_LIVE splits the world into 100 km × 100 km Zarr datacubes in local projections: polar stereographic (EPSG:3413) for Greenland. A GeoJSON catalogue lists every cube’s footprint and URL, so finding the right one is a point-in-polygon test with geopandas. `glacier_velocity_its_live.py` ```python import geopandas as gpd import requests import xarray as xr from shapely.geometry import Point LON, LAT = -49.14, 69.10 # a point on the main trunk of Jakobshavn Isbrae, ~25 km above the front CATALOG = "https://its-live-data.s3.amazonaws.com/datacubes/catalog_v02.json" cat_path = DATA / "catalog_v02.json" if not cat_path.exists(): cat_path.write_bytes(requests.get(CATALOG, timeout=300).content) catalog = gpd.read_file(cat_path) hit = catalog[catalog.contains(Point(LON, LAT))].iloc[0] print(f"catalogue: {len(catalog)} datacubes") print(f"cube: {hit.zarr_url}") print(f"composite: {hit.composite_zarr_url}") cube = xr.open_dataset(hit.zarr_url, engine="zarr", chunks={}) # lazy: nothing downloaded yet comp = xr.open_dataset(hit.composite_zarr_url, engine="zarr", chunks={}) ``` `Output — datacube lookup` ```text catalogue: 3160 datacubes cube: https://its-live-data.s3.amazonaws.com/datacubes/v2-updated-october2024/N60W040/ITS_LIVE_vel_EPSG3413_G0120_X-150000_Y-2250000.zarr composite: https://its-live-data.s3.amazonaws.com/composites/annual/v2-updated-september2025/N60W040/ITS_LIVE_velocity_EPSG3413_120m_X-150000_Y-2250000.zarr EPSG:3413, 80568 image pairs datacube: 80568 image pairs on a 833 x 833 grid of 120 m pixels, chunks (20000, 10, 10) mid_date range: 1984-10-21 to 2025-05-28 composites: 42 annual mosaics, 1984-2025 ``` `chunks={}` keeps everything lazy: opening the cube reads only its metadata. Note the chunk shape, **(20000, 10, 10)**. Each chunk holds 20,000 dates for a 10 × 10 pixel block, so the cube is optimised for **time series at a point**: reading every image pair at one pixel touches only five chunks and took about a second. Reading a full map for one date would touch thousands of chunks. That is why maps come from the composites, which are chunked the other way. ## Step 2: map glacier speed from the annual composites The composite store has one mosaic per calendar year (`v`, `vx`, `vy`, `v_error`) and a fitted climatology for 2014–2024: the mean speed `v0`, the components `vx0` and `vy0`, and the amplitude (`v_amp`) and phase (`v_phase`) of a seasonal sine wave. These are small enough to load whole. `glacier_velocity_its_live.py` ```python v0 = comp.v0.load() # 2014-2024 climatological speed (time intercept 2018-01-01) vx0, vy0 = comp.vx0.load(), comp.vy0.load() v_amp, v_phase = comp.v_amp.load(), comp.v_phase.load() landice = comp.landice.load() x, y = comp.x.values, comp.y.values print(f"fastest climatological speed in the cube: {float(v0.max()):,.0f} m/yr") print(f"pixels with v0 > 1 km/yr: {int((v0 > 1000).sum()):,} ({int((v0 > 1000).sum()) * 0.12 ** 2:,.0f} km2)") ``` `Output — annual composites` ```text fastest climatological speed in the cube: 15,195 m/yr pixels with v0 > 1 km/yr: 41,503 (598 km2) ``` A logarithmic colour scale is essential: the ice sheet interior moves at tens of metres a year and the trunk at over 10,000. ![Map of the 2014–2024 ITS_LIVE climatological speed around Jakobshavn Isbræ on a logarithmic colour scale from 10 to over 10,000 m/yr. The fast trunk is a narrow bright channel running west into the fjord, fed by a broad fan of slower ice; the white flowline runs from the ice edge east-south-east for 89 km with points A, B, C and D. An inset map of Greenland marks the cube on the west coast.](https://www.technicalgis.com/_astro/01-jakobshavn-velocity-map.D9yVi5Wp_Z1H8LQG.webp) _The 2014–2024 ITS_LIVE climatology. Speed exceeds 15 km/yr near the ice edge. The white line is the flowline traced in Step 3, and A–D are the points used for time series at 5, 15, 30 and 60 km upstream._ Annual maps show how much the glacier changes from year to year. 2013 was close to the peak of the post-2003 acceleration; 2017 was in the middle of the slowdown. `glacier_velocity_its_live.py` ```python MAP_YEARS = [2013, 2017, 2024] yr = {y_: comp.v.sel(time=f"{y_}-01-01").load() for y_ in MAP_YEARS} im2 = axs[1, 1].imshow(yr[2017] - yr[2013], extent=ext, cmap="RdBu_r", norm=TwoSlopeNorm(0, -4000, 4000)) ``` ![Four map panels zoomed on the lower 50 km of Jakobshavn Isbræ: annual mean speed in 2013, 2017 and 2024 on a log scale, and the difference 2017 minus 2013, which is dark blue (slower by over 4,000 m/yr) along the whole fast trunk and its northern branch and pale elsewhere.](https://www.technicalgis.com/_astro/02-annual-velocity-change.D7TtYCRy_ZNgs7r.webp) _Annual mean speed in 2013, 2017 and 2024, and the change from 2013 to 2017. The slowdown is confined to the fast trunk: it was more than 4 km/yr near the ice edge, and almost nothing on the slow ice either side._ ## Step 3: trace a flowline and extract speed profiles A flowline follows the ice from the ice sheet to the front. We trace it by stepping along the direction of the climatological velocity vectors with a fourth-order Runge–Kutta integrator: upstream until the ice slows below 200 m/yr or we leave the cube, and downstream until the ice reaches the end of the composite’s land-ice mask. Beyond the mask edge is the ice mélange in the fjord, which moves but is not glacier. `glacier_velocity_its_live.py` ```python from scipy import ndimage from scipy.interpolate import RegularGridInterpolator # Smooth the vectors a little, then interpolate them anywhere on the grid. valid = ndimage.gaussian_filter(np.isfinite(vx0.values).astype(float), 2) sx = ndimage.gaussian_filter(np.nan_to_num(vx0.values), 2) / np.maximum(valid, 1e-6) sy = ndimage.gaussian_filter(np.nan_to_num(vy0.values), 2) / np.maximum(valid, 1e-6) sx[valid < 0.5] = np.nan sy[valid < 0.5] = np.nan grid = (y[::-1], x) # RegularGridInterpolator needs ascending axes ix = RegularGridInterpolator(grid, sx[::-1], bounds_error=False, fill_value=np.nan) iy = RegularGridInterpolator(grid, sy[::-1], bounds_error=False, fill_value=np.nan) spd = RegularGridInterpolator(grid, v0.values[::-1], bounds_error=False, fill_value=np.nan) def unit(p): u_, v_ = ix([p[::-1]])[0], iy([p[::-1]])[0] n = np.hypot(u_, v_) return np.array([u_, v_]) / n if np.isfinite(n) and n > 0 else None def trace(start, sign, step=120.0, max_km=150, min_speed=200.0): pts_ = [np.asarray(start, float)] for _ in range(int(max_km * 1000 / step)): p = pts_[-1] k1 = unit(p) k2 = unit(p + sign * step / 2 * k1) if k1 is not None else None k3 = unit(p + sign * step / 2 * k2) if k2 is not None else None k4 = unit(p + sign * step * k3) if k3 is not None else None if k4 is None: break q = p + sign * step / 6 * (k1 + 2 * k2 + 2 * k3 + k4) if not (x[0] < q[0] < x[-1] and y[-1] < q[1] < y[0]) or not spd([q[::-1]])[0] >= min_speed: break pts_.append(q) return np.array(pts_) lon_to = Transformer.from_crs("EPSG:4326", f"EPSG:{hit.epsg}", always_xy=True) seed = np.array(lon_to.transform(LON, LAT)) line = np.vstack([trace(seed, -1)[::-1], trace(seed, +1)[1:]]) # End the flowline where the static land-ice mask ends: beyond it is the ice melange in the fjord. on_ice = landice.sel(x=xr.DataArray(line[:, 0]), y=xr.DataArray(line[:, 1]), method="nearest").values == 1 line = line[: np.flatnonzero(on_ice)[-1] + 1] ``` Where you put the seed matters. Jakobshavn has a northern and a southern branch that join about 15 km above the front. A seed on the lower trunk, traced upstream, follows whichever branch the vectors favour, so check the line on the map. Our seed, on the trunk about 25 km above the front, follows the main southern branch. Then sample every annual composite along the line with vectorised (pointwise) indexing, and drop values whose error is more than 25% of the speed: `glacier_velocity_its_live.py` ```python lx, ly = xr.DataArray(line[:, 0], dims="d"), xr.DataArray(line[:, 1], dims="d") prof = comp.v.sel(x=lx, y=ly, method="nearest").load() prof_err = comp.v_error.sel(x=lx, y=ly, method="nearest").load() prof = prof.where(prof_err < 0.25 * prof) # drop poorly constrained annual values ``` `Output — flowline and annual speeds (m/yr)` ```text flowline: 743 vertices, 89.0 km long; ice-mask edge at x=-196465, y=-2273015; upstream end at x=-114215, y=-2299962 point A: 5 km upstream, x=-191493, y=-2273824, v0=11,810 m/yr point B: 15 km upstream, x=-182234, y=-2276829, v0=10,223 m/yr point C: 30 km upstream, x=-169666, y=-2281882, v0=3,829 m/yr point D: 60 km upstream, x=-142347, y=-2292969, v0=438 m/yr A (5 km) B (15 km) C (30 km) D (60 km) max on line year 1985 6208.0 5527.0 1871.0 NaN 6472.0 1995 5885.0 4435.0 1822.0 NaN 6028.0 2000 10143.0 7583.0 2605.0 NaN 10447.0 2003 14233.0 10616.0 3274.0 363.0 14790.0 2007 15276.0 10554.0 3530.0 376.0 15703.0 2013 14481.0 13767.0 4819.0 438.0 16721.0 2015 13355.0 12647.0 4663.0 455.0 15825.0 2017 13482.0 9451.0 3776.0 452.0 15230.0 2018 11263.0 7997.0 3385.0 431.0 12189.0 2019 9569.0 10523.0 3692.0 419.0 11479.0 2021 13731.0 13190.0 4919.0 429.0 15364.0 2022 10052.0 10866.0 4105.0 423.0 12365.0 2024 10489.0 10436.0 4032.0 420.0 13482.0 2025 11719.0 11291.0 4138.0 415.0 13597.0 ``` ![Left: speed against distance upstream of the flowline end for eight years from 1995 to 2024, plus the 2014–2024 climatology. All curves are below 1,000 m/yr beyond 40 km and rise steeply towards the front; 1995 peaks near 6,000 m/yr, 2007 and 2013 reach 15,000–16,700 m/yr, and 2019 and 2024 sit near 10,000–13,000 m/yr. Right: annual speed at points A to D from 1985 to 2025; A and B about double between 1997 and 2004, peak around 2013, dip in 2018–2019 and recover in 2021.](https://www.technicalgis.com/_astro/03-flowline-profiles.kYcFhlTg_Z2j1tff.webp) _Left: speed along the flowline in selected years. The acceleration reached 40 km inland. Right: annual speed at points A–D. Gaps before 2013 are years in which the composite had no reliable value at that point._ The profiles tell the story of the glacier. In 1995 the lower trunk moved at about 6 km/yr. After the floating tongue broke up, point A sped up to over 14 km/yr by 2003 and 15.3 km/yr by 2007, and the speed-up spread more than 30 km inland: point C went from 1,822 m/yr in 1995 to 4,819 m/yr in 2013. Between 2015 and 2019, point A slowed by 28%, from 13,355 to 9,569 m/yr, and it sped up again in 2021. Point D, 60 km inland, barely changed. ## Step 4: image-pair time series and error filtering For anything shorter than a year you need the image-pair datacube. Select the nearest pixel to each point and load the variables you need. This is where the time-first chunking pays off. `glacier_velocity_its_live.py` ```python VARS = ["v", "v_error", "date_dt", "mission_img1", "satellite_img1"] series = {} for k, p in pts.items(): s = cube[VARS].sel(x=p[0], y=p[1], method="nearest").load().to_dataframe() s = s.dropna(subset=["v"]).sort_index() # mid_date is not sorted in the cube s["sensor"] = s.mission_img1.astype(str) + s.satellite_img1.astype(str) series[k] = s ``` `Output — image pairs at the points` ```text point A: 2,320 image pairs with a value, 1985-03-21 to 2025-05-10 point B: 4,422 image pairs with a value, 1985-04-06 to 2025-05-22 point C: 5,304 image pairs with a value, 1985-03-29 to 2025-05-22 point D: 7,593 image pairs with a value, 1985-06-22 to 2025-05-22 ``` The raw series are noisy. We apply three filters, to pairs from 2014 onwards (when Sentinel-1, Sentinel-2 and Landsat 8 give a dense record): 1. **Time separation**: keep pairs at most 48 days apart. They resolve the seasonal cycle, and they are where feature tracking works best on fast ice. 2. **Error**: keep pairs whose `v_error` is at most 15% of the speed. 3. **Outliers**: keep pairs within 3 median absolute deviations (MAD) of a centred 45-day rolling median, with a floor of 2% of the median. `glacier_velocity_its_live.py` ```python def clean(s: pd.DataFrame, max_dt=48, max_rel_err=0.15, window="45D", n_mad=3.0) -> pd.DataFrame: """Flag short, well-constrained pairs whose speed is close to the local rolling median.""" s = s[s.index >= "2014-01-01"].copy() s["keep_dt"] = s.date_dt <= max_dt s["keep_err"] = s.v_error <= max_rel_err * s.v s1 = s[s.keep_dt & s.keep_err] med = s1.v.rolling(window, center=True, min_periods=5).median() mad = (s1.v - med).abs().rolling(window, center=True, min_periods=5).median() * 1.4826 ok = (s1.v - med).abs() <= n_mad * np.maximum(mad, 0.02 * med) s["keep_mad"] = False s.loc[s1.index[ok.values], "keep_mad"] = True s["keep"] = s.keep_dt & s.keep_err & s.keep_mad return s filt = {k: clean(s) for k, s in series.items()} ``` `Output — filtering` ```text pairs since 2014 dt <= 48 d + error <= 15% + rolling MAD kept % A 1648.0 1493.0 1493.0 1267.0 76.9 B 3873.0 2870.0 2870.0 2784.0 71.9 C 4341.0 4289.0 4191.0 4083.0 94.1 D 7395.0 3027.0 1864.0 1742.0 23.6 sensors among kept pairs at point B: sensor S1A 677 S1B 590 S2B 451 S2A 434 L8 396 L9 236 ``` ![Top: every image-pair speed at point B since 2014, with kept pairs coloured by sensor (Sentinel-1A and 1B, Sentinel-2A and 2B, Landsat 8 and 9) tracing a clear seasonal cycle between about 7,500 and 17,000 m/yr, and rejected pairs in grey, most of them in a band near 700 m/yr. Bottom: days between images on a log scale; kept pairs lie between 5 and 48 days, and nearly all rejected pairs are between about 200 and 500 days apart.](https://www.technicalgis.com/_astro/04-image-pair-filtering.Cbuz_KU2_Z2i5Jx6.webp) _Point B, 15 km above the ice edge. The rejected pairs (grey) are almost all long pairs, 200–500 days apart, that report about 700 m/yr. These are failed matches, not slow ice. The kept pairs from six sensors agree with each other and trace the seasonal cycle._ The figure shows why filtering on `date_dt` matters so much on fast glaciers. Over 200 days, ice moving at 10 km/yr travels more than 5 km. The surface changes beyond recognition and the offset exceeds the search window, so the correlation peak comes from something else and the result is far too slow. A median over all pairs would be badly biased. On slow ice the trade-off reverses: at point D (438 m/yr) the relative-error test removed most pairs, because a `v_error` of a few tens of m/yr is more than 15% of the speed. For slow ice, use an absolute error limit instead. ## Step 5: seasonal cycle and multi-year change With clean pairs, compute monthly medians at each point. For the seasonal cycle, express each month relative to a centred 12-month mean, so that the multi-year trend drops out, then take the median for each calendar month. `glacier_velocity_its_live.py` ```python for k, f in filt.items(): g = f[f.keep] m = g.v.resample("MS").median() trend = m.rolling(13, center=True, min_periods=7).mean() rel = (100 * (m / trend - 1)).groupby(m.index.month).median() ``` `Output — seasonal cycle and multi-year change` ```text km pairs peak month min month seasonal range % 2014 median 2017 median 2024 median point A 5 1267 9 3 54.6 NaN 10025.0 9457.0 B 15 2784 8 4 49.1 11400.0 9114.0 9001.0 C 30 4083 10 5 13.8 4574.0 4000.0 3976.0 D 60 1742 6 11 12.8 444.5 475.0 445.0 ITS_LIVE climatological seasonal fit (2014-2024) at the points: A: v0 11,810 m/yr, amplitude 2,104 m/yr (18%), peak on day 259 (~16 Sep) B: v0 10,223 m/yr, amplitude 2,698 m/yr (26%), peak on day 249 (~06 Sep) C: v0 3,829 m/yr, amplitude 239 m/yr (6%), peak on day 267 (~24 Sep) D: v0 438 m/yr, amplitude 2 m/yr (0%), peak on day 172 (~21 Jun) ``` ![Left: filtered image-pair speeds and monthly median lines at points A to D from 2014 to 2025. A and B swing every year between lows near 8,000 m/yr in spring and peaks of 12,000–19,000 m/yr in late summer, with the highest peaks in 2021; C varies between about 3,500 and 5,000 m/yr and D is flat near 450 m/yr. Right: median departure from the 12-month mean by calendar month; A and B are about 15–19% below average in February to April and 30–37% above in August and September, while C and D vary by less than 10%.](https://www.technicalgis.com/_astro/05-seasonal-multiyear-speed.BSWSxcmC_Z1B8Uqz.webp) _Left: filtered pairs (dots) and monthly medians (lines). Right: the median seasonal cycle. The lower trunk is slowest in March–April and fastest in August–September, a range of about 50%. The cycle fades upstream._ The lower trunk is fastest in late summer and slowest in spring, a range of 55% at point A and 49% at point B. ITS\_LIVE’s own sine-wave fit agrees on the timing (peaks on about 16 and 6 September) but gives a smaller amplitude, 18–26% of the mean. That is expected: Jakobshavn’s cycle is not a sine wave. It speeds up sharply in early summer, when the front retreats and the mélange that buttresses it breaks up, and slows gradually through the winter as the front advances again. Upstream, the cycle is much weaker: 6% in the ITS\_LIVE fit at point C, and nothing measurable at point D. ## Step 6: strain rates and flow direction Spatial derivatives of the velocity field give the **strain rate**: how fast the ice stretches or compresses. Derivatives amplify noise, so smooth the climatological components first. We use a Gaussian with σ = 5 pixels (600 m), a few times the ice thickness. The **longitudinal strain rate** is the stretching rate along the flow direction. The **effective strain rate** is the second invariant of the strain-rate tensor, with the vertical component taken from incompressibility. It is the quantity that enters Glen’s flow law for ice. `glacier_velocity_its_live.py` ```python dx, dy = float(x[1] - x[0]), float(y[1] - y[0]) # dy is negative: y decreases down the rows w = ndimage.gaussian_filter(np.isfinite(vx0.values).astype(float), 5) u = ndimage.gaussian_filter(np.nan_to_num(vx0.values), 5) / np.maximum(w, 1e-6) v = ndimage.gaussian_filter(np.nan_to_num(vy0.values), 5) / np.maximum(w, 1e-6) mask = (w > 0.9) & (landice.values == 1) dudy, dudx = np.gradient(u, dy, dx) dvdy, dvdx = np.gradient(v, dy, dx) exy = 0.5 * (dudy + dvdx) speed = np.hypot(u, v) e_long = (u ** 2 * dudx + u * v * (dudy + dvdx) + v ** 2 * dvdy) / speed ** 2 # along-flow strain rate, 1/yr e_eff = np.sqrt(0.5 * (dudx ** 2 + dvdy ** 2 + (dudx + dvdy) ** 2) + exy ** 2) # effective strain rate, 1/yr ``` `Output — strain rates` ```text point A: longitudinal strain rate +0.015 /yr, effective 0.394 /yr point B: longitudinal strain rate +0.720 /yr, effective 0.871 /yr point C: longitudinal strain rate +0.291 /yr, effective 0.950 /yr point D: longitudinal strain rate +0.009 /yr, effective 0.015 /yr max longitudinal stretching on the flowline: 1.12 /yr at 12.4 km max compression on the flowline: -1.09 /yr at 0.0 km ``` For flow direction, plot unit vectors with `quiver`, so that arrow length does not swamp the slow ice: `glacier_velocity_its_live.py` ```python st = 25 # one arrow every 25 pixels = 3 km X, Y = np.meshgrid(x[::st], y[::st]) U, V = u[::st, ::st], v[::st, ::st] S = np.hypot(U, V) q = (S > 100) & mask[::st, ::st] axs[1].quiver(X[q], Y[q], (U / S)[q], (V / S)[q], color="white", scale=32, width=0.003, headwidth=4) ``` ![Left: map of longitudinal strain rate, from blue (compression) to red (stretching), saturating at ±0.6 per year. A band of strong stretching runs along the trunk from about 25 km to 8 km above the ice edge and up the northern branch, and there is strong compression near the ice edge. Right: speed on a log scale with white unit arrows every 3 km, converging from a wide area of the ice sheet into the trunk and the fjord.](https://www.technicalgis.com/_astro/06-strain-rate-flow-direction.D_E7jiBg_1rmbfJ.webp) _Left: longitudinal strain rate from the 2014–2024 climatology. Strain rates exceed 1 per year where the ice accelerates into the trunk, and compressed near the ice edge. Right: flow direction, with unit arrows every 3 km. The basin funnels ice from a wide area into a channel only a few kilometres wide._ Ice entering the trunk at point B is stretched at 0.72 per year: a 100 m length of ice grows by about 20 cm every day. The maximum, 1.12 per year, is 12.4 km above the ice edge, where the trunk is accelerating hardest. Rates like these open the dense crevasse fields that make the lower glacier impassable. Close to the ice edge the flow is compressive (−1.09 per year), because the front is buttressed by the mélange and the fjord walls. Point C has modest stretching along the flow (0.29 per year) but a high effective strain rate (0.95 per year), so other components, mainly shear, dominate there. ## Step 7: save the results `glacier_velocity_its_live.py` ```python out = pd.DataFrame({"dist_km": dist_km, "x": line[:, 0], "y": line[:, 1], "v0": spd(line[:, ::-1]), "e_long": along}) out = out.join(pd.DataFrame(prof.values.T, columns=[f"v_{y_}" for y_ in years])) out.to_csv(OUT / "flowline_profiles.csv", index=False, float_format="%.3f") pd.concat({k: f[f.keep][["v", "v_error", "date_dt", "sensor"]] for k, f in filt.items()}, names=["point"]).to_csv( OUT / "point_timeseries_filtered.csv", float_format="%.2f") ``` `flowline_profiles.csv` has one row per flowline vertex, with its distance, coordinates in EPSG:3413, climatological speed, strain rate and the annual speed for every year from 1984 to 2025. `point_timeseries_filtered.csv` has every kept image pair at points A–D. To use the flowline in QGIS, write `x` and `y` as points with geopandas in EPSG:3413. ## Pitfalls and limitations - **ITS\_LIVE measures surface velocity.** It says nothing directly about the bed or the ice flux. To get discharge you need ice thickness, for example from BedMachine Greenland. - **Failed matches on fast ice.** Long image pairs, pairs across the melt season, and pairs over fresh snow can report speeds that are far too low. Filter by `date_dt` and look at the scatter before you trust a mean. - **Relative error filters penalise slow ice.** At point D only 24% of pairs survived the 15% test. Use an absolute error limit for slow-moving ice, or the `v_error` weighting used in the composites. - **Uneven sampling.** Sentinel-1B failed in December 2021, and optical sensors see nothing in the polar night. Monthly medians, averaged month by month, give fairer annual figures than the median of all pairs. - **Partial years and a moving front.** The 2025 composite was built from data up to May 2025 only. The ice edge also moves by kilometres each season, so a pixel near the front may be glacier in one image and mélange or open water in the next. We stopped the flowline at the static land-ice mask for that reason. - **Smoothing controls strain rates.** Peak strain rates depend strongly on the smoothing length: less smoothing gives higher, noisier peaks. Report the smoothing with any strain-rate figure. - **Datacube versions.** ITS\_LIVE reprocesses and updates its cubes. The catalogue pointed to `v2-updated-october2024` cubes with data up to May 2025. Re-running later may return more pairs and slightly different numbers. ## Frequently asked questions ### How do I get glacier velocity data in Python? Use NASA ITS\_LIVE. Download the public catalogue (`catalog_v02.json`), find the datacube whose footprint contains your glacier, and open its `zarr_url` with `xr.open_dataset(url, engine="zarr")`. No account is needed. For maps, open the matching annual composite (`composite_zarr_url`). The `itslive` Python package wraps the same steps if you prefer a higher-level interface. ### What is the difference between the ITS\_LIVE datacube and the composites? The datacube holds every individual image-pair velocity field: tens of thousands of layers, each an average over its own time separation, chunked for fast time series at a point. The annual composites are one error-weighted mosaic per year, plus a fitted 2014–2024 climatology with a seasonal amplitude and phase. Use composites for maps and profiles, and the datacube for seasonal and short-term change. ### How does feature tracking measure glacier velocity? It matches small image chips between two images of the same place taken days to months apart. The best match, found by normalised cross-correlation and refined to sub-pixel precision, gives the displacement. Divided by the time between the images, that is the velocity. ITS\_LIVE uses the autoRIFT algorithm on Landsat, Sentinel-2 and Sentinel-1 pairs. ### Why do some ITS\_LIVE velocities look far too slow? Usually because matching failed. On fast glaciers, image pairs that are hundreds of days apart see the surface change beyond recognition, so the correlation locks onto something else. At Jakobshavn such pairs report about 700 m/yr where the ice moves at about 10,000 m/yr. Filter by time separation (`date_dt`), error (`v_error`) and a rolling-median test. ### How fast does Jakobshavn Isbræ flow? In the 2014–2024 ITS\_LIVE climatology, the fastest ice in the cube moves at 15,195 m/yr, about 42 m per day. Annual speed 5 km above the ice edge rose from 6,208 m/yr in 1985 to 15,276 m/yr in 2007, fell to 9,569 m/yr in 2019 and was 10,489 m/yr in 2024. In summer, filtered image pairs 15 km upstream exceed 17,000 m/yr. ### Can I use ITS\_LIVE for mountain glaciers such as the Aletsch? Yes. The datacubes cover glaciers worldwide, in UTM projections outside the polar regions, and the same code works if you change `LON` and `LAT`. Mountain glaciers are slower and narrower than Jakobshavn, so use an absolute error threshold, check that the 120 m grid resolves the glacier, and expect more gaps in winter from snow and shadow. ## Next steps - Map snow, the other half of the cryosphere, in [snow cover mapping with NDSI in Python](https://www.technicalgis.com/tutorials/snow-cover-ndsi-sentinel-2-python/). - Measure millimetre-scale ground movement with radar interferometry in [InSAR ground deformation in Python](https://www.technicalgis.com/tutorials/insar-ground-deformation-python/). Feature tracking is the large-displacement counterpart of InSAR. - Compute slope and aspect of the ice surface from a DEM with [raster terrain analysis with rasterio](https://www.technicalgis.com/tutorials/raster-terrain-analysis-rasterio/). - Build gridded time series from STAC instead of Zarr with [NDVI time series in Python with Sentinel-2 and odc-stac](https://www.technicalgis.com/tutorials/ndvi-time-series-odc-stac/). - Watch the calving front itself: the water-edge methods in [coastal erosion and shoreline change from Sentinel-2](https://www.technicalgis.com/tutorials/coastal-erosion-shoreline-sentinel-2/) and [Sentinel-1 SAR flood mapping](https://www.technicalgis.com/tutorials/sentinel-1-sar-flood-mapping/) work on the fjord too. - Put the flowline on an interactive map with [Folium](https://www.technicalgis.com/tutorials/interactive-web-maps-folium/), or browse more open data in the [TechnicalGIS data catalogue](https://www.technicalgis.com/data-sources/). ## Data and licences - **Glacier velocities**: velocity data generated using auto-RIFT (Gardner et al., 2018; Lei et al., 2021) and provided by the NASA MEaSUREs ITS\_LIVE project (Gardner et al., [doi:10.5067/JQ6337239C96](https://doi.org/10.5067/JQ6337239C96)). ITS\_LIVE image-pair datacube `ITS_LIVE_vel_EPSG3413_G0120_X-150000_Y-2250000.zarr` (v2, updated October 2024, data to 28 May 2025) and annual composite `ITS_LIVE_velocity_EPSG3413_120m_X-150000_Y-2250000.zarr` (v2, updated September 2025), read from the public `its-live-data` AWS S3 bucket on 21 October 2026. ITS\_LIVE data are openly available with no restrictions on use; the attribution follows the project’s recommended wording at [its-live.jpl.nasa.gov](https://its-live.jpl.nasa.gov/). - **Underlying imagery**: Landsat courtesy of the U.S. Geological Survey; contains modified Copernicus Sentinel-1 and Sentinel-2 data, processed by ITS\_LIVE. - **Country outlines**: Natural Earth, public domain. - References: Gardner, A. et al., _MEaSUREs ITS\_LIVE Regional Glacier and Ice Sheet Surface Velocities, Version 2_ (NSIDC-0776), NASA National Snow and Ice Data Center DAAC, [doi:10.5067/JQ6337239C96](https://doi.org/10.5067/JQ6337239C96); Gardner, A. S. et al. (2018), “Increased West Antarctic and unchanged East Antarctic ice discharge over the last 7 years”, _The Cryosphere_ 12(2), 521–547, [doi:10.5194/tc-12-521-2018](https://doi.org/10.5194/tc-12-521-2018); Lei, Y., Gardner, A. and Agram, P. (2021), “Autonomous Repeat Image Feature Tracking (autoRIFT) and its application for tracking ice displacement”, _Remote Sensing_ 13(4), 749; Joughin, I., Abdalati, W. and Fahnestock, M. (2004), “Large fluctuations in speed on Greenland’s Jakobshavn Isbræ glacier”, _Nature_ 432, 608–610; Khazendar, A. et al. (2019), “Interruption of two decades of Jakobshavn Isbrae acceleration and thinning as regional ocean cools”, _Nature Geoscience_ 12, 277–283. - Code: [`glacier_velocity_its_live.py`](https://www.technicalgis.com/code/glacier-velocity-its-live-python/), tested with Python 3.12.0, xarray 2026.9.0, zarr 3.4.1, dask 2026.8.0, numpy 2.5.4, pandas 3.0.6, scipy 1.18.1, geopandas 1.2.0 and matplotlib 3.11.2. ## Key terms - **ITS_LIVE**: Inter-mission Time Series of Land Ice Velocity and Elevation, a NASA MEaSUREs project at JPL that measures glacier and ice-sheet surface velocity from Landsat 4–9, Sentinel-1 and Sentinel-2 image pairs and publishes it openly. - **Feature tracking**: Measuring how far the surface moved between two images by finding where a small patch (a chip) of the first image best matches the second, usually by normalised cross-correlation. Displacement divided by the time between the images gives the velocity. - **autoRIFT**: Autonomous Repeat Image Feature Tracking (Lei et al., 2021), the open-source feature-tracking code that ITS\_LIVE uses for both optical and radar image pairs. - **Datacube**: Here, a Zarr store holding every image-pair velocity field for one 100 km tile, stacked along a mid\_date dimension, with per-pair metadata such as the time separation (date\_dt), error estimates and sensor. - **Flowline**: A line that follows the direction of ice flow, traced here by stepping along the velocity vectors from a seed point. Profiles along it show how speed changes from the ice sheet to the calving front. - **Strain rate**: The spatial rate of change of velocity, in 1/yr. The longitudinal strain rate along flow is positive where the ice stretches (and crevasses open) and negative where it is compressed. ## Data for this tutorial - [NASA ITS\_LIVE](https://its-live.jpl.nasa.gov/) · Image-pair velocity datacubes and annual composites as Zarr on AWS S3 - [ITS\_LIVE datacube catalogue](https://its-live-data.s3.amazonaws.com/datacubes/catalog_v02.json) · GeoJSON footprint and URL of every datacube - [Natural Earth](https://www.technicalgis.com/data-sources/#src-natural-earth) · Greenland outline for the locator inset Licences and formats for each source are listed in the [GIS data sources catalogue](https://www.technicalgis.com/data-sources/). --- Source: https://www.technicalgis.com/tutorials/glacier-velocity-its-live-python/