Introduction

We are enclosing the database used in our research titled "Concentration and Geospatial Modelling of Health Development Offices' Accessibility for the Total and Elderly Populations in Hungary", along with our statistical calculations. For the sake of reproducibility, further information can be found in the file Short_Description_of_Data_Analysis.pdf and Statistical_formulas.pdf

The sharing of data is part of our aim to strengthen the base of our scientific research. As of March 7, 2024, the detailed submission and analysis of our research findings to a scientific journal has not yet been completed.

The dataset was expanded on 23rd September 2024 to include SPSS statistical analysis data, a heatmap, and buffer zone analysis around the Health Development Offices (HDOs) created in QGIS software.

Short Description of Data Analysis and Attached Files (datasets):

Our research utilised data from 2022, serving as the basis for statistical standardisation. The 2022 Hungarian census provided an objective basis for our analysis, with age group data available at the county level from the Hungarian Central Statistical Office (KSH) website. The 2022 demographic data provided an accurate picture compared to the data available from the 2023 microcensus. The used calculation is based on our standardisation of the 2022 data. For xlsx files, we used MS Excel 2019 (version: 1808, build: 10406.20006) with the SOLVER add-in.

Hungarian Central Statistical Office served as the data source for population by age group, county, and regions: https://www.ksh.hu/stadat_files/nep/hu/nep0035.html, (accessed 04 Jan. 2024.) with data recorded in MS Excel in the Data_of_demography.xlsx file.

In 2022, 108 Health Development Offices (HDOs) were operational, and it's noteworthy that no developments have occurred in this area since 2022. The availability of these offices and the demographic data from the Central Statistical Office in Hungary are considered public interest data, freely usable for research purposes without requiring permission.

The contact details for the Health Development Offices were sourced from the following page (Hungarian National Population Centre (NNK)): https://www.nnk.gov.hu/index.php/efi (n=107). The Semmelweis University Health Development Centre was not listed by NNK, hence it was separately recorded as the 108th HDO. More information about the office can be found here: https://semmelweis.hu/egeszsegfejlesztes/en/ (n=1). (accessed 05 Dec. 2023.)

Geocoordinates were determined using Google Maps (N=108): https://www.google.com/maps. (accessed 02 Jan. 2024.) Recording of geocoordinates (latitude and longitude according to WGS 84 standard), address data (postal code, town name, street, and house number), and the name of each HDO was carried out in the: Geo_coordinates_and_names_of_Hungarian_Health_Development_Offices.csv file.

The foundational software for geospatial modelling and display (QGIS 3.34), an open-source software, can be downloaded from:

https://qgis.org/en/site/forusers/download.html. (accessed 04 Jan. 2024.)

The HDOs_GeoCoordinates.gpkg QGIS project file contains Hungary's administrative map and the recorded addresses of the HDOs from the

Geo_coordinates_and_names_of_Hungarian_Health_Development_Offices.csv file,

imported via .csv file.

The OpenStreetMap tileset is directly accessible from www.openstreetmap.org in QGIS. (accessed 04 Jan. 2024.)

The Hungarian county administrative boundaries were downloaded from the following website: https://data2.openstreetmap.hu/hatarok/index.php?admin=6 (accessed 04 Jan. 2024.)

HDO_Buffers.gpkg is a QGIS project file that includes the administrative map of Hungary, the county boundaries, as well as the HDO offices and their corresponding buffer zones with a radius of 7.5 km.

Heatmap.gpkg is a QGIS project file that includes the administrative map of Hungary, the county boundaries, as well as the HDO offices and their corresponding heatmap (Kernel Density Estimation).

A brief description of the statistical formulas applied is included in the Statistical_formulas.pdf.

Recording of our base data for statistical concentration and diversification measurement was done using MS Excel 2019 (version: 1808, build: 10406.20006) in .xlsx format.

• Aggregated number of HDOs by county: Number_of_HDOs.xlsx
• Standardised data (Number of HDOs per 100,000 residents): Standardized_data.xlsx
• Calculation of the Lorenz curve: Lorenz_curve.xlsx
• Calculation of the Gini index: Gini_Index.xlsx
• Calculation of the LQ index: LQ_Index.xlsx
• Calculation of the Herfindahl-Hirschman Index: Herfindahl_Hirschman_Index.xlsx
• Calculation of the Entropy index: Entropy_Index.xlsx
• Regression and correlation analysis calculation: Regression_correlation.xlsx

Using the SPSS 29.0.1.0 program, we performed the following statistical calculations with the databases Data_HDOs_population_without_outliers.sav and Data_HDOs_population.sav:

