Accessing EMF public spatial products with R
2026-03-11
Víctor Granda-García
Ecosystem Modelling Facility - CREAF
Introduction
Along with the main activities we at the EMF work on, including the development and maintenance of the web apps in the LFC, we generate several data products. This products are made publicly available for anyone to use them.
Datasets available
-
Daily continuous meteorological data (500 m²) for peninsular Spain and Balearic Islands. This data is interpolated from meteorological stations using the
meteolandR package.
Data format isgpkg, containing points in a 500x500 grid in the ETRS89 / UTM zone 30N coordinate reference system (points represent the cell center).
Each file contains ~2 million rows and an avergage file size of 630MB. -
Daily continuous modelled forest water balance for peninsular Spain and Balearic Islands. Data has been generated by the MEDFATE model from soil, climate and forest inputs.
Data format isgpkg, containing points in a 500x500 grid in the ETRS89 / UTM zone 30N coordinate reference system (points represent the cell center).
Each file contains ~760000 rows and an average file size of 175MB
Accessing the data
The simplest way of using the data is downloading the desired files from the
public reposiory (https://data-emf.creaf.cat/public/gpkg/). But because the size
of the files, this can be time, memory and disk storage consuming.
Given the formats offered, in this case we can make use of the capabilities
of the GDAL system library and its
implementation in R by the sf package to avoid this costs.
Simple example
If we only need to access one day, the simplest way is to use sf to read the
file directly from the server.
1# libraries needed
2library(sf)
Linking to GEOS 3.12.2, GDAL 3.11.4, PROJ 9.4.1; sf_use_s2() is TRUE
1library(dplyr)
Attaching package: 'dplyr'
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
1library(mapSpain)
2library(glue)
3
4# retrieving 25th April 2025
5data_20250425 <- st_read(
6 "https://data-emf.creaf.cat/public/gpkg/daily_interpolated_meteo/20250425.gpkg"
7)
Reading layer `20250425' from data source
`https://data-emf.creaf.cat/public/gpkg/daily_interpolated_meteo/20250425.gpkg'
using driver `GPKG'
Simple feature collection with 1988410 features and 24 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -13237.5 ymin: 3988462 xmax: 1125762 ymax: 4858462
Projected CRS: ETRS89 / UTM zone 30N
1# explore result
2data_20250425
Simple feature collection with 1988410 features and 24 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -13237.5 ymin: 3988462 xmax: 1125762 ymax: 4858462
Projected CRS: ETRS89 / UTM zone 30N
First 10 features:
elevation aspect slope TPI partition dates DOY
1 44.1015 309.44631 3.5304249 16.290937 15 2025-04-25 115
2 41.2945 267.64572 1.8394422 0.837438 15 2025-04-25 115
3 29.4370 327.55463 0.9359342 13.133563 15 2025-04-25 115
4 11.0610 65.39845 1.8948245 -5.665062 15 2025-04-25 115
5 37.0690 231.37968 2.3435920 22.689000 15 2025-04-25 115
6 43.3040 193.34568 2.3342511 16.185438 15 2025-04-25 115
7 97.3895 326.60791 4.2136355 58.339062 15 2025-04-25 115
8 39.9230 239.78689 1.2984658 25.387063 15 2025-04-25 115
9 30.6355 90.18747 1.6958913 16.813000 15 2025-04-25 115
10 45.4990 334.25661 5.7842953 -27.030876 15 2025-04-25 115
MeanTemperature MinTemperature MaxTemperature Precipitation
1 14.85347 12.47154 16.40212 1.170740
2 14.87700 12.51194 16.41468 1.171610
3 15.00206 12.77808 16.44801 1.205206
4 15.11709 12.98610 16.50259 1.144229
5 14.91845 12.58773 16.43380 1.172206
6 14.87069 12.49993 16.41208 1.185449
7 14.57880 11.99755 16.25703 1.484663
8 14.93944 12.66118 16.42069 1.244125
9 14.99840 12.76619 16.44970 1.208068
10 14.87326 12.57591 16.36691 1.335326
MeanRelativeHumidity MinRelativeHumidity MaxRelativeHumidity Radiation
1 87.33471 79.08665 100.00000 10.45841
2 87.25151 79.06795 100.00000 10.59436
3 86.56603 78.91359 100.00000 10.29172
4 85.96987 78.67865 98.71368 10.06622
5 87.09487 79.04080 100.00000 10.61602
6 87.33051 79.12048 100.00000 10.75071
7 88.52361 79.48637 100.00000 10.85287
8 86.92086 79.05558 100.00000 10.51207
9 86.63630 78.95054 100.00000 10.35829
10 86.75567 78.83911 100.00000 10.18524
WindSpeed WindDirection PET ThermalAmplitude day month year
1 4.400000 NA 2.174045 3.930580 25 4 2025
2 4.400000 NA 2.203128 3.902745 25 4 2025
3 4.400000 NA 2.171483 3.669932 25 4 2025
4 4.400000 NA 2.176963 3.516493 25 4 2025
5 4.414037 NA 2.215269 3.846070 25 4 2025
6 4.419019 NA 2.231693 3.912151 25 4 2025
7 4.400000 NA 2.193226 4.259476 25 4 2025
8 4.400000 NA 2.198797 3.759505 25 4 2025
9 4.400000 NA 2.182352 3.683511 25 4 2025
10 4.400000 NA 2.130263 3.790991 25 4 2025
geom_hex geom_text
1 010100000000000000e86bf040000000a087615241 POINT (67262.5 4818462)
2 010100000000000000e86bf040000000a00a615241 POINT (67262.5 4817962)
3 010100000000000000d066ee40000000a08d605241 POINT (62262.5 4817462)
4 01010000000000000050a5ee40000000a08d605241 POINT (62762.5 4817462)
5 010100000000000000a84cf040000000a08d605241 POINT (66762.5 4817462)
6 010100000000000000e86bf040000000a08d605241 POINT (67262.5 4817462)
7 010100000000000000d0efec40000000a010605241 POINT (59262.5 4816962)
8 010100000000000000d066ee40000000a010605241 POINT (62262.5 4816962)
9 01010000000000000050a5ee40000000a010605241 POINT (62762.5 4816962)
10 01010000000000000050b1ec40000000a0935f5241 POINT (58762.5 4816462)
geom
1 POINT (67262.5 4818462)
2 POINT (67262.5 4817962)
3 POINT (62262.5 4817462)
4 POINT (62762.5 4817462)
5 POINT (66762.5 4817462)
6 POINT (67262.5 4817462)
7 POINT (59262.5 4816962)
8 POINT (62262.5 4816962)
9 POINT (62762.5 4816962)
10 POINT (58762.5 4816462)
But this is not so different from downloading the data and reading it
locally, as the default behaviour of st_read is to download the file to a
temporal folder and read it from there.
It’s also not so efficient in terms of
time, as we need to wait to download the file and read all the data, which also
means storing all data in memory.
Filtering by polygon
If we don’t need all the data extent (peninsular Spain and Balearic Islands), but
an smaller area instead, we can take advantage of the wkt_filter argument in
st_read.
Given a polygon (in Well Known Text format), st_read uses GDAL to
access and retrieve only the part of the files that intersects with that
polygon. This is useful because reduces transfer size (only a part of the file
is read and transfered), time and memory used.
For this, before reading the remote file, we need to get the polygon of the area
we want. In the following example we’ll use the
mapSpain R package to obtain the
polygon corresponding to Posada de Valdeón (Léon) municipality and use it to
retrieve meteorological data for that polygon from the same file as before:
1posada_wkt <- mapSpain::esp_get_munic(munic = "Posada de Valdeón") |>
2 # ensure the poly is in the same CRS that the data
3 sf::st_transform(crs = 25830) |>
4 # get the polygon geometry
5 dplyr::pull("geometry") |>
6 # transform to text (Well Known Text, wkt)
7 sf::st_as_text()
! The file to be downloaded has size 74.6 Mb.
1# retrieving 25th April 2025 only for Posada
2data_posada_20250425 <- st_read(
3 "https://data-emf.creaf.cat/public/gpkg/daily_interpolated_meteo/20250425.gpkg",
4 wkt_filter = posada_wkt
5)
Reading layer `20250425' from data source
`https://data-emf.creaf.cat/public/gpkg/daily_interpolated_meteo/20250425.gpkg'
using driver `GPKG'
Simple feature collection with 660 features and 24 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 337762.5 ymin: 4772962 xmax: 350262.5 ymax: 4788962
Projected CRS: ETRS89 / UTM zone 30N
1# explore result
2data_posada_20250425
Simple feature collection with 660 features and 24 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 337762.5 ymin: 4772962 xmax: 350262.5 ymax: 4788962
Projected CRS: ETRS89 / UTM zone 30N
First 10 features:
elevation aspect slope TPI partition dates DOY
1 1375.875 283.27479 6.528737 -49.38017 19 2025-04-25 115
2 1434.002 288.97189 9.650236 -52.99203 19 2025-04-25 115
3 1551.576 343.19092 5.846972 36.80266 19 2025-04-25 115
4 1490.074 344.87441 9.266371 -66.80790 19 2025-04-25 115
5 1636.578 311.98249 13.085693 51.79424 19 2025-04-25 115
6 1683.808 311.22276 12.595003 11.19392 19 2025-04-25 115
7 1771.919 337.50742 9.771273 50.60201 19 2025-04-25 115
8 1754.044 32.32764 11.101343 36.49842 19 2025-04-25 115
9 1652.695 357.09260 7.898497 -77.63956 19 2025-04-25 115
10 1722.074 285.44047 12.642417 -48.54472 19 2025-04-25 115
MeanTemperature MinTemperature MaxTemperature Precipitation
1 13.76025 4.815477 19.57583 6.528750
2 13.48020 4.581549 19.26579 6.435703
3 12.90728 4.096716 18.63560 6.249397
4 13.19269 4.359613 18.93565 6.348672
5 12.45904 3.745736 18.12413 6.119528
6 12.20375 3.546991 17.83208 6.048183
7 11.73559 3.170681 17.30419 5.915649
8 11.80104 3.246387 17.36298 5.944145
9 12.30387 3.685165 17.90746 6.100725
10 11.92162 3.384304 17.47229 5.995743
MeanRelativeHumidity MinRelativeHumidity MaxRelativeHumidity Radiation
1 82.50033 57.01226 100 19.39115
2 83.72868 57.92193 100 19.18948
3 86.61378 60.03243 100 19.01024
4 84.70466 58.70404 100 18.40383
5 88.54516 61.53508 100 18.39183
6 89.70303 62.43815 100 18.51165
7 92.15792 64.30271 100 18.57665
8 91.39795 63.81123 100 18.49715
9 88.06260 61.40946 100 18.76046
10 89.93778 62.85809 100 19.35640
WindSpeed WindDirection PET ThermalAmplitude day month year
1 2.580353 NA 3.498776 14.76035 25 4 2025
2 2.616275 NA 3.430047 14.68424 25 4 2025
3 2.645677 NA 3.322786 14.53889 25 4 2025
4 2.669256 NA 3.255552 14.57603 25 4 2025
5 2.687456 NA 3.155307 14.37839 25 4 2025
6 2.700068 NA 3.145606 14.28509 25 4 2025
7 2.698855 NA 3.093759 14.13351 25 4 2025
8 2.688307 NA 3.088102 14.11659 25 4 2025
9 2.673474 NA 3.206316 14.22229 25 4 2025
10 2.656094 NA 3.261536 14.08799 25 4 2025
geom_hex geom_text
1 0101000000000000008a9d1441000000a08f365241 POINT (337762.5 4774462)
2 0101000000000000005aa51441000000a08f365241 POINT (338262.5 4774462)
3 0101000000000000002aad1441000000a08f365241 POINT (338762.5 4774462)
4 010100000000000000fab41441000000a08f365241 POINT (339262.5 4774462)
5 010100000000000000cabc1441000000a08f365241 POINT (339762.5 4774462)
6 0101000000000000009ac41441000000a08f365241 POINT (340262.5 4774462)
7 0101000000000000006acc1441000000a08f365241 POINT (340762.5 4774462)
8 0101000000000000003ad41441000000a08f365241 POINT (341262.5 4774462)
9 0101000000000000000adc1441000000a08f365241 POINT (341762.5 4774462)
10 010100000000000000dae31441000000a08f365241 POINT (342262.5 4774462)
geom
1 POINT (337762.5 4774462)
2 POINT (338262.5 4774462)
3 POINT (338762.5 4774462)
4 POINT (339262.5 4774462)
5 POINT (339762.5 4774462)
6 POINT (340262.5 4774462)
7 POINT (340762.5 4774462)
8 POINT (341262.5 4774462)
9 POINT (341762.5 4774462)
10 POINT (342262.5 4774462)
Important
Be sure that the polygon is in the same Coordinate Reference System (CRS) as the data (EPSG:25380).
wkt_filterdon’t check this and assumes the filter polygon is in the same CRS, so most probably it will return no data when a different CRS is used.
Filtering by other variables (with queries)
Another option we have to reduce the data read into memory is retrieving only
the data we are really interested on. Besides the spatial filter we saw in the
previous section, GDAL also exposes an SQL dialect that can be
used to query only the parts of the data we need. This means we can only read
the columns and rows that comply with the supplied query.
We will use the query argument in st_read to make use of this GDAL feature
and read only temperatures data in Posada de Valdeón for points above
1500m asl:
1# creating the query text
2query_text <- glue::glue(
3 "SELECT
4 geom, elevation, MeanTemperature, MinTemperature, MaxTemperature
5 FROM '20250425'
6 WHERE elevation >= 1500;"
7)
8
9# retrieving 25th April 2025 temperatures only for Posada >= 1500m
10data_posada_1500m_20250425 <- st_read(
11 "https://data-emf.creaf.cat/public/gpkg/daily_interpolated_meteo/20250425.gpkg",
12 wkt_filter = posada_wkt,
13 query = query_text
14)
Reading query `SELECT
geom, elevation, MeanTemperature, MinTemperature, MaxTemperature
FROM '20250425'
WHERE elevation >= 1500;'
from data source `https://data-emf.creaf.cat/public/gpkg/daily_interpolated_meteo/20250425.gpkg'
using driver `GPKG'
Simple feature collection with 356 features and 4 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 337762.5 ymin: 4772962 xmax: 350262.5 ymax: 4787962
Projected CRS: ETRS89 / UTM zone 30N
1# explore result
2data_posada_1500m_20250425
Simple feature collection with 356 features and 4 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 337762.5 ymin: 4772962 xmax: 350262.5 ymax: 4787962
Projected CRS: ETRS89 / UTM zone 30N
First 10 features:
elevation MeanTemperature MinTemperature MaxTemperature
1 1714.730 12.25840 3.757355 17.78548
2 1994.477 11.24679 2.686823 16.81218
3 1865.795 11.71265 3.174516 17.26384
4 1704.162 12.31471 3.793404 17.85496
5 1515.876 13.03160 4.520778 18.56504
6 1845.979 11.54822 3.206440 16.97176
7 2124.161 10.74204 2.180587 16.30839
8 2043.219 11.02119 2.482460 16.57277
9 1973.810 11.26110 2.741995 16.79992
10 1552.001 12.89499 4.374182 18.43493
geom
1 POINT (344262.5 4787962)
2 POINT (342262.5 4787462)
3 POINT (342762.5 4787462)
4 POINT (343262.5 4787462)
5 POINT (343762.5 4787462)
6 POINT (348262.5 4787462)
7 POINT (341762.5 4786962)
8 POINT (342262.5 4786962)
9 POINT (342762.5 4786962)
10 POINT (343262.5 4786962)
Let’s see the query with more detail:
-
SELECT geom, elevation, MeanTemperature, MinTemperature, MaxTemperatureselects the specified data columns. If we want to retrieve all columns, we can useSELECT *, were*means all columns. -
FROM '20250425'the table, in this case, the name of the file without the extension. -
WHERE elevation >= 1500the filters we apply. If we apply more than one filter we need to use logical connectors, likeWHERE elevation >= 1500 AND MinTemperature < 3
If both arguments, query and wkt_filter are used, the geometry column must
be always selected (SELECT geom in this example) or the wkt_filter will not
work.
Tip
For more info about the GDAL SQL dialect, visit https://gdal.org/en/stable/user/ogr_sql_dialect.html
Temporal series (loop)
All the previous examples are using only one file (day). If we want to read more than one, then the simplest way is to make a loop for all the dates we want to work with
1# create the dates to read, but in the format of the files in the server
2dates_to_read <-
3 seq(as.Date("2025-04-25"), as.Date("2025-04-27"), by = "day") |>
4 as.character() |>
5 stringr::str_remove_all("-")
6
7# loop by dates
8data_posada_days <- purrr::map(
9 dates_to_read,
10 \(day) {
11 # query for the day
12 query_text <- glue::glue(
13 "SELECT
14 geom, dates, elevation, MeanTemperature, MinTemperature, MaxTemperature
15 FROM '{day}'
16 WHERE elevation >= 1500;"
17 )
18 # read the posada data of that day
19 st_read(
20 glue::glue("https://data-emf.creaf.cat/public/gpkg/daily_interpolated_meteo/{day}.gpkg"),
21 wkt_filter = posada_wkt,
22 query = query_text
23 )
24 }
25) |>
26 # bind the results in a dataframe
27 purrr::list_rbind()
Reading query `SELECT
geom, dates, elevation, MeanTemperature, MinTemperature, MaxTemperature
FROM '20250425'
WHERE elevation >= 1500;'
from data source `https://data-emf.creaf.cat/public/gpkg/daily_interpolated_meteo/20250425.gpkg'
using driver `GPKG'
Simple feature collection with 356 features and 5 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 337762.5 ymin: 4772962 xmax: 350262.5 ymax: 4787962
Projected CRS: ETRS89 / UTM zone 30N
Reading query `SELECT
geom, dates, elevation, MeanTemperature, MinTemperature, MaxTemperature
FROM '20250426'
WHERE elevation >= 1500;'
from data source `https://data-emf.creaf.cat/public/gpkg/daily_interpolated_meteo/20250426.gpkg'
using driver `GPKG'
Simple feature collection with 356 features and 5 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 337762.5 ymin: 4772962 xmax: 350262.5 ymax: 4787962
Projected CRS: ETRS89 / UTM zone 30N
Reading query `SELECT
geom, dates, elevation, MeanTemperature, MinTemperature, MaxTemperature
FROM '20250427'
WHERE elevation >= 1500;'
from data source `https://data-emf.creaf.cat/public/gpkg/daily_interpolated_meteo/20250427.gpkg'
using driver `GPKG'
Simple feature collection with 356 features and 5 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 337762.5 ymin: 4772962 xmax: 350262.5 ymax: 4787962
Projected CRS: ETRS89 / UTM zone 30N
1# explore result
2head(data_posada_days)
dates elevation MeanTemperature MinTemperature MaxTemperature
1 2025-04-25 1714.730 12.25840 3.757355 17.78548
2 2025-04-25 1994.477 11.24679 2.686823 16.81218
3 2025-04-25 1865.795 11.71265 3.174516 17.26384
4 2025-04-25 1704.162 12.31471 3.793404 17.85496
5 2025-04-25 1515.876 13.03160 4.520778 18.56504
6 2025-04-25 1845.979 11.54822 3.206440 16.97176
geom
1 POINT (344262.5 4787962)
2 POINT (342262.5 4787462)
3 POINT (342762.5 4787462)
4 POINT (343262.5 4787462)
5 POINT (343762.5 4787462)
6 POINT (348262.5 4787462)
