Dataset DOI: 10.5061/dryad.r2280gbq7

Authors: Urquiaga Flores, Erickson Giomar; Kessler, Michael; Körner, Christian; Bader, Maaike Y.

Description of the data and file structure

We present two sets of validation points (table format, csv), three vector datasets delineating our study area (GeoPackages), and four global raster datasets of treeline and upper forest limit elevations (GeoTIFF). Of the rasters, three have a 30-arcsecond resolution and one a 0.27-degree resolution, covering the world’s mountain regions or selected areas within these. The tables have a Comma Separated Values (.csv) format, the vector data a Geopackage format, and the raster data a Geotif format (coordinate system: +proj=longlat +datum=WGS84 +no_defs, EPSG:4326, datum: WGS84).

Files and variables

Authors: Urquiaga Flores, Erickson Giomar; Kessler, Michael; Körner, Christian; Bader, Maaike Y.

Access this dataset on Dryad: Dataset DOI: 10.5061/dryad.r2280gbq7

Sharing and contact information

These data are free to view, edit, and publish in modified forms for scientific use as long as reference is made to the source of the data (this data repository and the research paper Urquiaga Flores et al 2026, OR, for “TREELIM-coordinates_Paulsen_Koerner2014.csv” the research paper Paulsen & Körner 2014), and information about the modifications made is given. The current authors assume no responsibility for any damage caused by the use of these datasets or any derivatives.

Description of the data and file structure

We present two sets of validation points (table format, csv), three vector datasets delineating our study area (geopackages), and four global raster datasets of treeline and upper forest limit elevations (geotif). Of the rasters, three have a 30-arcsecond resolution and one a 0.27-degree resolution, covering the world’s mountain regions or selected areas within these. The tables have a Comma Separated Values (.csv) format, the vector data a Geopackage format, and the raster data a Geotif format (coordinate system: +proj=longlat +datum=WGS84 +no_defs, EPSG:4326, datum: WGS84).

Overview of files:

• TREELIM-coordinates_Paulsen_Koerner2014.csv (376 training points)
• Validation_points_1241_Urquiaga.csv (1241 validation points)
• Regions_new_2025.gpkg (20 polygons of major global regions)
• Study-area-mountains_GMBApolygons.gpkg (486 polygons of mountain ranges)
• Buffer_30km_alpine-250m.gpkg (182 polygons of the alpine zone, defined as treeline minus 250 m elevation, plus a 30-km buffer)
• Treeline_correction_surface.tif (raster used to correct modelled treeline elevation)
• Treeline_elevation_asl_corrected.tif (raster of corrected treeline elevation)
• UFL_pixels_dist_below_corr_tlh.tif (raster of upper forest limit pixels with distance below treeline)
• UFL_pixels_RESAMPLED-0p27degrees.tif (raster with average distance below treeline in 30-arcsecond upper forest limit pixels within 0.27-degree pixels)

Details about the files:

The two tables contain:

1) ”TREELIM-coordinates_Paulsen_Koerner2014.csv”: The geographical location and elevation of 376 points of uppermost tree patches in mountains worldwide, used to train the original TREELIM model by Paulsen & Körner (2014): 
Attributes in this table are:
STAONR: Standortnummer = site ID
STAOTYPE: site type, with two options: GOOGLE_EARTH or LOGGERSTAO (the latter indicating that datalogger recordings of soil temperature are available – data with the original authors)
LATITUDE: in decimal degrees
LONGITUDE: in decimal degrees
ELEVATION: in m a.s.l.
EXPOSURE: discrete (N-S etc) (text)
CONTINENT: name of continent (text)

