This is the GIS data and imagery used for analyses in the article
Sixty-seven years of land-use change in southern Costa Rica by Zahawi
et al. currently in revision at PLOS One.

This study required the orthorectification of historic aerial photographs, as well as forest cover mapping and landscape analysis of 320 km² around the Las Cruces Biological Station in San Vito de Coto Brus, Costa Rica. The imagery and GIS data generated were used to account for forest cover change over five different time periods from 1947 to 2014.

The datasets supplied include GIS files for:

• Extent of the study area (shapefile).
• Forest cover mapped for each time period (geotiff).
• Imagery of the mosaics generated with the orthorectified historic aerial photographs (geotiff).
• Age in studied time periods of the current forest patches (shapefile).
• Connectivity lines inside the studied area (shapefiles).

All files are in Costa Rica Transverse Mercator 2005 (CRTM05) projected coordinate reference system. For transformation between coordinate systems please refer to http://epsg.io/5367

Aerial photographs for the years 1947, 1960, 1980 and 1997 were acquired from the Organization for Tropical Studies GIS Lab and the Instituto Geográfico Nacional of Costa Rica. The orthorectification process was done first on the 1997 set of images and used the current 1:50,000 and 1:25,000 Costa Rican cartography to identify geographical reference points. The set of 1997 orthophotos was used as a reference set to orthorectify remaining years with the exception of 1947 images. The orthorectification process and all other geospatial analyses were done on the CRTM05 spatial reference system and the resulting orthophotos had a 2m cell size. The largest Root Mean Square error (RMSE) of the orthorectification of these three time slices of aerial photographs was 15 m.

Given the lack of information on flight parameters, and the expansive forest coverage in 1947 photographs, images were georeferenced and built into a mosaic using river basins and the few forest clearings that had a similar shape in the 1960 flyover. The 1947 set of images did not cover the whole study area, having empty areas without photographs that represented ˜12.1% of the analysis extent. Nonetheless, these areas were classified as forested given that forest was present in these same areas in the 1960 imagery.

Forest mapping was done by visual interpretation of orthophotos and Google imagery. The areas were considered forested if tree crowns were easily identified when viewing the images at a scale of 1:10,000. In areas where it was difficult to discern the type of land cover, a scale of 1:5,000 was used. This was done to eliminate agroforestry systems such as shaded coffee areas (with trees planted in rows) or very early stages of forest regeneration from the forest land-cover class. The analysis was done only in areas that were cloud free in the five time slices. This resulted in the elimination of 134 ha (~0.4%) from of the original area outlined above. Polygons were drawn over the different areas using QGIS and were transformed into raster files of 10 m cell size.

