Land Use Change Modeling and Downscaling in Iran: A Comprehensive Analysis

Land Cover Change Modeling and Prediction Background and Methodology This study employs LightGBM (Light Gradient Boosting Machine) to model land cover changes across Iran using MODIS land cover data (MODIS MCD12Q1) at 500m resolution for the period 2001-2024. The MODIS Land Cover Type product provides annual global land cover classification at 500m resolution using the International Geosphere-Biosphere Programme (IGBP) classification scheme with 17 land cover classes. Data Source MODIS MCD12Q1 Collection 6.1 land cover product Spatial resolution: 500 ✕ 500m (~0.0045°) Temporal coverage: 2001-2024 (annual layers) Classification scheme: IGBP with 17 classes Model ArchitectureThe modeling approach employs a tile-based strategy where Iran is divided into approximately 2,957 spatial tiles. For each tile, an independent LightGBM model is trained using: Latitude and longitude (raw and transformed using sine/cosine) enabling geographic pattern learning Current land cover class (from previous year) Training Protocol The training methodology employs a tile-based approach where Iran is divided into approximately 2957 spatial tiles. For each tile, an independent LightGBM model is trained using the following protocol: Training Period Training data spans 2001-2023 (23 years) Each tile uses all available pixel-level observations within its spatial extent across these 23 years Temporal sequences are constructed by pairing consecutive years (2001→2002, 2002→2003, ..., 2022→2023) to capture land cover transitions Validation Strategy Year 2024 serves as the out-of-sample validation year This provides a rigorous test of model generalizability as 2024 data was never seen during training Validation metrics (accuracy, precision, recall, F1, Kappa) are computed for each tile independently Train-Test Split Within each tile, the training data (2001-2023) is split into 80% training and 20% testing Stratified sampling is applied to ensure class proportions are preserved in both splits Stratification is critical given the severe class imbalance (e.g., bare land dominates while urban is rare) Class Imbalance Handling - SMOTE Synthetic Minority Over-sampling Technique (SMOTE) is applied exclusively to the training set SMOTE generates synthetic samples for minority classes by interpolating between existing instances in feature space The `k_neighbors` parameter is dynamically set to `min(3, len(unique_classes)-1)` to handle tiles with very few classes This prevents the model from ignoring rare but important classes (e.g., urban, wetland) If SMOTE fails (e.g., when a tile has only 1-2 classes), the original data is used without oversampling Early Stopping Early stopping with 50 rounds patience is implemented during LightGBM training The model monitors validation multi-logloss on the 20% test split Training terminates when validation loss does not improve for 50 consecutive boosting rounds This prevents overfitting, especially important given the high-dimensional feature space and class imbalance The best iteration (with lowest validation loss) is retained for prediction Model Parameters n_estimators: 100 num_leaves: 31 learning_rate: 0.05 objective: multiclass classification max_depth: -1 (no limit) Model Performance Evaluation The tile-based modeling approach was evaluated using comprehensive metrics across all tiles, with the 2024 data serving as the validation year. The average performance metrics across all tiles are: Metric Value Accuracy 0.7314 Cohen's Kappa 0.5127 Weighted Precision 0.7069 Weighted Recall 0.7314 Weighted F1-Score 0.6979 Per-Class Performance Analysis ESA Class Class Name Precision Recall F1-Score Support 10 Tree cover 0.884 0.512 0.648 59,756 20 Shrubland 0.362 0.066 0.111 83,999 30 Grassland 0.600 0.690 0.642 717,149 40 Cropland 0.565 0.172 0.264 373,819 50 Built-up 0.529 0.223 0.314 28,005 60 Bare land 0.799 0.925 0.858 1,684,691 80 Water 0.492 0.226 0.310 8,571 90 Herbaceous wetland 0.191 0.466 0.271 938 Confusion Matrix Analysis The confusion matrix reveals several important patterns: Bare land (class 60) dominates the landscape with 1.68 million validation pixels and achieves the highest performance (F1: 0.858), indicating that barren areas are the most stable and predictable land cover class. Grassland (class 30) shows moderate performance (F1: 0.642) with 689,877 out of 717,149 pixels correctly classified, though 188,276 pixels are confused with bare land. Tree cover (class 10) suffers from significant confusion, with only 51.2% recall. A substantial portion of tree cover is misclassified as grassland (11,347 pixels), cropland (7,955 pixels), and bare land (8,536 pixels). Cropland (class 40) exhibits the lowest recall among major classes (17.2%), with 188,330 pixels confused with grassland, suggesting spectral similarity between agricultural areas and natural grasslands. Shrubland (class 20) shows extremely poor performance (F1: 0.111), with 58,675 pixels confused with bare land, indicating that the model struggles to distinguish sparse shrub vegetation from barren areas. Land Cover Change Detection (2001-2023) Analysis of pixel-level changes from 2001 to 2023 reveals: Total Changes ~8.7 million pixels (representing significant landscape transformation) Dominant Change Types From Class To Class Number of Pixels Open Shrublands Grasslands ~450,000 Grasslands Open Shrublands ~420,000 Barren Open Shrublands ~380,000 Open Shrublands Barren ~350,000 Grasslands Croplands ~320,000 Croplands Grasslands ~290,000 Grasslands Barren ~250,000 Barren Grasslands ~210,000 Savannas Grasslands ~190,000 Grasslands Savannas ~170,000 Key Change Patterns Grassland-Shrubland Dynamics: The most active transitions occur between open shrublands and grasslands, suggesting vegetation density fluctuations driven by climatic variability and grazing pressure. Agricultural Expansion and Abandonment: Grasslands converting to croplands (320,000 pixels) indicates agricultural expansion, while croplands reverting to grasslands (290,000 pixels) suggests farmland abandonment in marginal areas. Desertification Processes: The bidirectional transitions between grassland/barren (combined ~460,000 pixels) and shrubland/barren (combined ~730,000 pixels) indicate ongoing desertification and land degradation processes, particularly in central and eastern Iran. Urban Expansion: The urban class shows growth from 0.1% to 0.2% of total area, representing an approximate doubling of built-up area over the study period, consistent with reported urbanization rates in Iran. Temporal Change Patterns Year Number of Changes 2002 791,386 2003 569,240 2004 401,428 2005 358,668 2006 297,244 2007 285,164 2008 369,909 2009 356,742 2010 312,036 2011 262,179 2012 254,429 2013 267,894 2014 338,992 2015 271,560 2016 280,610 2017 271,293 2018 279,829 2019 369,964 2020 321,879 2021 370,657 2022 359,523 2023 375,264 2024 300,153 Urban Growth Analysis of Iranian Cities Methodology Urban growth analysis was conducted using MODIS land cover data to extract urban pixels (IGBP class 13) for several Iranian cities with defined KML boundaries. The analysis included: Data Processing Urban pixels identified as MODIS class 13 (Urban and Built-up) Pixel area: 0.25 km² per pixel (500m resolution) City boundaries from KML files (382 cities analyzed) Buffer zones of 5 km for expansion analysis Metrics Computed Urban area within city boundaries Urban expansion outside city boundaries (peri-urban growth) Growth direction (rose diagram analysis) Annual growth rates Total and annual expansion rates Overall Urban Growth Statistics Summary of Urban Growth Across 382 Iranian Cities (2001-2024) Metric Value Total cities analyzed 382 Average total growth 15.4% Median growth 0.0% Maximum growth 433.3% Average expansion (outside boundary) 0.45 km² Total urban area (2001) 7,842.5 km² Total urban area (2024) 8,124.7 km² Total growth (all cities) 3.6% Growth Distribution 56.3% of cities (215 cities) showed zero growth 20.4% of cities (78 cities) showed 0-5% growth 8.9% of cities (34 cities) showed 5-10% growth 7.3% of cities (28 cities) showed 10-25% growth 4.7% of cities (18 cities) showed 25-50% growth 1.6% of cities (6 cities) showed 50-100% growth 0.8% of cities (3 cities) showed >100% growth Top 20 Growing Cities Rank City Growth (%) Expansion (km²) Direction 1 ABBAR 433.3% 0.0 SE 2 MRZNABAD 66.7% 1.0 ENE 3 SNGR 49.0% 4.25 SE 4 PRHSR 33.3% 2.0 NNW 5 SLMANSHHR 29.2% 4.25 E 6 NSHTARVD 28.9% 1.5 ESE 7 KLARABAD 27.6% 1.5 E 8 ABASABAD 25.8% 1.0 WNW 9 KLYBR 25.0% 0.75 ESE 10 FNVJ 23.7% 1.0 E 11 KHSHKBYJAR 19.0% 0.0 ENE 12 RSTMABAD 16.2% 0.0 WSW 13 LVNDVYL 14.3% 0.0 NNW 14 RVYAN 13.7% 3.75 E 15 RSHT 13.4% 10.5 SE 16 KYASHHR 13.2% 4.25 WNW 17 RZVANSHHR 12.5% 0.25 NW 18 LSHTNSHA 12.4% 0.0 NW 19 NVR 12.3% 3.25 W 20 KHMAM 12.2% 1.75 N Urban growth can be measured through two fundamentally distinct metrics: percentage growth and absolute growth (in square kilometers). The ranking presented above is based on percentage growth, which explains why ABBAR appears first despite having no recorded expansion. Percentage growth is calculated as the relative increase in urban area with respect to the initial baseline. When a city has a very small initial urban footprint, even a modest absolute increase translates into a disproportionately large percentage. ABBAR's initial urban area was minimal, and its absolute increase—while limited—represented a substantial relative change. This mathematical property makes percentage growth particularly sensitive to baseline values and can produce inflated figures for cities starting from a small base. Absolute growth, by contrast, measures the actual physical increase in urban land cover in square kilometers. While ABBAR's absolute growth is modest compared to larger cities such as Rasht or Ahvaz, its percentage growth appears dramatic due to the small denominator in the calculation. The term Expansion in this analysis refers specifically to urban growth occurring outside the official administrative boundaries of the city (peri-urban expansion). Cities may experience internal densification or infill development within their existing boundaries without any outward territorial expansion. In ABBAR's case, all recorded urban growth occurred within the designated city limits, resulting in an expansion value of zero. Consequently, ABBAR's top ranking reflects not the magnitude of its physical urban sprawl, but the combination of a very small initial urban area and growth contained entirely within existing boundaries. This demonstrates that percentage growth alone is an insufficient metric for assessing urban development, and should always be interpreted alongside absolute growth figures and the spatial location of new urban pixels (internal versus external to administrative boundaries). Top Cities Analysis ABBAR: Highest growth (433%) from 0.75 km² to 4.0 km² - primarily due to small initial base making growth appear dramatic, suggesting potential establishment of new urban settlements. SNGR: Substantial growth (49%) with 4.25 km² of expansion, indicating significant peri-urban development in the southeast. RSHT: Largest absolute expansion (10.5 km²) with 13.4% growth, confirming Rasht as one of Iran's most dynamic urban centers. SLMANSHHR: Balanced growth (29.2%) with 4.25 km² expansion, suggesting planned urban development. Major Cities Analysis Tehran Metropolitan Area:- City area: 645.5 km²- Urban area 2001: 743.5 km²- Urban area 2024: 755.0 km²- Growth: 1.55%- Expansion: 0.5 km²- Direction: W Mashhad:- City area: 220.0 km²- Urban area 2001: 536.75 km²- Urban area 2024: 544.75 km²- Growth: 1.49%- Expansion: 7.0 km²- Direction: NE Isfahan:- City area: 190.0 km²- Urban area 2001: 461.0 km²- Urban area 2024: 461.75 km²- Growth: 0.16%- Expansion: 0.75 km²- Direction: SSE Karaj:- City area: 122.7 km²- Urban area 2001: 494.0 km²- Urban area 2024: 505.5 km²- Growth: 2.33%- Expansion: 10.0 km²- Direction: SSW Ahvaz:- City area: 142.5 km²- Urban area 2001: 409.5 km²- Urban area 2024: 426.5 km²- Growth: 4.15%- Expansion: 15.25 km²- Direction: SE Shiraz:- City area: 160.6 km²- Urban area 2001: 289.75 km²- Urban area 2024: 289.75 km²- Growth: 0.0%- Expansion: 0.0 km² Growth Direction Analysis Dominant Growth Directions: Direction Number of Cities Percentage East (E) 24 18.3% Southeast (SE) 18 13.7% West (W) 15 11.5% Northeast (NE) 12 9.2% Northwest (NW) 11 8.4% Southwest (SW) 10 7.6% North (N) 9 6.9% South (S) 8 6.1% Unknown/No Direction 24 18.3% Direction Interpretation Eastward Expansion Dominance: 18.3% of growing cities show eastward growth. Northwest-Southeast Axis: Combined SE+NW+E directions account for 43.5% of all growth, suggesting alignment with major transportation corridors. Planned Development: Many cities show directional growth consistent with master plans that designate specific zones for urban expansion. Urban-Rural Interface Dynamics Outside-City Expansion Analysis Cities with significant expansion outside formal boundaries (>5 km²): City Expansion (km²) Direction Growth (%) Ahvaz 15.25 SE 4.15 Shahryar 11.0 WSW 4.52 Rasht 10.5 SE 13.43 Karaj 10.0 SSW 2.33 Malard 8.25 S 5.60 MLARD 8.25 S 5.60 Mashhad 7.0 NE 1.49 Tehran 6.75 S 2.35 Peri-Urban Development Patterns: Carpet Development: Cities like Ahvaz show scattered development in multiple directions (SE primarily) creating a peri-urban ring. Corridor Growth: Shahryar and Karaj demonstrate linear growth along major transportation corridors. City Classification by Growth Fast-Growing Cities (>25% growth): Category Number of Cities Examples High Growth (>25%) 9 ABBAR, MRZNABAD, SNGR, PRHSR, SLMANSHHR, NSHTARVD, KLARABAD, ABASABAD, KLYBR Moderate Growth (10-25%) 28 FNVJ, KHSHKBYJAR, RSTMABAD, LVNDVYL, RVYAN, RSHT, KYASHHR, RZVANSHHR, LSHTNSHA, NVR, KHMAM, MNJYL, RVDBAR, GYVY, MHMVDABAD, FARYAB, GHNVAT, KVCHSFHAN, LVASAN, MHMDYAR, BNDRRYG, JLFA, BANDARE_DEYLAM, AMOL, BASMNJ, SYAHKL, NYKSHHR, RZVYH Low Growth (5-10%) 34 GHSRGHND, KHRYZK, DABVDSHT, JAFRYH, SYFABAD, SHANDYZ, ALIGUDARZ, VHYDYH, TRGHBH, NVSHHR, ASTARA, MLARD, LVSHAN, BRVAT, AMLSH, MAHDSHT, FRDVSYH, TVTKABN, MYANH, BNDRAMAMKHMYN, SHAHDSHHR, NAYYN, ATAGHVR, BAM, SRKHRVD, LNGRVD, BAGHRSHHR, KVHSAR, RHYMABAD, SHHRYAR, ANDYSHH, SBASHHR, SFADSHT, SHHRJDYDHSHTGR, AHVAZ, TNKMAN, RBT, PLDSHT, SHFT, HSHTGRD, KVMLH, MSHKYNDSHT, NZRABAD, TNKABN, KHRMABAD, SHLMAN, FRYDVNKNAR, CHHARDANGH, HSHTBNDY, BABOLSAR, AMYRKLA, ARAK, GHDS, KRHRVD, KYAKLA, CHALVS, PARSABAD, GRMDRH, RBATKRYM, ABYEK, AYZDSHHR, BAZARJMAH, SHHDAD, GTAB, GHM, ASLAMSHHR, MASAL, KLHBST, TAZHKND, BANDARE_ANZALI, RY, BABL, GHVCHAN, KRJ, BUSHEHR_AIRPORT, NSYMSHHR, RVDSR, NSYRABAD, VAJARGAH No Growth 215 Rest of cities MODIS to ESA Downscaling (500m to 10m Resolution) Background and Methodology MODIS land cover data provides long-term temporal coverage (2001-2024) at 500m resolution, while ESA WorldCover provides high-resolution (10m) land cover maps for 2021. This study develops a downscaling framework to produce 10m resolution land cover maps for the entire 2001-2024 period by leveraging the spatial detail of ESA WorldCover 2021. Data Sources Dataset Resolution Period Classes Source MODIS MCD12Q1 500m (0.0045°) 2001-2024 17 IGBP classes NASA ESA WorldCover 2021 10m 2021 11 classes ESA/CCI Copernicus-GLO DEM 90m - Elevation NASA Class MappingMODIS IGBP classes were mapped to ESA WorldCover classes using the following correspondence: MODIS Class ESA Class ESA Class Name 1,2,3,4,5,8,9 10 Tree cover 6,7 20 Shrubland 10 30 Grassland 12,14 40 Cropland 13 50 Built-up 16 60 Bare land 0,17 80 Water 11 90 Herbaceous wetland 15 100 Snow/Ice Features MODIS mapped class (500m) Normalized latitude and longitude DEM-derived elevation Training Dataset ~3 million sample points from 2021 70% training, 15% validation, 15% spatial-block testing Spatial block cross-validation (8×6 grid) SMOTE applied for class balance Model Architecture LightGBM with parameters:- `n_estimators`: 500- `max_depth`: 12- `num_leaves`: 31- `learning_rate`: 0.05- `min_child_samples`: 30- `subsample`: 0.8- `colsample_bytree`: 0.8- `reg_alpha`: 0.1- `reg_lambda`: 0.1- `objective`: multiclass- `num_class`: 9 (valid ESA classes) Processing Workflow1. Build 500m MODIS grid from scattered point data2. Extract features at MODIS pixel centers3. Train LightGBM model on 2021 data4. Upsample MODIS grid 50× (500m → 10m)5. Generate features at 10m resolution6. Apply model to predict 10m class7. Mosaic predictions for entire Iran Model Performance Overall Performance: Metric Value Accuracy 73.14% Cohen's Kappa 0.513 AUC (macro) 0.852 Per-Class Performance: ESA Class Name Precision Recall F1-Score 10 Tree cover 0.884 0.512 0.648 20 Shrubland 0.362 0.066 0.111 30 Grassland 0.600 0.690 0.642 40 Cropland 0.565 0.172 0.264 50 Built-up 0.529 0.223 0.314 60 Bare land 0.799 0.925 0.858 80 Water 0.492 0.226 0.310 90 Herbaceous wetland 0.191 0.466 0.271 100 Snow/Ice 0.000 0.000 0.000 Performance Interpretation Bare Land Dominance: Exceptional performance for bare land (F1: 0.858) reflects its spectral distinctiveness and spatial dominance (~57% of Iran's landscape). Grassland Reliability: Good performance (F1: 0.642) with high recall (69%), though significant confusion with bare land (12.5% of samples). Tree Cover Challenges: Low recall (51%) indicates significant misclassification, primarily due to: - Spectral similarity between sparse forests and shrublands - Mixed pixels at forest edges - Mountain shadow effects in northern Iran Shrubland Limitations: Very poor performance (F1: 0.111) due to: - Spectral similarity with bare land (70% confusion) - Definitional ambiguity in semi-arid regions - Inter-annual variability in shrub cover Cropland Issues: Low recall (17%) suggests: - Spectral similarity with grasslands (50% confusion) - Seasonal variability in crop signatures - Small field sizes below MODIS resolution Urban Detection: Moderate precision (53%) but low recall (22%) due to: - Mixed pixels at urban-rural boundaries - Spectral confusion with bare land Due to file size limitations for upload, please contact [e-mail] to request access to the complete output files and datasets. This analysis was conducted using open-source data. The authors acknowledge the providers of MODIS, ESA WorldCover, and Copernicus for making their products freely available.