2) ”Validation_points_1241_Urquiaga.csv”: The geographical location and elevation of 1245 points representing the uppermost patches of trees in different mountain regions, thereby coming closest to the potential treeline. These points were manually selected based on Google Earth imagery. They were used for validation of the potential treeline elevation as predicted by the TREELIM model (Paulsen & Körner 2014) in the implementation by Karger et al. (2019). 
Attributes in this table are:
Fid: site ID number
Point_code: waypoint number
X: longitude in decimal degrees
Y: latitude in decimal degrees
Elev_GoogleEarth: terrain elevation in m above sea level (a.s.l.): elevation according to Google Earth (which is based on the Shuttle, Radar Topographic Mission, SRTM, 90-m rmajortion)
Elev_GMTED2, 010: terrain elevation in m a.s.l.: elevation according to Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010, 30-Arcsecond resolution)
Mountain.Range: name of mountain range
Country: name of country
CHELSA_sd_tlh_1979-2013_s: standard deviation of predicted treeline elevation between years for the Karger et al 2019 TREELIM implementation (m)
CHELSA_tlh_1979-2013_s: predicted treeline elevation above the terrain surface from the Karger et al 2019 TREELIM implementation (m)
CHELSA_asl_tlh_1979-2013: predicted treeline elevation a.s.l. from the Karger et al 2019 TREELIM implementation (m)
dist_below_CHELSA_tlh: distance of the Google-Earth-based elevation of the terrain at the point below the treeline elevation predicted by Karger et al 2019 (m; negative = point lies above the predicted treeline)
corrected_tlh_asl: treeline elevation a.s.l. after we correct the CHELSA-based TREELIM implementation (CHELSA_asl_tlh_1979-2013) areas with trees above the predicted treeline (m). For the correction, we used a B-spline interpolation of negative deviations (i.e,. treeline predicted too low) on the global grid, as explained in Urquiaga Flores et al 2026.
correction_in_m: our correction of treeline elevation (m)(=corrected_tlh_asl - CHELSA_asl_tlh_1979-2013)
for three points in this dataset, the elevation was corrected during processing, since the original point was found to be incorrect upon reinspection. We corrected the elevation of these points (not their geolocation) to the elevation where we found the highest nearby trees (Santa Marta: wpt593, 3930 instead of 4430m; Cocuy, wpt134: 4190m instead of 4545m) or to the elevation of the highest nearby validation point, which was clearly correctly located (NZ Alps: wpt1214, 1120 instead of 1181m).

The three vector layers represent:

1) Regions_new_2025.gpkg: Manually defined regions (n=20) within which treeline and UFL statistics were calculated in Urquiaga Flores et al.
The attributes are:
Fid / id = Region ID number
Region nam = Region name (text)

2) Study-area-mountains_GMBApolygons.gpkg: Selection of mountain ranges, as delineated by the Global Mountain Biodiversity Assessment mountains of the world, version 1.2 (Körner et al., 2017; https://doi.org/10.1007/s00035-016-0182-6). The selection was based on the minimum and maximum distance of the terrain in the range to the corrected potential treeline and the distance to the nearest validation point (see methods below).
The attributes are:
Fid = QGIS-assigned ID number
Name = name of mountain range (text)
Country = name of country or countries (text)

3) Buffer_30km_alpine-250m.gpkg: Alpine zone (i.e. area above the corrected treeline) plus area 250-elevational-m below it plus 30-km buffer (created in ArcGIS in km) around this area. Used to restrict the inclusion of UFL pixels to forest limits near a potential treeline point on the terrain.
The only attribute is:
Fid = QGIS-assigned ID number (meaningless)

The four raster layers represent:

1) The correction surface “Treeline_correction_surface.tif” used to correct the potential treeline elevations. This was based on a selection and spatial interpolation of the elevational differences (in m) between treeline elevation according to the TREELIM model as implemented by Karger et al. (2019) and the elevation of globally-distributed potential-treeline validation points based ;n the regionally highest tree occurrences, representing the observed potential treeline (see the dataset “Validation_points_1241_Urquiaga.csv”, plus, to complete the coverage, we used one training point in the south-east corner of Tierra del Fuego from the Paulsen & Körner 2014 dataset, see “TREELIM-coordinates_ Paulsen&Koerner2014.csv”).
This raster was used as correction surface for the predicted potential treeline elevation in the world’s mountain areas. As presented in Appendix S3 Figure S1b in Urquiaga Flores et al (2026).