Aim: Effective management decisions depend on knowledge of species distribution and habitat use. Maps generated from species distribution models are important in predicting previously unknown occurrences of protected species. However, if populations are seasonally dynamic or locally adapted, failing to consider population level differences could lead to erroneous determinations of occurrence probability and ineffective management. The study goal was to model the distribution of a species of special concern, Townsend’s big-eared bats (Corynorhinus townsendii), in California. We incorporate seasonal and spatial differences to estimate the distribution under current and future climate conditions. Methods: We built species distribution models using all records from statewide roost surveys and by subsetting data to seasonal colonies, representing different phenological stages, and to Environmental Protection Agency Level III Ecoregions to understand how environmental needs vary based on these factors. We projected species’ distribution for 2061-2080 in response to low and high emissions scenarios and calculated the expected range shifts. Results: The estimated distribution differed between the combined (full dataset) and phenologically-explicit models, while ecoregion-specific models were largely congruent with the combined model. Across the majority of models, precipitation was the most important variable predicting the presence of C. townsendii roosts. Under future climate scnearios, distribution of C. townsendii is expected to contract throughout the state, however suitable areas will expand within some ecoregions. Main conclusion: Comparison of phenologically-explicit models with combined models indicate the combined models better predict the extent of the known range of C. townsendii in California. However, life history-explicit models aid in understanding of different environmental needs and distribution of their major phenological stages. Differences between ecoregion-specific and statewide predictions of habitat contractions highlight the need to consider regional variation when forecasting species’ responses to climate change. These models can aid in directing seasonally explicit surveys and predicting regions most vulnerable under future climate conditions. Methods Study area and survey data The study area covers the U.S. state of California, which has steep environmental gradients that support an array of species (Dobrowski et al. 2011). Because California is ecologically diverse, with regions ranging from forested mountain ranges to deserts, we examined local environmental needs by modeling at both the state-wide and ecoregion scale, using U.S. Environmental Protection Agency (EPA) Level III ecoregion designations and there are thirteen Level III ecoregions in California (Table S1.1) (Griffith et al. 2016). Species occurrence data used in this study were from a statewide survey of C. townsendii in California conducted by Harris et al. (2019). Briefly, methods included field surveys from 2014-2017 following a modified bat survey protocol to create a stratified random sampling scheme. Corynorhinus townsendii presence at roost sites was based on visual bat sightings. From these survey efforts, we have visual occurrence data for 65 maternity roosts, 82 hibernation roosts (hibernacula), and 91 active-season non-maternity roosts (transition roosts) for a total of 238 occurrence records (Figure 1, Table S1.1). Ecogeographical factors We downloaded climatic variables from WorldClim 2.0 bioclimatic variables (Fick & Hijmans, 2017) at a resolution of 5 arcmin for broad-scale analysis and 30 arcsec for our ecoregion-specific analyses. To calculate elevation and slope, we used a digital elevation model (USGS 2022) in ArcGIS 10.8.1 (ESRI, 2006). The chosen set of environmental variables reflects knowledge on climatic conditions and habitat relevant to bat physiology, phenology, and life history (Rebelo et al. 2010, Razgour et al. 2011, Loeb and Winters 2013, Razgour 2015, Ancillotto et al. 2016). To trim the global environmental variables to the same extent (the state of California), we used the R package “raster” (Hijmans et al. 2022). We performed a correlation analysis on the raster layers using the “layerStats” function and removed variables with a Pearson’s coefficient > 0.7 (see Table 1 for final model variables). For future climate conditions, we selected three general circulation models (GCMs) based on previous species distribution models of temperate bat species (Razgour et al. 2019) [Hadley Centre Global Environment Model version 2 Earth Systems model (HadGEM3-GC31_LL; Webb, 2019), Institut Pierre-Simon Laplace Coupled Model 6th Assessment Low Resolution (IPSL-CM6A-LR; Boucher et al., 2018), and Max Planck Institute for Meteorology Earth System Model Low Resolution (MPI-ESM1-2-LR; Brovkin et al., 2019)] and two contrasting greenhouse concentration trajectories (Shared Socio-economic Pathways (SSPs): a steady decline pathway with CO2 concentrations of 360 ppmv (SSP1-2.6) and an increasing pathway with CO2 reaching around 2,000 ppmv (SSP5-8.5) (IPCC6). We modeled distribution for present conditions future (2061-2080) time periods. Because one aim of our study was to determine the consequences of changing climate, we changed only the climatic data when projecting future distributions, while keeping the other variables constant over time (elevation, slope). Species distribution modeling We generated distribution maps for total occurrences (maternity + hibernacula + transition, hereafter defined as “combined models”), maternity colonies , hibernacula, and transition roosts. To estimate the present and future habitat suitability for C. townsendii in California, we used the maximum entropy (MaxEnt) algorithm in the “dismo” R package (Hijmans et al. 2021) through the advanced computing resources provided by Texas A&M High Performance Research Computing. We chose MaxEnt to aid in the comparisons of state-wide and ecoregion-specific models as MaxEnt outperforms other approaches when using small datasets (as is the case in our ecoregion-specific models). We created 1,000 background points from random points in the environmental layers and performed a 5-fold cross validation approach, which divided the occurrence records into training (80%) and testing (20%) datasets. We assessed the performance of our models by measuring the area under the receiver operating characteristic curve (AUC; Hanley & McNeil, 1982), where values >0.5 indicate that the model is performing better than random, values 0.5-0.7 indicating poor performance, 0.7-0.9 moderate performance and values of 0.9-1 excellent performance (BCCVL, Hallgren et al., 2016). We also measured the maximum true skill statistic (TSS; Allouche, Tsoar, & Kadmon, 2006) to assess model performance. The maxTSS ranges from -1 to +1:values <0.4 indicate a model that performs no better than random, 0.4-0.55 indicates poor performance, (0.55-0.7) moderate performance, (0.7-0.85) good performance, and values >0.80 indicate excellent performance (Samadi et al. 2022). Final distribution maps were generated using all occurrence records for each region (rather than the training/testing subset), and the models were projected onto present and future climate conditions. Additionally, because the climatic conditions of the different ecoregions of California vary widely, we generated separate models for each ecoregion in an attempt to capture potential local effects of climate change. A general rule in species distribution modeling is that the occurrence points should be 10 times the number of predictors included in the model, meaning that we would need 50 occurrences in each ecoregion. One common way to overcome this limitation is through the ensemble of small models (ESMs) (Breiner et al. 2015., 2018; Virtanen et al. 2018; Scherrer et al. 2019; Song et al. 2019) included in ecospat R package (references). For our ESMs we implemented MaxEnt modeling, and the final ensemble model was created by averaging individual bivariate models by weighted performance (AUC > 0.5). We also used null model significance testing with to evaluate the performance of our ESMs (Raes and Ter Steege 2007). To perform null model testing we compared AUC scores from 100 null models using randomly generated presence locations equal to the number used in the developed distribution model. All ecoregion models outperformed the null expectation (p<0.002). Estimating range shifts For each of the three GCMs and each RCP scenario, we converted the probability distribution map into a binary map (0=unsuitable, 1=suitable) using the threshold that maximizes sensitivity and specificity (Liu et al. 2016). To create the final maps for each SSP scenario, we summed the three binary GCM layers and took a consensus approach, meaning climatically suitable areas were pixels where at least two of the three models predicted species presence (Araújo and New 2007, Piccioli Cappelli et al. 2021). We combined the future binary maps (fmap) and the present binary maps (pmap) following the formula fmap x 2 + pmap (from Huang et al., 2017) to produce maps with values of 0 (areas not suitable), 1 (areas that are suitable in the present but not the future), 2 (areas that are not suitable in the present but suitable in the future), and 3 (areas currently suitable that will remain suitable) using the raster calculator function in QGIS. We then calculated the total area of suitability, area of maintenance, area of expansion, and area of contraction for each binary model using the “BIOMOD_RangeSize” function in R package “biomod2” (Thuiller et al. 2021).

Reason for Selection Protected natural areas in urban environments provide urban residents a nearby place to connect with nature and offer refugia for some species. They help foster a conservation ethic by providing opportunities for people to connect with nature, and also support ecosystem services like offsetting heat island effects (Greene and Millward 2017, Simpson 1998), water filtration, stormwater retention, and more (Hoover and Hopton 2019). In addition, parks, greenspace, and greenways can help improve physical and psychological health in communities (Gies 2006). Urban park size complements the equitable access to potential parks indicator by capturing the value of existing parks.Input DataSoutheast Blueprint 2024 extentFWS National Realty Tracts, accessed 12-13-2023Protected Areas Database of the United States(PAD-US):PAD-US 3.0 national geodatabase -Combined Proclamation Marine Fee Designation Easement, accessed 12-6-20232020 Census Urban Areas from the Census Bureau’s urban-rural classification; download the data, read more about how urban areas were redefined following the 2020 censusOpenStreetMap data “multipolygons” layer, accessed 12-5-2023A polygon from this dataset is considered a beach if the value in the “natural” tag attribute is “beach”. Data for coastal states (VA, NC, SC, GA, FL, AL, MS, LA, TX) were downloaded in .pbf format and translated to an ESRI shapefile using R code. OpenStreetMap® is open data, licensed under theOpen Data Commons Open Database License (ODbL) by theOpenStreetMap Foundation (OSMF). Additional credit to OSM contributors. Read more onthe OSM copyright page.2021 National Land Cover Database (NLCD): Percentdevelopedimperviousness2023NOAA coastal relief model: volumes 2 (Southeast Atlantic), 3 (Florida and East Gulf of America), 4 (Central Gulf of America), and 5 (Western Gulf of America), accessed 3-27-2024Mapping StepsCreate a seamless vector layer to constrain the extent of the urban park size indicator to inland and nearshore marine areas <10 m in depth. The deep offshore areas of marine parks do not meet the intent of this indicator to capture nearby opportunities for urban residents to connect with nature. Shallow areas are more accessible for recreational activities like snorkeling, which typically has a maximum recommended depth of 12-15 meters. This step mirrors the approach taken in the Caribbean version of this indicator.Merge all coastal relief model rasters (.nc format) together using QGIS “create virtual raster”.Save merged raster to .tif and import into ArcPro.Reclassify the NOAA coastal relief model data to assign areas with an elevation of land to -10 m a value of 1. Assign all other areas (deep marine) a value of 0.Convert the raster produced above to vector using the “RasterToPolygon” tool.Clip to 2024 subregions using “Pairwise Clip” tool.Break apart multipart polygons using “Multipart to single parts” tool.Hand-edit to remove deep marine polygon.Dissolve the resulting data layer.This produces a seamless polygon defining land and shallow marine areas.Clip the Census urban area layer to the bounding box of NoData surrounding the extent of Southeast Blueprint 2024.Clip PAD-US 3.0 to the bounding box of NoData surrounding the extent of Southeast Blueprint 2024.Remove the following areas from PAD-US 3.0, which are outside the scope of this indicator to represent parks:All School Trust Lands in Oklahoma and Mississippi (Loc Des = “School Lands” or “School Trust Lands”). These extensive lands are leased out and are not open to the public.All tribal and military lands (“Des_Tp” = "TRIBL" or “Des_Tp” = "MIL"). Generally, these lands are not intended for public recreational use.All BOEM marine lease blocks (“Own_Name” = "BOEM"). These Outer Continental Shelf lease blocks do not represent actively protected marine parks, but serve as the “legal definition for BOEM offshore boundary coordinates...for leasing and administrative purposes” (BOEM).All lands designated as “proclamation” (“Des_Tp” = "PROC"). These typically represent the approved boundary of public lands, within which land protection is authorized to occur, but not all lands within the proclamation boundary are necessarily currently in a conserved status.Retain only selected attribute fields from PAD-US to get rid of irrelevant attributes.Merged the filtered PAD-US layer produced above with the OSM beaches and FWS National Realty Tracts to produce a combined protected areas dataset.The resulting merged data layer contains overlapping polygons. To remove overlapping polygons, use the Dissolve function.Clip the resulting data layer to the inland and nearshore extent.Process all multipart polygons (e.g., separate parcels within a National Wildlife Refuge) to single parts (referred to in Arc software as an “explode”).Select all polygons that intersect the Census urban extent within 0.5 miles. We chose 0.5 miles to represent a reasonable walking distance based on input and feedback from park access experts. Assuming a moderate intensity walking pace of 3 miles per hour, as defined by the U.S. Department of Health and Human Service’s physical activity guidelines, the 0.5 mi distance also corresponds to the 10-minute walk threshold used in the equitable access to potential parks indicator.Dissolve all the park polygons that were selected in the previous step.Process all multipart polygons to single parts (“explode”) again.Add a unique ID to the selected parks. This value will be used in a later step to join the parks to their buffers.Create a 0.5 mi (805 m) buffer ring around each park using the multiring plugin in QGIS. Ensure that “dissolve buffers” is disabled so that a single 0.5 mi buffer is created for each park.Assess the amount of overlap between the buffered park and the Census urban area using “overlap analysis”. This step is necessary to identify parks that do not intersect the urban area, but which lie within an urban matrix (e.g., Umstead Park in Raleigh, NC and Davidson-Arabia Mountain Nature Preserve in Atlanta, GA). This step creates a table that is joined back to the park polygons using the UniqueID.Remove parks that had ≤10% overlap with the urban areas when buffered. This excludes mostly non-urban parks that do not meet the intent of this indicator to capture parks that provide nearby access for urban residents. Note: The 10% threshold is a judgement call based on testing which known urban parks and urban National Wildlife Refuges are captured at different overlap cutoffs and is intended to be as inclusive as possible.Calculate the GIS acres of each remaining park unit using the Add Geometry Attributes function.Buffer the selected parks by 15 m. Buffering prevents very small and narrow parks from being left out of the indicator when the polygons are converted to raster.Reclassify the parks based on their area into the 7 classes seen in the final indicator values below. These thresholds were informed by park classification guidelines from the National Recreation and Park Association, which classify neighborhood parks as 5-10 acres, community parks as 30-50 acres, and large urban parks as optimally 75+ acres (Mertes and Hall 1995).Assess the impervious surface composition of each park using the NLCD 2021 impervious layer and the Zonal Statistics “MEAN” function. Retain only the mean percent impervious value for each park.Extract only parks with a mean impervious pixel value <80%. This step excludes parks that do not meet the intent of the indicator to capture opportunities to connect with nature and offer refugia for species (e.g., the Superdome in New Orleans, LA, the Astrodome in Houston, TX, and City Plaza in Raleigh, NC).Extract again to the inland and nearshore extent.Export the final vector file to a shapefile and import to ArcGIS Pro.Convert the resulting polygons to raster using the ArcPy Feature to Raster function and the area class field.Assign a value of 0 to all other pixels in the Southeast Blueprint 2024 extent not already identified as an urban park in the mapping steps above. Zero values are intended to help users better understand the extent of this indicator and make it perform better in online tools.Use the land and shallow marine layer and “extract by mask” tool to save the final version of this indicator.Add color and legend to raster attribute table.As a final step, clip to the spatial extent of Southeast Blueprint 2024.Note: For more details on the mapping steps, code used to create this layer is available in theSoutheast Blueprint Data Downloadunder > 6_Code. Final indicator valuesIndicator values are assigned as follows:6= 75+ acre urban park5= 50 to <75 acre urban park4= 30 to <50 acre urban park3= 10 to <30 acre urban park2=5 to <10acreurbanpark1 = <5 acre urban park0 = Not identified as an urban parkKnown IssuesThis indicator does not include park amenities that influence how well the park serves people and should not be the only tool used for parks and recreation planning. Park standards should be determined at a local level to account for various community issues, values, needs, and available resources.This indicator includes some protected areas that are not open to the public and not typically thought of as “parks”, like mitigation lands, private easements, and private golf courses. While we experimented with excluding them using the public access attribute in PAD, due to numerous inaccuracies, this inadvertently removed protected lands that are known to be publicly accessible. As a result, we erred on the side of including the non-publicly accessible lands.The NLCD percent impervious layer contains classification inaccuracies. As a result, this indicator may exclude parks that are mostly natural because they are misclassified as mostly impervious. Conversely, this indicator may include parks that are mostly impervious because they are misclassified as mostly

This dataset contains open vector data for railways, forests and power lines, as well an open digital elevation model (DEM) for a small area around a sample forest range in Europe (Germany, Upper Bavaria, Kochel Forest Range, some 70 km south of München, at the edge of Bavarian Alps). The purpose of this dataset is to provide a documented sample dataset in order to demonstrate geospatial preprocessing at FOSS4G2019 based on open data and software. This sample has been produced based on several existing open data sources (detailed below), therefore documenting the sources for obtaining some data needed for computations related to forest accessibility and wood harvesting. For example, they can be used with the open methodology and QGIS plugin Seilaplan for optimising the geometric layout cable roads or with additional open software for computing the forest accessibility for wood harvesting. The vector data (railways, forests and power lines) was extracted from OpenStreetMap (data copyrighted OpenStreetMap contributors and available from https://www.openstreetmap.org). The railways and forests were downloaded and extracted on 18.05.2019 using the open sources QGIS (https://www.qgis.org) with the QuickOSM plugin, while the power lines were downloaded a couple of days later on 23.05.2019. Additional notes for vector data: Please note that OpenStreeMap data extracts such as forests, roads and railways (except power lines) can also be downloaded in a GIS friendly format (Shapefile) from http://download.geofabrik.de/ or using the QGIS built-in download function for OpenStreetMap data. The most efficient way to retrieve specific OSM tags (such as power=line) is to use the QuickOSM plugin for QGIS (using the Overpass API - https://wiki.openstreetmap.org/wiki/Overpass_API) or directly using overpass turbo (https://overpass-turbo.eu/). Finally, the digitised perimeter of the sample forest range is also made available for reproducibility purposes, although any perimeter or area can be digitised freely using the QGIS editing toolbar. The DEM was originally adapted and modified also with QGIS (https://www.qgis.org) based on the elevation data available from two different sources, by reprojecting and downsampling datasets to 25m then selecting, for each individual raster cell, the elevation value that was closer to the average. These two different elevation sources are: - Copernicus Land Monitoring Service - EU-DEM v.1.1 (TILE ID E40N20, downloaded from https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1.1; this original DEM was produced by the Copernicus Land Monitoring Service “with funding by the European Union” based on SRTM and ASTER GDEM) - Digitales Geländemodell 50 m Gitterweite (https://opendata.bayern.de/detailansicht/datensatz/digitales-gelaendemodell-50-m-gitterweite/), produced by the Bayerische Vermessungsverwaltung – www.geodaten.bayern.de –and downloaded from http://www.geodaten.bayern.de/opendata/DGM50/dgm50_epsg4258.tif This methodology was chosen as a way of performing a basic quality check, by comparing the EU-DEM v.1.1 derived from globally available DEM data (such as SRTM) with more authoritative data for the randomly selected region, since using authoritative data is preferred (if open and available). For other sample regions, where authoritative open data is not available, such comparisons cannot longer be performed. Additional notes DEM: a very good DEM open data source for Germany is the open data set collected and resampled by Sonny (sonnyy7@gmail.com) and made available on the Austrian Open Data Portal http://data.opendataportal.at/dataset/dtm-germany. In order to simplify end-to-end reproducibility of the paper planned for FOSS4G2019, we use and distribute an adapted (reprojected and resampled to 25 meters) sample of the above mentioned dataset for the selected forest range. This sample dataset is accompanied by software in Python, as a Jupiter Notebook that generates harmonized output rasters with the same extent from the input data. The extent is given by the polygon vector dataset (Perimeter). These output rasters, such as obstacles, aspect, slope, forest cover, can serve as input data for later computations related to forest accessibility and wood harvesting questions. The obstacles output is obtained by transforming line vector datasets (railway lines, high voltage power lines) to raster. Aspect and slope are both derived from the sample digital elevation model.

NDVI Analysis for Kiambu county. The NDVI was computed using Landsat 8 bands: NDVI=(Band 5(Nir)-Band 4(red))/(Band 5+Band4). The NDVI image was classified into four vegetation health classes. Each class was assigned a specific color using the Qgis symbology tools. The output was a ndvi Kiambu map.Land use Land cover Mapping. Combine all the satellite bands using the band set in the SCP.Landsat image bands were stacked using the Build Virtual Raster tool in QGIS.This formed a multispectral image for further analysis. Clip the image to area of extent. The study area was clipped the map canvas extent. Supervised Classification for LULC Map was performed using the Semi-Automatic Classification Plugin (SCP):Training samples were created for land cover classes: water, vegetation, built-up, bare land.Classification method used: Maximum LikelihoodA classified LULC map was produced, representing each pixel with its respective class.Output: clipped_stack.tif (study area image)

Reason for SelectionHardbottom provides an anchor for important seafloor habitats such as deep-sea corals, plants, and sponges. Hardbottom is also sometimes associated with chemosynthetic communities that form around cold seeps or hydrothermal vents. In these unique ecosystems, micro-organisms that convert chemicals into energy form the base of complex food webs (Love et al. 2013). Hardbottom and associated species provide important habitat structure for many fish and invertebrates (NOAA 2018). Hardbottom areas serve as fish nursery, spawning, and foraging grounds, supporting commercially valuable fisheries like snapper and grouper (NCDEQ 2016).According to Dunn and Halpin (2009), “hardbottom habitats support high levels of biodiversity and are frequently used as a surrogate for it in marine spatial planning.” Artificial reefs arealso known to provide additional habitat that is quickly colonized to provide a suite of ecosystem services commonly associated with naturally occurring hardbottom (Wu et al. 2019). We did not include active oil and gas structures as human-created hardbottom. Although they provide habitat, because of their temporary nature, risk of contamination, and contributions to climate change, they do not have the same level of conservation value as other artificial structures.Input DataSoutheast Blueprint 2024 extentSoutheast Blueprint 2024 subregionsCoral & hardbottomusSEABED Gulf of America sediments, accessed 12-14-2023; download the data; view and read more about the data on the National Oceanic and Atmospheric Administration (NOAA) Gulf Data Atlas (select Physical --> Marine geology --> 1. Dominant bottom types and habitats)Bureau of Ocean Energy Management (BOEM) Gulf of America, seismic water bottom anomalies, accessed 12-20-2023The Nature Conservancy’s (TNC)South Atlantic Bight Marine Assessment(SABMA); chapter 3 of the final report provides more detail on the seafloor habitats analysisNOAA deep-sea coral and sponge locations, accessed 12-20-2023 on the NOAA Deep-Sea Coral & Sponge Map PortalFlorida coral and hardbottom habitats, accessed 12-19-2023Shipwrecks & artificial reefsNOAA wrecks and obstructions layer, accessed 12-12-2023 on the Marine CadastreLouisiana Department of Wildlife and Fisheries (LDWF) Artificial Reefs: Inshore Artificial Reefs, Nearshore Artificial Reefs, Offshore and Deepwater Artificial Reefs (Google Earth/KML files), accessed 12-19-2023Texas Parks and Wildlife Department (TPWD) Artificial Reefs, accessed 12-19-2023; download the data fromThe Artificial Reefs Interactive Mapping Application(direct download from interactive mapping application)Mississippi Department of Marine Resources (MDMR) Artificial Reef Bureau: Inshore Reefs, Offshore Reefs, Rigs to Reef (lat/long coordinates), accessed 12-19-2023Alabama Department of Conservation and Natural Resources (ADCNR) Artificial Reefs: Master Alabama Public Reefs v2023 (.xls), accessed 12-19-2023Florida Fish and Wildlife Conservation Commission (FWC):Artificial Reefs in Florida(.xlsx), accessed 12-19-2023Defining inland extent & split with AtlanticMarine Ecoregions Level III from the Commission for Environmental Cooperation North American Environmental Atlas, accessed 12-8-20212023 NOAA coastal relief model: volumes 2 (Southeast Atlantic), 3 (Florida and East Gulf of America), 4 (Central Gulf of America), and 5 (Western Gulf of America), accessed 3-27-2024National Oceanic and Atmospheric Administration (NOAA)Characterizing Spatial Distributions of Deep-sea Corals and Hardbottom Habitats in the U.S. Southeast Atlantic;read the final report; data shared prior to official release on 2-4-2022 by Matt Poti with the NOAA National Centers for Coastal Ocean Science (NCCOS) (matthew.poti@noaa.gov)Predictive Modeling and Mapping of Hardbottom Seafloor Habitats off the Southeast U.S: unpublished NOAA data anddraft final report entitled Assessment of Benthic Habitats for Fisheries Managementprovided on 1-28-2021 by Matt Poti with NOAA NCCOS (matthew.poti@noaa.gov)Mapping StepsNote: Most of the mapping steps were accomplished using the graphical modeler in QGIS 3.34. Individual models were created to combine data sources and assign ranked values. These models were combined in a single model to assemble all the data sources and create a summary raster. Create a seamless vector layer to constrain the extent of the Atlantic coral and hardbottom indicator to marine and estuarine areas <1 m in elevation. This defines how far inland it extends.Merge together all coastal relief model rasters (.nc format) using the create virtual raster tool in QGIS.Save the merged raster to .tif format and import it into ArcPro.Reclassify the NOAA coastal relief model data to assign a value of 1 to areas from deep marine to 1 m elevation. Assign all other areas (land) a value of 0.Convert the raster produced above to vector using the raster to polygon tool.Clip to the 2024 Blueprint subregions using the pairwise clip tool.Hand-edit to remove terrestrial polygons (one large terrestrial polygon and the Delmarva peninsula).Dissolve the resulting data layer to produce a seamless polygon defining marine and estuarine areas <1 m in elevation.Hand-edit to select all but the main marine polygon and delete.Define the extent of the Gulf version of this indicator to separate it from the Atlantic. This split reflects the extent of the different datasets available to represent coral and hardbottom habitat in the Atlantic and Gulf, rather than a meaningful ecological transition.Use the select tool to select the Florida Keys class from the Level III marine ecoregions (“NAME_L3 = "Florida Keys"“).Buffer the “Florida Keys” Level III marine ecoregion by 2 km to extend it far enough inland to intersect the inland edge of the <1 m elevation layer.Reclassify the two NOAA Atlantic hardbottom suitability datasets to give all non-NoData pixels a value of 0. Combine the reclassified hardbottom suitability datasets to define the total extent of these data. Convert the raster extent to vector and dissolve to create a polygon representing the extent of both NOAA hardbottom datasets.Union the buffered ecoregion with the combined NOAA extent polygon created above. Add a field and use it to dissolve the unioned polygons into one polygon. This leaves some holes inside the polygon, so use the eliminate polygon part tool to fill in those holes, then convert the polygon to a line.Hand-edit to extract the resulting line between the Gulf and Atlantic.Hand-edit to use this line to split the <1 m elevation layer created earlier in the mapping steps to create the separation between the Gulf and Atlantic extent.From the BOEM seismic water bottom anomaly data, extract the following shapefiles: anomaly_confirmed_relic_patchreefs.shp, anomaly_Cretaceous.shp, anomaly_relic_patchreefs.shp, seep_anomaly_confirmed_buried_carbonate.shp, seep_anomaly_confirmed_carbonate.shp, seep_anomaly_confirmed_organisms.shp, seep_anomaly_positives.shp, seep_anomaly_positives_confirmed_gas.shp, seep_anomaly_positives_confirmed_oil.shp, seep_anomaly_positives_possible_oil.shp, seep_anomaly_confirmed_corals.shp, seep_anomaly_confirmed_hydrate.shp.To create a class of confirmed BOEM features, merge anomaly_confirmed_relic_patchreefs.shp, seep_anomaly_confirmed_organisms.shp, seep_anomaly_confirmed_corals.shp, and seep_anomaly_confirmed_hydrate.shp and assign a value of 6.To create a class of predicted BOEM features, merge the remaining extracted shapefiles and assign a value of 3.From usSEABED sediments data, use the field “gom_domnc” to extract polygons: rock (dominant and subdominant) receives a value of 2 and gravel (dominant and subdominant) receives a value of 1.From the wrecks database, extract locations having “high” and “medium” confidence (positionQuality = “high” and positionQuality = “medium”). Buffer these locations by 150 m and assign a value of 4. The buffer distance used here, and later for coral locations, follows guidance from the Army Corps of Engineers for setbacks around artificial reefs and fish havens (Riley et al. 2021).Merge artificial reef point locations from FL, AL, MS and TX. Buffer these locations by 150 m. Merge this file with the three LA artificial reef polygons and assign a value of 5.From the NOAA deep-sea coral and sponge point locations, select all points. Buffer the point locations by 150 m and assign a value of 7.From the FWC coral and hardbottom dataset polygon locations, fix geometries, reproject to EPSG=5070, then assign coral reefs a value of 7, hardbottom a value of 6, hardbottom with seagrass a value of 6, and probable hardbottom a value of 3. Hand-edit to remove an erroneous hardbottom polygon off of Matagorda Island, TX, resulting from a mistake by Sheridan and Caldwell (2002) when they digitized a DOI sediment map. This error is documented on page 6 of the Gulf of Mexico Fishery Management Council’s5-Year Review of the Final Generic Amendment Number 3.From the TNC SABMA data, fix geometries and reproject to EPSG=5070, then select all polygons with TEXT_DESC = "01. mapped hard bottom area" and assign a value of 6.Union all of the above vector datasets together—except the vector for class 6 that combines the SABMA and FL data—and assign final indicator values. Class 6 had to be handled separately due to some unexpected GIS processing issues. For overlapping polygons, this value will represent the maximum value at a given location.Clip the unioned polygon dataset to the buffered marine subregions.Convert both the unioned polygon dataset and the separate vector layer for class 6 using GDAL “rasterize”.Fill NoData cells in both rasters with zeroes and, using Extract by Mask, mask the resulting raster with the Gulf indicator extent. Adding zero values helps users better understand the extent of this indicator and to make this indicator layer perform better in online tools.Use the raster calculator to evaluate the maximum value among

Developed by SOLARGIS and provided by the Global Solar Atlas (GSA), this data resource contains air temperature at 2m above ground level in °C covering the globe. Data is provided in a geographic spatial reference (EPSG:4326). The resolution (pixel size) of solar resource data (GHI, DIF, GTI, DNI) is 9 arcsec (nominally 250 m), PVOUT and TEMP 30 arcsec (nominally 1 km) and OPTA 2 arcmin (nominally 4 km). The data is hyperlinked under 'resources' with the following characteristics: TEMP GISdata (GeoTIFF) Data format: GEOTIFF File size : 121.03 MB There are two temporal representation of solar resource and PVOUT data available: • Longterm yearly/monthly average of daily totals (LTAym_AvgDailyTotals) • Longterm average of yearly/monthly totals (LTAym_YearlyMonthlyTotals) Both type of data are equivalent, you can select the summarization of your preference. The relation between datasets is described by simple equations: • LTAy_YearlyTotals = LTAy_DailyTotals * 365.25 • LTAy_MonthlyTotals = LTAy_DailyTotals * Number_of_Days_In_The_Month For individual country or regional data downloads please see: https://globalsolaratlas.info/download (use the drop-down menu to select country or region of interest) For data provided in AAIGrid please see: https://globalsolaratlas.info/download/world. For more information and terms of use, please, read metadata, provided in PDF and XML format for each data layer in a download file. For other data formats, resolution or time aggregation, please, visit Solargis website. Data can be used for visualization, further processing, and geo-analysis in all mainstream GIS software with raster data processing capabilities (such as open source QGIS, commercial ESRI ArcGIS products and others).

Developed by SOLARGIS and provided by the Global Solar Atlas (GSA), this data resource contains direct normal irradiation (DNI) in kWh/m² covering the globe. Data is provided in a geographic spatial reference (EPSG:4326). The resolution (pixel size) of solar resource data (GHI, DIF, GTI, DNI) is 9 arcsec (nominally 250 m), PVOUT and TEMP 30 arcsec (nominally 1 km) and OPTA 2 arcmin (nominally 4 km). The data is hyperlinked under 'resources' with the following characteristics: DNI LTAy_AvgDailyTotals (GeoTIFF) Data format: GEOTIFF File size : 343.99 MB There are two temporal representation of solar resource and PVOUT data available: • Longterm yearly/monthly average of daily totals (LTAym_AvgDailyTotals) • Longterm average of yearly/monthly totals (LTAym_YearlyMonthlyTotals) Both type of data are equivalent, you can select the summarization of your preference. The relation between datasets is described by simple equations: • LTAy_YearlyTotals = LTAy_DailyTotals * 365.25 • LTAy_MonthlyTotals = LTAy_DailyTotals * Number_of_Days_In_The_Month For individual country or regional data downloads please see: https://globalsolaratlas.info/download (use the drop-down menu to select country or region of interest) For data provided in AAIGrid please see: https://globalsolaratlas.info/download/world. For more information and terms of use, please, read metadata, provided in PDF and XML format for each data layer in a download file. For other data formats, resolution or time aggregation, please, visit Solargis website. Data can be used for visualization, further processing, and geo-analysis in all mainstream GIS software with raster data processing capabilities (such as open source QGIS, commercial ESRI ArcGIS products and others).

Reason for SelectionProtected natural areas in urban environments provide urban residents a nearby place to connect with nature and offer refugia for some species. Because beaches in Puerto Rico and the U.S. Virgin Islands are open to the public, beaches also provide important outdoor recreation opportunities for urban residents, so we include beaches as parks in this indicator.Input DataSoutheast Blueprint 2023 subregions: CaribbeanSoutheast Blueprint 2023 extentNational Oceanic and Atmospheric Administration (NOAA) National Centers for Environmental Information (NCEI) Coastal Relief Model, accessed 11-22-2022Protected Areas Database of the United States (PAD-US) 3.0: VI, PR, and Marine Combined Fee EasementPuerto Rico Protected Natural Areas 2018 (December 2018 update): Terrestrial and marine protected areas (PACAT2018_areas_protegidasPR_TERRESTRES_07052019.shp, PACAT2018_areas_protegidasPR_MARINAS_07052019.shp) 2020 Census Urban Areas from the Census Bureau’s urban-rural classification; download the data, read more about how urban areas were redefined following the 2020 censusOpenStreetMap data “multipolygons” layer, accessed 3-14-2023A polygon from this dataset is considered a park if the “leisure” tag attribute is either “park” or “nature_reserve”, and considered a beach if the value in the “natural” tag attribute is “beach”. OpenStreetMap describes leisure areas as “places people go in their spare time” and natural areas as “a wide variety of physical geography, geological and landcover features”. Data were downloaded in .pbf format and translated ton an ESRI shapefile using R code. OpenStreetMap® is open data, licensed under the Open Data Commons Open Database License (ODbL) by the OpenStreetMap Foundation (OSMF). Additional credit to OSM contributors. Read more on the OSM copyright page. TNC Lands - Public Layer, accessed 3-8-2023U.S. Virgin Islands beaches layer (separate vector layers for St. Croix, St. Thomas, and St. John) provided by Joe Dwyer with Lynker/the NOAA Caribbean Climate Adaptation Program on 3-3-2023 (contact jdwyer@lynker.com for more information)Mapping StepsMost mapping steps were completed using QGIS (v 3.22) Graphical Modeler.Fix geometry errors in the PAD-US PR data using Fix Geometry. This must be done before any analysis is possible.Merge the terrestrial PR and VI PAD-US layers.Use the NOAA coastal relief model to restrict marine parks (marine polygons from PAD-US and Puerto Rico Protected Natural Areas) to areas shallower than 10 m in depth. The deep offshore areas of marine parks do not meet the intent of this indicator to capture nearby opportunities for urban residents to connect with nature.Merge into one layer the resulting shallow marine parks from marine PAD-US and the Puerto Rico Protected Natural Areas along with the combined terrestrial PAD-US parks, OpenStreetMap, TNC Lands, and USVI beaches. Omit from the Puerto Rico Protected Areas layer the “Zona de Conservación del Carso”, which has some policy protections and conservation incentives but is not formally protected.Fix geometry errors in the resulting merged layer using Fix Geometry.Intersect the resulting fixed file with the Caribbean Blueprint subregion.Process all multipart polygons to single parts (referred to in Arc software as an “explode”). This helps the indicator capture, as much as possible, the discrete units of a protected area that serve urban residents.Clip the Census urban area to the Caribbean Blueprint subregion.Select all polygons that intersect the Census urban extent within 1.2 miles (1,931 m). The 1.2 mi threshold is consistent with the average walking trip on a summer day (U.S. DOT 2002) used to define the walking distance threshold used in the greenways and trails indicator. Note: this is further than the 0.5 mi distance used in the continental version of the indicator. We extended it to capture East Bay and Point Udall based on feedback from the local conservation community about the importance of the park for outdoor recreation.Dissolve all the park polygons that were selected in the previous step.Process all multipart polygons to single parts (“explode”) again.Add a unique ID to the selected parks. This value will be used to join the parks to their buffers.Create a 1.2 mi (1,931 m) buffer ring around each park using the multiring buffer plugin in QGIS. Ensure that “dissolve buffers” is disabled so that a single 1.2 mi buffer is created for each park.Assess the amount of overlap between the buffered park and the Census urban area using overlap analysis. This step is necessary to identify parks that do not intersect the urban area, but which lie within an urban matrix. This step creates a table that is joined back to the park polygons using the UniqueID.Remove parks that had ≤2% overlap with the urban areas when buffered. This excludes mostly non-urban parks that do not meet the intent of this indicator to capture parks that provide nearby access for urban residents. Note: In the continental version of this indicator, we used a threshold of 10%. In the Caribbean version, we lowered this to 2% in order to capture small parks that dropped out of the indicator when we extended the buffer distance to 1.2 miles.Calculate the GIS acres of each remaining park unit using the Add Geometry Attributes function.Join the buffer attribute table to the previously selected parks, retaining only the parks that exceeded the 2% urban area overlap threshold while buffered. Buffer the selected parks by 15 m. Buffering prevents very small parks and narrow beaches from being left out of the indicator when the polygons are converted to raster.Reclassify the polygons into 7 classes, seen in the final indicator values below. These thresholds were informed by park classification guidelines from the National Recreation and Park Association, which classify neighborhood parks as 5-10 acres, community parks as 30-50 acres, and large urban parks as optimally 75+ acres (Mertes and Hall 1995).Export the final vector file to a shapefile and import to ArcGIS Pro.Convert the resulting polygons to raster using the ArcPy Polygon to Raster function. Assign values to the pixels in the resulting raster based on the polygon class sizes of the contiguous park areas.Clip to the Caribbean Blueprint 2023 subregion.As a final step, clip to the spatial extent of Southeast Blueprint 2023. Note: For more details on the mapping steps, code used to create this layer is available in the Southeast Blueprint Data Download under > 6_Code. Final indicator valuesIndicator values are assigned as follows:6 = 75+ acre urban park5 = >50 to <75 acre urban park4 = 30 to <50 acre urban park3 = 10 to <30 acre urban park2 = 5 to <10 acre urban park1 = <5 acre urban park0 = Not identified as an urban parkKnown IssuesThis indicator does not include park amenities that influence how well the park serves people and should not be the only tool used for parks and recreation planning. Park standards should be determined at a local level to account for various community issues, values, needs, and available resources. This indicator includes some protected areas that are not open to the public and not typically thought of as “parks”, like mitigation lands, private easements, and private golf courses. While we experimented with excluding them using the public access attribute in PAD, due to numerous inaccuracies, this inadvertently removed protected lands that are known to be publicly accessible. As a result, we erred on the side of including the non-publicly accessible lands.This indicator includes parks and beaches from OpenStreetMap, which is a crowdsourced dataset. While members of the OpenStreetMap community often verify map features to check for accuracy and completeness, there is the potential for spatial errors (e.g., misrepresenting the boundary of a park) or incorrect tags (e.g., labelling an area as a park that is not actually a park). However, using a crowdsourced dataset gives on-the-ground experts, Blueprint users, and community members the power to fix errors and add new parks to improve the accuracy and coverage of this indicator in the future.Other Things to Keep in MindThis indicator calculates the area of each park using the park polygons from the source data. However, simply converting those park polygons to raster results in some small parks and narrow beaches being left out of the indicator. To capture those areas, we buffered parks and beaches by 15 m and applied the original area calculation to the larger buffered polygon, so as not to inflate the area by including the buffer. As a result, when the buffered polygons are rasterized, the final indicator has some areas of adjacent pixels that receive different scores. While these pixels may appear to be part of one contiguous park or suite of parks, they are scored differently because the park polygons themselves are not actually contiguous. The Caribbean version of this indicator uses a slightly different methodology than the continental Southeast version. It includes parks within a 1.2 mi distance from the Census urban area, compared to 0.5 mi in the continental Southeast. We extended it to capture East Bay and Point Udall based on feedback from the local conservation community about the importance of the park for outdoor recreation. Similarly, this indicator uses a 2% threshold of overlap between buffered parks and the Census urban areas, compared to a 10% threshold in the continental Southeast. This helped capture small parks that dropped out of the indicator when we extended the buffer distance to 1.2 miles. Finally, the Caribbean version does not use the impervious surface cutoff applied in the continental Southeast because the landcover data available in the Caribbean does not assess percent impervious in a comparable way.Disclaimer: Comparing with Older Indicator VersionsThere are numerous problems with using Southeast Blueprint