Environmental data for a field site

How to download snow, climate and vegetation data for a field site (SPASS, SLF, MeteoSwiss, ERA5-Land, MODIS).

Author

Arpat Ozgul

Last updated

October 20, 2026

These are the free sources we used for the snow-vole site (downloaded September 2026). Each can be repeated for any site from its coordinates. Four of them cover Switzerland only: SPASS, the SLF 250 m grid, SLF stations and MeteoSwiss stations. ERA5-Land and MODIS cover the world.

The code is in R; Python equivalents are in the folded boxes. All R code was run on 20 October 2026, except the MODIS file search, which NASA’s search service could not answer that day.

Before you start

One point per site, in WGS84, used for every source. The Swiss grids use LV95 (EPSG:2056):

library(sf)
site <- st_sfc(st_point(c(9.572, 46.799)), crs = 4326)   # longitude, latitude
xy   <- st_coordinates(st_transform(site, 2056))          # LV95 easting, northing

Packages: sf, ncdf4, terra, ecmwfr, MODISTools, jsonlite.

Two free accounts:

NASA Earthdata for the MODIS snow files: generate a token and save it as one line in a private file, e.g. ~/.earthdata_token.
Copernicus Climate Data Store for ERA5-Land: copy your API key into ~/.Renviron as ecmwfr_PAT=your-key and restart R. Then accept the terms on the dataset page; requests fail until you do.

Never put a key or token in a script.

Overview

Source	What	Resolution	Account
SPASS (EnviDat)	modelled daily snow water equivalent and depth, 1961–2023	1 km	none
SLF FSM2OSHD (Zenodo)	modelled daily snow depth, SWE, snow-cover fraction, 2016–2025	250 m	none
SLF stations	measured daily snow depth, new snow	points	none
MeteoSwiss stations	measured daily temperature, precipitation, snow	points	none
ERA5-Land	hourly temperature, precipitation, snow, radiation, soil	~9 km	CDS
MODIS MOD13Q1	16-day NDVI	250 m	none (via MODISTools)
MODIS MOD10A1F / MYD10A1F	daily snow index (Terra / Aqua)	500 m	Earthdata

Snow grids

SPASS, 1 km

Two ZIPs from EnviDat, doi:10.16904/envidat.580, each about 1.8 GB and one 9 GB NetCDF once unzipped. Read only your site’s cell:

library(ncdf4)
options(timeout = 3600)
url <- "https://www.envidat.ch/dataset/da0c9ee4-ff02-4848-bc65-cbcaae343636/resource/688d1caf-ed63-4cb8-8579-8c017bb7458c/download/hsclqmd_ch01h.swiss.lv95_wy_1962_2023.zip"
download.file(url, "spass_hs.zip", mode = "wb")
unzip("spass_hs.zip")

nc <- nc_open("HSCLQMD_ch01h.swiss.lv95_WY_1962_2023.nc")
E <- ncvar_get(nc, "E"); N <- ncvar_get(nc, "N")
i <- which.min(abs(E - xy[1])); j <- which.min(abs(N - xy[2]))       # nearest cell
depth <- ncvar_get(nc, "HSCLQMD", start = c(i, j, 1), count = c(1, 1, -1))
date  <- as.Date("1900-01-01") + ncvar_get(nc, "time")
nc_close(nc)
spass <- data.frame(date, snow_depth_m = depth)

Snow water equivalent works the same way: the SWE ZIP, file SWECLQMD_ch01h.swiss.lv95_WY_1962_2023.nc, variable SWECLQMD (in metres).

SWE drops to exactly 0 on day 162 in most years. This is a reset in the product, not melt, so use snow depth ≥ 5 cm to date snow-off.
The depth file has no unit attribute; the values look like metres.

SLF 250 m

One ZIP per water year (2016–2025, about 0.6 GB each) from Zenodo, doi:10.5281/zenodo.17313889, each with 12 monthly NetCDFs:

url <- "https://zenodo.org/api/records/17313889/files/OSHD_DATA_WY2016.zip/content"
download.file(url, "OSHD_DATA_WY2016.zip", mode = "wb")
unzip("OSHD_DATA_WY2016.zip")

read_month <- function(f) {
  nc <- nc_open(f)
  e <- ncvar_get(nc, "easting"); n <- ncvar_get(nc, "northing")
  i <- which.min(abs(e - xy[1])); j <- which.min(abs(n - xy[2]))
  out <- data.frame(time = ncvar_get(nc, "time"))
  for (v in c("hs", "swe", "scf"))        # dimensions here: time, northing, easting
    out[[v]] <- as.numeric(ncvar_get(nc, v, start = c(1, j, i), count = c(-1, 1, 1)))
  nc_close(nc); out
}
slf250 <- do.call(rbind, lapply(list.files(pattern = "^OSHD_DATA_.*nc$"), read_month))

hs is snow depth (m), swe kg/m² (= mm of water), scf snow-cover fraction. Each monthly file counts time in days from its own start (e.g. days since 2016-01-01T06:00:00 UTC+01:00), so convert dates per file with ncatt_get(nc, "time", "units"). This is a model: at our site its snow-off differed from the satellite by about a month.

Python

