DOI: 10.6084/m9.figshare.29540903.v1
title: Spatio-Temporal Changes in Habitat Type and Quality in Hong Kong (1973-2022)

description: <b>Supplementary materials used in the following studies:</b>Kwong, I. H. Y., Lai, D. Y. F., Wong, F. K. K., & Fung, T. (Manuscript submitted for publication). Integrating five decades of Landsat imagery for territory-wide habitat mapping and change detection in a subtropical metropolitan city.Kwong, I. H. Y. (2025). <i>Spatio-Temporal Changes in Habitat Type and Quality in Hong Kong Using a 50-Year Archive of Remote Sensing Imagery</i> [Doctoral thesis, Department of Geography and Resource Management, The Chinese University of Hong Kong].Kwong, I. H. Y., Lai, D. Y. F., Wong, F. K. K., & Fung, T. (2025). Spatial variations in forest succession rates revealed from multi-temporal habitat maps using Landsat imagery in subtropical Hong Kong. European Geosciences Union (EGU) General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025. https://doi.org/10.5194/egusphere-egu25-2667. [Poster Presentation: https://presentations.copernicus.org/EGU25/EGU25-2667_presentation-h291057.pdf]<i>Disclaimer: All datasets described here are for reference only. No express or implied warranty or representation is given to the accuracy or completeness of the data or its appropriateness for use in any particular circumstances.</i>GIS mapping results:All raster layers (GeoTiff format) have a pixel size of 30 m covering the 1117-km2 terrestrial area in Hong Kong in this study (Hong Kong 1980 Grid coordinate system). The time period of 1973–2022 was divided into 10 five-year periods in the mapping process.<b>HabitatMapHK_6class_yyyy-yyyy.tif</b>: Raster data showing the 6 habitat classes mapped in this study. Pixel values range from 1 to 6 representing woodland, shrubland, grassland, barren land, built-up area, and water respectively.<b>HabitatMapHK_EstimatedArea.csv</b>: Area coverage (km2) of different habitat classes, as well as their confidence intervals, as mapped in this study.<b>HabitatMapHK_6class_ArcGISsymbology.lyrx</b>: Used to apply the suggested symbology in ArcGIS Pro.<b>ClassificationProbability_yyyy-yyyy.tif</b>: The probability values belonging to each class for every pixel. They were the intermediate products generated from the classification workflow and used to determine the final class with the highest probability and compute the forest index in this study. The sum of probabilities for all six classes is equal to 1. A scale factor of 10000 was applied to the GeoTiff files for storage convenience.<b>HabitatMapHK_8class_yyyy-yyyy.tif</b>: Based on the 6-class outputs, two more classes are added in this product, including wetland (pixel value 7) and plantation (pixel value 8), to serve as inputs for the habitat quality model.<b>HabitatQualityHK_yyyy-yyyy.tif</b>: Habitat quality maps produced in this study. The pixel value is a continuous variable ranging from 0 to 1, with 1 meaning the highest habitat quality.GIS supplementary data:All datasets were collected and compiled from January to June 2024 and represent the conditions at that time.Environmental Raster:<b>DistanceFromCoast.tif</b>: Geometric distance (m) from the coastline.<b>Elevation.tif</b>: Terrain height (m) from a LiDAR-based digital terrain model.<b>Hillfire_10periods.tif</b>: Hill fires occurred in each five-year period, based on burn-area products by Chan et al. (2023) and manual digitisation for early years.<b>Insolation.tif</b>: Annual amount of incoming solar radiation (kWh/m2) computed using SAGA GIS.<b>Landslide_10periods.tif</b>: Landslides occurred in each five-year period, based on the Enhanced Natural Terrain Landslide Inventory (Dias et al., 2009).<b>Northness.tif</b>: Terrain aspect from 1 (due north) to -1 (due south) computed from the DTM.<b>Precipitation.tif</b>: Annual precipitation (mm) (average between 1991-2020) from Hong Kong Observatory.<b>Slope.tif</b>: Steepness (°) of the ground surface computed from the DTM.<b>SoilCEC.tif</b>: Cation exchange capacity (CEC) (mmol/kg) of topsoil from Luo et al. (2007).<b>SoilOrganicMatter.tif</b>: Organic matter content (%) of topsoil from Luo et al. (2007).<b>Temperature.tif</b>: Annual mean temperature (°C) from Morgan and Guénard (2019).<b>TopographicWetnessIndex.tif</b>: Amount of water accumulation due to topographic effects computed using SAGA GIS.<b>Typhoon_10periods.tif</b>: Wind speed (km/h) estimated from WindNinja based on maximum hourly mean wind records associated with typhoon events in each five-year period.<b>WindSpeed.tif</b>: Mean wind speed (km/h) estimated from WindNinja based on monthly prevailing wind records.Human Activities:<b>BuiltupAreas_10periods_shp.zip</b>: Shapefile (polygons) of built-up areas, with attributes on the years of construction (estimated from topographic maps) and density (high and low). It was used as a threat factor in habitat quality mapping and variables in habitat changes.<b>CountryParksProtectedAreas_shp.zip</b>: Shapefile (polygons) of protected areas (Country Parks, Special Areas, etc.), with attributes on the years of designation and revision. It was used as a protection factor in habitat quality mapping and variables in habitat changes.<b>PollutionSource_shp.zip</b>: Shapefile (polygons) of pollution sources (landfills, power stations, and incineration plants), with attributes on the years of construction and closure. It was used as a threat factor in habitat quality mapping.<b>Roads_10periods_shp.zip</b>: Shapefile (polylines) of roads, with attributes on the years of construction (estimated from topographic maps) and type (main and secondary). It was used as a threat factor in habitat quality mapping.Mapping Reference:<b>ForestIndex_FieldCollectedReferenceData.csv</b>: Field survey records of habitat types which were used to evaluate the forest index variable in this study.<b>HabitatMapHK_FieldCollectedReferenceData.csv</b>: Field survey records of habitat types which were used to assess the habitat mapping results in this study.<b>HabitatMapHK_OfficeInterpretedReferenceData.csv</b>: Reference points where the habitat class in each period was determined through visual interpretation of the aerial photographs and other historical records. The points were used for both training and validation of the habitat maps in this study.<b>HabitatQualityHK_FieldSurveyedEcologicalValue2008.csv</b>: Field survey records of ecological values in 2008 which were used to evaluate the habitat quality maps in this study.<b>LandsatHK_CrossSensorCalibrationPoints.csv</b>: Selected points that were assumed to remain unchanged over time and used to cross-calibrate different Landsat sensors in this study.<b>LandsatHK_ImageMetadata.csv</b>: Metadata of the Landsat imagery (1,100 downloaded scenes and 607 valid scenes after pre-processing) acquired and processed in this study.<b>Plantation_1975_1990_2008_2019.tif</b>: Pixels that were identified as plantations on four existing maps in different years (1975, 1990, 2008, 2019), as represented by the four layers contained in this raster file respectively. These pixels were used to help extract plantation class on the habitat map (when producing habitat quality) and denote areas with plantation activities (when modelling habitat changes) in this study.<b>SpeciesObsHK_SpeciesChecklist.csv</b>: A species checklist of 7 taxa in Hong Kong (Plants, Butterflies, Birds, Reptiles, Dragonflies, Amphibians, Mammals) compiled from AFCD, Hong Kong Biodiversity Information Hub, and other secondary sources. Species of conservation concern are identified based on local assessments (Corlett et al., 2000; Fellowes et al., 2002), environmental protection laws, and national and global assessments. The checklist was used to match with the iNaturalist observation data to compute biodiversity metrics at grid levels and evaluate habitat quality maps in this study.<b>SpeciesObsHK_SynonymList.csv</b>: A list of species name synonyms for matching names used in iNaturalist and other secondary sources with the species checklist. It was used to pre-process the iNaturalist observation data and unify the species names from different records in this study.Analysis scripts:Part 1: Mapping Vegetation Habitats from a Satellite Image Time-Series<b>P1_01_SearchAndDownloadFromGEE.ipynb</b>: Query and download all available Landsat 1-9 imagery covering the study area using Google Earth Engine. Atmospheric correction is performed if necessary.<b>P1_02_Preprocess_part1.py</b>: Some basic pre-processing steps after downloading the images from cloud platform to local computer, such as mosaicking adjacent scenes and reprojecting to local coordinate system.<b>P1_03_TopographicCorrection.R</b>: SCS+C topographic correction based on terrain slope, aspect, sun azimuth and sun elevation angles.<b>P1_04_CrossSensorCal.R</b>: Cross-calibration of different Landsat sensors based on pseudo-invariant features, followed by computing variables for image classification.<b>P1_05_ImageComposite.R</b>: Create image composites (median and standard deviation statistics) by combining all imagery acquired in the same period.<b>P1_06_ExtractPixelValue.R</b>: Extract pixel values at the locations of reference points.<b>P1_07_TrainingDataStat.R</b>: Summarise the characteristics of pixel values (e.g., spectral reflectance) of each habitat class and Landsat sensor.<b>P1_08_TrainRFModel.R</b>: Train the Random Forest model, fuse probability outputs from each image, evaluate the model accuracies with cross-validation, and create the final model for classifying the entire dataset.<b>P1_09_TestProcedures.R</b>: Modify the classification procedures and re-run the Random Forest models to evaluate their impacts on the classification accuracies.<b>P1_10_ApplyModel.R</b>: Apply the Random Forest model and fusion steps to all images to create the habitat map for each period.<b>P1_11_AreaCoverage.R</b>: Obtain the area coverage of each class on the habitat map as well as the confidence interval of the area estimates.<b>P1_12_CompareFieldData.R</b>: Assess the accuracies of the habitat maps by overlaying with field-collected points and LiDAR height information at different times.<b>P1_13_SurvivalAnalysis.R</b>: Analyse the number of years required for transitioning between vegetation classes as well as the correlations between transition times and environmental variables.Part 2: Computing Habitat Quality Maps with Reference to Habitat Type Information<b>P2_01_ModifyInputMap.R</b>: Add two more classes (wetland and plantation) to the habitat map to serve as inputs for the habitat quality model.<b>P2_02_SensitivityParameter.R</b>: Compute the sensitivity parameter based on the statistical relationship between habitat changes and nearby threat factors to serve as inputs for the habitat quality model.<b>P2_03_RunInVESTModel.R</b>: Execute the InVEST habitat quality model with both default and optimal half saturation constant values.<b>P2_04_CompareFieldSurvey.R</b>: Evaluate the habitat quality mapping results with field-surveyed ecological value data.<b>P2_05_CompareSpeciesObservation.R</b>: Pre-process and match the iNaturalist observation data, compute biodiversity metrics at grid levels, and compare with the estimated habitat quality values.<b>P2_06_DegradationRarity.R</b>: Create degradation and rarity outputs from the InVEST habitat quality model and summarise their patterns over time.Part 3: Modelling Habitat Changes in Relation to Natural and Human Factors<b>P3_01_CreateEnvRaster.R</b>: Create raster layers representing the environmental variables to serve as model inputs in this part of the study.<b>P3_02_ForestIndex.R</b>: Create forest index raster by combining the probability layers generated from the Random Forest model and evaluate it with habitat classes and field-surveyed points.<b>P3_03_ExtractPixelValue.R</b>: Extract the pixel values of the environmental variables for each of the 10 periods at the locations of randomly selected points.<b>P3_04_EFA.R</b>: Conduct exploratory factor analysis for the variables related to the natural environment and reduce them to three factors.<b>P3_05_SEM.R</b>: Conduct structural equation modelling based on the hypothesised model and computed variables.<b>P3_06_ModelValidation.R</b>: Test the assumptions of individual paths in the structural equation model and compute metrics with cross-validation.<b>P3_07_YearChange.R</b>: Examine the model interaction with time and path coefficients that vary against time.<i>Last updated in September 2025</i>