2 The elevation of the corrected potential treeline in m above sea level (“Treeline_elevation_asl_corrected.tif”) in the world’s mountain areas. This is the treeline elevation based on the TREELIM model (Karger et al., 2019) corrected with the correction surface (see above). As presented in Appendix S3 Figure S2b in Urquiaga Flores et al (2026).

3a) The elevational distance of the UFL from the corrected potential treeline at a 30-arcsecond resolution (ca. 1 km at the equator). Positive values indicate that the UFL lies above, and negative values that the UFL lies below the potential treeline elevation. (“UFL_pixels_ dist_below_corr_tlh.tif”)

3b) The elevational distance of the UFL from the corrected potential treeline resampled to a 0.27-degree (=972-arcscond, ca. 30 km at the equator) resolution using the mean elevational distance of all 30-arcsecond UFL pixels within the 0.27-degree pixel. Negative values indicate that the UFL lies below the potential treeline. As presented in figure 7 in Urquiaga Flores et al (2026).(“UFL_pixels_RESAMPLED-0p27degrees.tif”)

Below, under Methods, we present more details on the steps taken to generate these data.

Methods

We used QGIS 3.22.10 and SAGAGIS 7.9 or higher (Conrad et al., 2015).

We used the following open-access databases:

· The digital elevation model (DEM) GMTED2010 (Danielson & Gesch, 2011; https://doi.org/10.3133/ofr20111073)

· The global forest cover from Hansen et al. (2013; http://earthenginepartners.appspot.com/science-2013-global-forest)

· The TREELIM model as implemented by Karger et al. (2019; https://envicloud.wsl.ch/#/?prefix=chelsa%2Fchelsa_V1 under climatologies/bio – tlh = treeline height)

· The Global Mountain Biodiversity Assessment (GMBA) mountains of the world, version 1.2 (Körner et al., 2017; https://doi.org/10.1007/s00035-016-0182-6)

· The Lakes and Ocean polygons from the Natural Earth database, version 4.1.0.0 (www.naturalearthdata.com)

And we used the 1241 validation points that we selected along all global mountain areas (presented in this dataset, “Validation_points_1241_Urquiaga.csv”), as well as the 376 points previously collected by Paulsen & Körner (2014) (presented in this dataset, “TREELIM-coordinates_ Paulsen&Koerner2014.csv”).

The following steps were taken to determine the geographical patterns in the accuracy of the TREELIM model (i.e. to validate it by determining the elevational distance between the model and the observations), to generate the corrected treeline elevation prediction, to determine the location and elevation of the upper forest limit (UFL) and the vertical distance of the UFL from the corrected treeline elevation:

Selection of validation points
The validation points were manually selected based on high-resolution Google Earth imagery in 2016-2017. Points were selected where tree cover reaches its maximum elevation within a region (representing either the natural closed forest line or the highest forest patches above the closed forest) and where there are no obvious signs of natural or human disturbance (indicated by nearby settlements and land-use, fire scars, unnaturally straight patterns in the vegetation cover, etc.). A similar approach was followed by Paulsen and Körner (2014) when collecting their 376 points used for developing the TREELIM model. They made their points available to us and they are presented here (with permission).
The results are available in these datasets:
Our (Urquiaga Flores et al 2026) validation points: “Validation_points_1241_Urquiaga.csv”
Paulsen & Körner (2014) validation points: “TREELIM-coordinates_ Paulsen&Koerner2014.csv”

TREELIM validation and correction of potential treeline elevation

1. To validate the predicted global treeline elevation from the CHELSA-based implementation of the TREELIM model (Karger et al., 2019) with observed treeline elevations, we used the 1241 validation points “Validation_points_1241_Urquiaga.csv” and subtracted the predicted potential treeline elevation at that location from their elevation. We used the more precise Google-Earth-based elevation rather than the coarser GMTED2010T-based elevation, since the 1-km resolution of the latter implied a larger deviation in the steep mountain terrain.
2. To create the correction surface, only positive differences (i.e. where the treeline prediction was lower than observed trees) were used, because trees above predicted treeline indicate that the prediction must be wrong, while and absence of trees at treeline (i.e. the highest trees are found below treeline) is quite possible, since many reasons can prevent trees from reaching treeline. Negative differences were therefore set to 0. In addition, we created “pseudopoints” covering all areas less than 2 degrees latitude and 4 degrees longitude from validation points, in a regular grid with points every 2 degree latitude and 4 degree longitude. This served to draw the interpolation surface to 0 in those regions, so as not to extend the correction surface to areas with no information about treeline elevation. We selected the validation points with the regionally (within 1 degree) largest deviations (i.e. with highest tree patches), which were assumed to demonstrate the regional potential of trees to grow at high elevation. We then used a multilevel B-spline interpolation in SAGA GIS, with 10 levels.
   From the resulting correction surface, we subtracted 250 m (to allow some uncertainty related to the DEM quality and positional accuracy of the points). We set <0 values to 0, thus creating a conservative correction surface, correcting only areas with strong underestimations of potential treeline elevation. We then checked all the areas with a correction >0 to see if any were based on a single high-elevation point among lower-elevation points and removed some small correction areas in the Brooks Range, the Scandes in Norway, the Ural Mountains, and Nepal.
   The result is available in this dataset: “Treeline_correction_surface.tif”
3. Then, this correction surface was added to the treeline elevations predicted by the CHELSA-based TREELIM implementation to obtain treeline elevations that were at least as high as the observed trees. The result is available in this dataset: “Treeline_elevation_asl_corrected.tif”

Definition of regions:
To summarize the distances of upper forest limits to treeline in different regions of the world, we manually created polygons representing major regions. These polygons do not cover all parts of the world, but only those with the mountains included in our analyses. They have little use beyond our analysis but are provided to allow replication of our work.
The result is available in this dataset: “Regions_new_2025.gpkg”

Selection of study area mountains
Based on the GMBA Mountain Inventory V 1.2 World

1. We added two Polygons (Chiapas and the Southern Alps of NZ from the supplementary large-scale set)(copy-paste between data sets in QGIS)
2. We deleted mountain ranges with no elevations less than 250m from the corrected potential treeline: assigned minimum distance to potential treeline to each polygon, selected and deleted polygons that have no point less than 250m from potential treeline (min_diff >250).
3. We also deleted ranges with no point below treeline (just some arctic ranges) (max_diff<0)
4. We also deleted ranges that were located far from mountain ranges with validation points, even if close to or even exceeding treeline, e.g. Appalachians, Baffin Island, Greenland, Iceland, Drakensberg, the Atlas Mountains.
   The result is available in this dataset: Study-area-“mountains_GMBApolygons.gpkg”

Creation of 30-km buffer around alpine zone – 250m
We wanted to restrict the analysis of upper forest limits to forest edges that have some level of spatial connection to terrain above or close to potential treeline elevation. We therefore selected all pixels (in the 30-arcsecond resolution) in which the terrain elevation was higher than the corrected potential treeline elevation minus 250m (to allow for some model uncertainty). We turned this area into a polygon layer and created a 30-km buffer (in ArcGIS, where, unlike in QGIS, buffers can be defined in m even in a projection in degrees) around the polygons (including the polygons themselves).
The result is available in this dataset: “Buffer_30km_alpine-250m.gpkg”

Upper forest limit (UFL) position and elevation

1. The 30-m resolution forest presence/absence data from Hansen et al., (2013) data was converted to percentage cover at a 1-km resolution using the mean value method in SAGA (Conrad et al., 2015).
2. To determine a forest-cover threshold to define the treeline, we used a subset of the validation points detailed above (463 points also used in Karger et al., 2019) and determined the mean forest cover at these treeline positions, which was 29%.
3. In the 1-km forest-cover raster, all cells with a forest cover ≥ 29% (forest) were converted to nodata, those <29% (no-forest) to 1.
4. We created a one-pixel buffer around the no-forest area. This was done with the intention that the buffer enters what is defined as forest area and establishes the forest edges (not yet the upper forest limit at this stage), with a value of -1 for each forest-limit cell.
5. Reclassify the raster with values of (1) for all areas except the forest edges. The value of (-1) is exclusive for forest edges.
6. We multiplied this forest-edge layer (map with values of 1 everywhere and -1 only for forest edges) with the DEM GMTED2010. This resulted in a layer with DEM values, positive for all areas except for the forest edges which will have negative DEM values. We will call this layer “DEMplusminus”.

This map could already be used to compare forest-edge elevation to potential treeline elevation, BUT this map includes all forest edges, also lower limits, edges bordering water bodies, and edges of forest gaps. We therefore got rid of all downward-facing forest edges, keeping only those that complied with these criteria: 1) The forest edge pixel should have a direct neighbour (considering 4 neighbours) without forest at a higher elevation; and 2) The forest-edge pixel should not have any pixel with forest in a radius of 2 pixels (considering 12 neighbours) at a higher elevation than itself. Forest borders bordering water bodies were also eliminated.

UFL selection

Condition 1: a direct non-forest neighbour at higher elevation

1. We prepared a forest / no-forest (“FNF1”) map; in the 1-km forest-cover raster, all cells with a forest cover ≥ 29% (forest) as well as all oceans were converted to 0, the rest with <29% (no-forest) to 1.
2. We multiplied FNF with the DEM GMTED2010 to obtain a layer with values of (0) in all forest areas including the oceans, and DEM values for non-forest areas. We will call this layer “DEMnon-forest”.
3. Here the aim was to have negative DEM values (-) on the forest edges, positive DEM values in the non-forest areas, and (0) values in the forest areas. To achieve this, two layers were loaded, DEMplusminus and DEMnon-forest, and the following operation was performed: ifelse(g1≥0,g2, g1), where g1= DEMplusminus and g2= DEMnon-forest. We will call the resulting layer “DEMplusminus_non-forest”.
4. Then we checked for higher-elevation non-forest among the direct neighbours of the forest edges. With DEMplusminus_non-forest a Maximum-Value Focal Statistic, with a ring type kernel with a radius of 1 pixel (i.e. only the 4 direct neighbours), was performed in SAGAGIS. What it does is assign to the focal pixel the maximum value of the neighbouring pixels. We will call the resulting layer “DEMmax_non-forest”. Since only non-forest pixels had elevation values, this allowed verifying that a non-forest pixel was located directly above the forest edge pixels, as should be the case with upper forest edges.
5. To do this check, the two previous layers were subtracted, taking the absolute value for DEMplusminus_non-forest. The operation was: [(DEMmax_non-forest - abs(DEMplusminus_non-forest)]. The result was a map with positive values in upper forest edges (i.e. the maximum non-forest elevation near the forest edge was higher than the edge itself) and negative values in other forest edges, so only edges with positive values were kept.

Condition 2: no forest neighbours at higher elevation

1. We prepared another forest / no-forest (“FNF2”) map; now in the 1-km forest-cover raster, all cells with a forest cover ≥ 29% (forest) were converted to 1, the rest with <29% (no-forest) to 0.
2. We multiplied FNF2 with the DEM GMTED2010 to obtain a layer with values of (0) in all non-forest areas, and DEM values for forest areas. We will call this layer “DEMforest”.
3. Now the aim was to have negative DEM values (-) on the forest edges, values of (0) in the non-forest areas, and positive DEM values (+) in the forest areas. To achieve this, two layers were loaded, DEMplusminus and DEMforest, and the following operation was performed: ifelse(g1>0,g2, g1), where g1= DEMplusminus and g2= DEMforest. We will call the resulting layer “DEMplusminus_forest”.
4. Now we checked for higher-elevation forest around forest edges. With DEMplusminus_forest a Maximum-Value Focal Statistic, with a ring type kernel with a radius of 2 pixels (i.e. including 12 neighbours), was performed in SAGAGIS. What it does is assign to the focal pixel the maximum value of the neighboring pixels. We will call the resulting layer “DEMmax_forest”. Since only forest pixels have elevation values, this allowed checking whether any forest is located within 2 pixels distance above the forest edge pixels, as would be the case with lower forest edges.
5. To do this check, the two previous layers were subtracted, taking the absolute value for DEMplusminus_forest. The operation was: [(DEMmax_forest - abs(DEMplusminus_forest)]. The result was a map with positive values in lower forest edges (i.e. the maximum forest elevation near the forest edge was higher than at the edge) and negative values in upper forest edges, so only edges with negative values were kept.
6. Finally, forest edges bordering water bodies not filtered out by the previous steps were removed using polygons of lakes and ocean from the Natural Earth database.
7. Then, the analysis was restricted to high-enough mountain areas, as presented here in “mountains_GMBApolygons.gpkg”.
8. In addition, a buffer of 30km was created around the intersection of the potential treeline elevation (minus 250m to allow for model uncertainty) with the Earth surface, i.e. the position of the potential treeline (minus 250m) in the mountain landscape (dataset provided here: “Buffer_30km_alpine-250m.gpkg”). Forest edges outside this buffer were not included as UFL, since they are too far away from mountains high enough to support a real treeline and thus do not inform us about the conservation status of subalpine or alpine ecosystems. Often, such edges bordered agricultural fields or other human land uses.

Distance of UFL to potential treeline

1. Finally, for all UFL pixels, the elevation of the pixel (derived from the DEM) was subtracted from the corrected potential treeline elevation (“Corrected_TREELIM_model_in_m_asl.tif”) to obtain the elevational distance between UFL and potential treeline.
   The result is available in this dataset: “Elevational_distances_UFL_vs_Potential_Treeline_1km.tif ”
2. For presentation (Figure 6 in Urquiaga-Flores et al 2025), this map was resampled from the 30 arcsecond to a 0.27-degree (ca 30 m at the equator) resolution as the average elevational distance of the UFL pixels within the larger pixels.
   The result is available in this dataset: “UFL_pixels_RESAMPLED-0p27degrees.tif”

References

1. Conrad, O., B. Bechtel, M. Bock, H. Dietrich, E. Fischer, L. Gerlitz, J. Wehberg, V. Wichmann, and J. Böhner. 2015. System for Automated Geoscientific Analyses (SAGA) v. 2.1.4. Geoscientific Model Development 8:1991–2007. https://doi.org/10.5194/gmd-8-1991-2015
2. Danielson, J.J., Gesch, D. B. 2011. Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010). U.S. Geological Survey Open-File Report 2011-1073 2010:26.
3. Hansen, M. C., P. V. Potapov, R. Moore, M. Hancher, S. A. Turubanova, A. Tyukavina, D. Thau, S. V. Stehman, S. J. Goetz, T. R. Loveland, A. Kommareddy, A. Egorov, L. Chini, C. O. Justice, and J. R. G. Townshend. 2013. High-Resolution Global Maps of 21st-Century Forest Cover Change. Science 342:850–853. https://doi.org/10.1126/science.1244693
4. Karger, D. N., O. Conrad, J. Böhner, T. Kawohl, H. Kreft, R. W. Soria-Auza, N. E. Zimmermann, H. P. Linder, and M. Kessler. 2017. Climatologies at high resolution for the Earth’s land surface areas. Scientific Data 4:170122. https://doi.org/10.1038/sdata.2017.122
5. Karger, D. N., M. Kessler, O. Conrad, P. Weigelt, H. Kreft, C. König, and N. E. Zimmermann. 2019. Why tree lines are lower on islands—Climatic and biogeographic effects hold the answer. Global Ecology and Biogeography 28:839–850. https://doi.org/10.1111/geb.12897
6. Körner, C., W. Jetz, J. Paulsen, D. Payne, K. Rudmann-Maurer, and E. M. Spehn. 2017. A global inventory of mountains for bio-geographical applications. Alpine Botany 127:1–15. https://doi.org/10.1007/s00035-016-0182-6
7. Paulsen, J., and C. Körner. 2014. A climate-based model to predict potential treeline position around the globe. Alpine Botany 124:1–12. https://doi.org/10.1007/s00035-014-0124-0
8. Urquiaga Flores, E.G., M. Kessler, and M.Y. Bader. 2026. Global patterns of treeline predictability and the gap between the forest and the treeline. Ecography https://doi.org/10.1002/ecog.08873