import xarray as xr
from pyproj import Transformer
e, n = Transformer.from_crs(4326, 2056, always_xy=True).transform(9.572, 46.799)
spass = xr.open_dataset("HSCLQMD_ch01h.swiss.lv95_WY_1962_2023.nc")["HSCLQMD"].sel(E=e, N=n, method="nearest")
slf = xr.open_dataset("OSHD_DATA_2016-01.nc")[["hs", "swe", "scf"]].sel(easting=e, northing=n, method="nearest")

Stations

Stations measure, but not at your site: choose them by distance and elevation, and use them as references.

SLF snow stations:

imis <- read.csv("https://measurement-data.slf.ch/imis/stations.csv")         # automatic stations
plots <- read.csv("https://measurement-data.slf.ch/study-plot/stations.csv")  # manual study plots
st <- st_as_sf(imis, coords = c("lon", "lat"), crs = 4326)
imis$distance_km <- as.numeric(st_distance(st, site)) / 1000
head(imis[order(imis$distance_km), c("station_code", "label", "elevation", "distance_km")])

dav2 <- read.csv("https://measurement-data.slf.ch/imis/data/daily_snow_values/daily_snow_by_station/daily_snow_DAV2.csv")

Manual study plots: https://measurement-data.slf.ch/study-plot/data/by_station/CODE.csv. Some listed stations return an error; log what failed.

MeteoSwiss (networks smn, automatic, and nime, precipitation and snow):

read_ms <- function(u) read.csv2(u, fileEncoding = "Windows-1252")
smn <- read_ms("https://data.geo.admin.ch/ch.meteoschweiz.ogd-smn/ogd-smn_meta_stations.csv")
dav <- read_ms("https://data.geo.admin.ch/ch.meteoschweiz.ogd-smn/dav/ogd-smn_dav_d_historical.csv")

Useful columns: tre200d0 temperature, rre150d0 precipitation, htoautd0 snow depth, hns000d0 new snow. Station codes are lower case in the URL.

ERA5-Land

Hourly values at the grid point nearest your site, one year per request (a month took a few minutes):

library(ecmwfr)          # reads the key from ecmwfr_PAT in ~/.Renviron
for (y in 2005:2025) {
  wf_request(list(
    dataset_short_name = "reanalysis-era5-land-timeseries",
    variable = c("2m_temperature", "total_precipitation", "snow_depth"),
    location = list(longitude = 9.572, latitude = 46.799),
    date = sprintf("%d-01-01/%d-12-31", y, y),
    data_format = "netcdf", target = sprintf("era5_land_%d.zip", y)), path = "era5")
}

Each ZIP holds one NetCDF per variable (e.g. t2m in kelvin). Times are UTC. The grid point’s elevation is not your site’s, so mountain temperatures need care.

Python

import cdsapi
c = cdsapi.Client()   # key in ~/.cdsapirc
c.retrieve("reanalysis-era5-land-timeseries", {
    "variable": ["2m_temperature", "total_precipitation", "snow_depth"],
    "location": {"latitude": 46.799, "longitude": 9.572},
    "date": "2020-01-01/2020-12-31", "data_format": "netcdf"}, "era5_land_2020.zip")

MODIS

NDVI (MOD13Q1) comes directly as a table, no account needed:

library(MODISTools)
ndvi <- mt_subset(product = "MOD13Q1", band = c("250m_16_days_NDVI", "250m_16_days_composite_day_of_the_year",
                  "250m_16_days_pixel_reliability"), lat = 46.799, lon = 9.572,
                  start = "2005-01-01", end = "2025-12-31", km_lr = 0, km_ab = 0, internal = TRUE)

Multiply NDVI by its scale (0.0001). Use each composite’s own day of year, and the reliability band: alpine composites are often snow- or cloud-contaminated.

Snow cover (MOD10A1F, MYD10A1F) is not in MODISTools. For a few sites, NASA’s AppEEARS point extraction gives a CSV without code. For the full files, search NASA’s catalogue and download with your token:

q <- paste0("https://cmr.earthdata.nasa.gov/search/granules.json?collection_concept_id=C3028765772-NSIDC_CPRD",
            "&bounding_box=9.56,46.79,9.58,46.81&temporal=2020-01-01T00:00:00Z,2020-12-31T23:59:59Z&page_size=2000")
hits  <- jsonlite::fromJSON(q)$feed$entry
links <- sapply(hits$links, function(l) grep("^https.*[.]hdf$", l$href, value = TRUE)[1])
token <- readLines("~/.earthdata_token")
for (u in links) download.file(u, basename(u), mode = "wb", headers = c(Authorization = paste("Bearer", token)))

library(terra)
r  <- rast(basename(links[1]), subds = "CGF_NDSI_Snow_Cover")
pt <- project(vect(cbind(9.572, 46.799), crs = "EPSG:4326"), crs(r))
extract(r, pt)

Aqua (MYD10A1F) is C3091257892-NSIDC_CPRD.

The snow value (0–100) is an index, not percent cover; values above 100 are codes for cloud, water or missing.
Cloud-gap-filled days can carry snow from earlier clear days; check the Cloud_Persistence layer.
Season dates need a rule (for example, the first run of 14 snow-free days). Decide it before looking at the demographic data.

Good practice

Keep the original downloads unchanged, in a dated folder per source, on OneDrive (not Git).
Record for each file the source URL and download date; for grids, the cell used and its distance from the site.
Keep download and extraction scripts; never overwrite originals.
Note units and quirks as you find them.