• Regression curve estimation with elderly population and number of HDOs, excluding outlier values (Types of analyzed equations: Linear, Logarithmic, Inverse, Quadratic, Cubic, Compound, Power, S, Growth, Exponential, Logistic, with summary and ANOVA analysis table): Curve_estimation_elderly_without_outlier.spv
• Pearson correlation table between the total population, elderly population, and number of HDOs per county, excluding outlier values such as Budapest and Pest County: Pearson_Correlation_populations_HDOs_number_without_outliers.spv.
• Dot diagram including total population and number of HDOs per county, excluding outlier values such as Budapest and Pest Counties: Dot_HDO_total_population_without_outliers.spv.
• Dot diagram including elderly (64<) population and number of HDOs per county, excluding outlier values such as Budapest and Pest Counties: Dot_HDO_elderly_population_without_outliers.spv
• Regression curve estimation with total population and number of HDOs, excluding outlier values (Types of analyzed equations: Linear, Logarithmic, Inverse, Quadratic, Cubic, Compound, Power, S, Growth, Exponential, Logistic, with summary and ANOVA analysis table): Curve_estimation_without_outlier.spv
• Dot diagram including elderly (64<) population and number of HDOs per county: Dot_HDO_elderly_population.spv
• Dot diagram including total population and number of HDOs per county: Dot_HDO_total_population.spv
• Pearson correlation table between the total population, elderly population, and number of HDOs per county: Pearson_Correlation_populations_HDOs_number.spv
• Regression curve estimation with total population and number of HDOs, (Types of analyzed equations: Linear, Logarithmic, Inverse, Quadratic, Cubic, Compound, Power, S, Growth, Exponential, Logistic, with summary and ANOVA analysis table): Curve_estimation_total_population.spv

For easier readability, the files have been provided in both SPV and PDF formats.

The translation of these supplementary files into English was completed on 23rd Sept. 2024.

If you have any further questions regarding the dataset, please contact the corresponding author: domjan.peter@phd.semmelweis.hu

Urban Tree Canopy Assessment. This was created using the Urban Tree Canopy Syracuse 2010 (All Layers) file HERE.The data for this map was created using LIDAR and other spatial analysis tools to identify and measure tree canopy in the landscape. This was a collaboration between the US Forest Service Northern Research Station (USFS), the University of Vermont Spatial Laboratory, and SUNY ESF. Because the full map is too large to be viewed in ArcGIS Online, this has been reduced to a vector tile layer to allow it to be viewed online. To download and view the shapefiles and all of the layers, you can download the data HERE and view this in either ArcGIS Pro or QGIS.Data DictionaryDescription source  USDA Forest ServiceList of values  Value 1 Description Tree CanopyValue 2 Description Grass/ShrubValue 3 Description Bare SoilValue 4 Description WaterValue 5 Description BuildingsValue 6 Description Roads/RailroadsValue 7 Description Other PavedField Class Alias Class Data type String Width 20Geometric objects  Feature class name landcover_2010_syracusecity Object type  complex Object count 7ArcGIS Feature Class Properties Feature class name landcover_2010_syracusecity Feature type  Simple Geometry type Polygon Has topology FALSE Feature count 7 Spatial index TRUE Linear referencing  FALSEDistributionAvailable format  Name ShapefileTransfer options  Transfer size 163.805Description Downloadable DataFieldsDetails for object landcover_2010_syracusecityType Feature Class Row count  7 Definition  UTCField FIDAlias FID Data type OID Width  4 Precision 0 Scale 0Field descriptionInternal feature number.Description source ESRIDescription of valueSequential unique whole numbers that are automatically generated.Field ShapeAlias Shape Data type Geometry Width 0 Precision 0 Scale 0Field description Feature geometry.Description source  ESRIDescription of values Coordinates defining the features.Field CodeAlias Code Data type Number Width 4Overview Description  Metadata DetailsMetadata language  English Metadata character set utf8 - 8 bit UCS Transfer FormatScope of the data described by the metadata  dataset Scope name  datasetLast update 2011-06-02ArcGIS metadata properties Metadata format ArcGIS 1.0 Metadata style North American Profile of ISO19115 2003Created in ArcGIS for the item 2011-06-02 16:48:35 Last modified in ArcGIS for the item 2011-06-02 16:44:43Automatic updates Have been performed Yes Last update 2011-06-02 16:44:43Item location history  Item copied or moved 2011-06-02 16:48:35 From T:\TestSites\NY\Syracuse\Temp\landcover_2010_syracusecity To \T7500\F$\Export\LandCover_2010_SyracuseCity\landcover_2010_syracusecity

The dataset presented here contains GIS (Geographic Information System) data relevant for infrastructure projects, specifically development of DC grid.

The data has been sourced from a range of sources and proceeded to be compatible for combination and analysis. For each data set included, the conversion of the original data into a 500m grid in EPSG 31468 projection with the raster calculator implemented in QGIS has been carried out.

• Elevation

The elevation data originate from the EEA [1]. The resolution of the original data corresponds to about 25m and is available in raster format. The bilinear resampling method was used to determine the values of the new grid fields.

• Landscape quality assessment

Uniform assessment per grid field between 0 (low) and 10 (high) derived from the original data set of [2] with a total of information on 44 primary land use types based on the assessment criteria of [3]. The data set used here is the CLC2012 which refers to the reference year 2012.

• Population density

Source of population density is a raster grid provided by the [4] which expresses the number of people per pixel with a resolution of 250m x 250m. The base year for population data is 2015.

• Protected areas

Reduction of data from [5] to continental Europe and change of projection.

• Slope

Derived from elevation data based on [1]. The slope is calculated by the tilt angle for each raster cell in degrees based on the first-order derivation.

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).