We introduce landlensdb, an open-source Python package for managing proximity sensing imagery, including action cameras, 360° cameras, and UAVs, using PostgreSQL/PostGIS. It automates metadata extraction, corrects geolocation errors via road network snapping, and enables scalable spatial-temporal queries and visualization for large-scale geotagged image datasets.

Biological invasions represent one of the most significant threats to global biodiversity and agricultural systems, causing substantial ecological and economic damage worldwide. Among emerging invasive pests in East Asia, Aromia bungii, commonly known as the red-necked longhorn beetle, has become a serious threat in Japan. The species attacks several Prunus species, including ornamental cherry trees (Cerasus spp.), peach (Prunus persica), and plum (Prunus salicina). Because cherry trees play an important ecological and cultural role in Japan, the spread of this invasive beetle has raised growing concerns for landscape management and biodiversity conservation.
Since its first detection in Aichi Prefecture in 2012, A. bungii has expanded rapidly across urban and peri-urban areas. Understanding its spatial pattern is therefore essential for effective monitoring and early intervention. Spatial analysis can be applied to these processes. However, such analyses fundamentally depend on how spatial data are defined, including the geometry of the spatial grid, which can influence the results.
Thus, to examine how grid shape influences spatial analysis results, this study evaluates spatial autocorrelation measures using different tessellations of invasive species occurrence data. Specifically, we compared rectangular and hexagonal grids for analysing spatial patterns in A. bungii occurrence records and density of rivers in Saitama Prefecture, Japan.

Volunteer-based occurrence data for A. bungii were compiled from field surveys across Saitama Prefecture from 2017 to 2024, yielding 2,349 confirmed presence records. Records were classified as confirmed presences if either adult beetle observations (Adult-yes = 1) or evidence of tree damage (Tree_damage = 1) was recorded. All records were georeferenced using latitude–longitude coordinates (WGS84) and then reprojected to UTM Zone 54N (EPSG:32654) to ensure metric accuracy in spatial calculations. We assembled the occurrence data into predefined grid cells to explore spatial patterns across the study area and to enable consistent spatial aggregation and neighbourhood-based analyses.

Two grid tessellation schemes were constructed over the study area. The rectangular grid consisted of 4,013 cells at a 1 km × 1 km resolution. Each rectangular cell overlaid with the occurence point of location Aromia bungii presence in Saitama. This configuration corresponds to rook contiguity, where only four directly adjacent neighbours are considered. To evaluate the effect of diagonal bias, the same rectangular grid was also analysed using queen contiguity, which includes eight neighbouring cells by incorporating both direct and diagonal neighbours. The hexagonal grid was generated using Shapely and Geopandas with hierarchical spatial indexing system, producing hexagonal cells with an almost equivalent spatial resolution to the rectangular grid. To keep the overall grid coverage comparable to the rectangular representation, 4,007 hexagonal cells were generated to cover the study area. The rectangular grid had a side length of 1.0 km, whereas the hexagonal grid had a slightly larger side length of 1.074 km. Each hexagonal cell has a circumradius of 620m and a cell area 0.999km square approximately equivalent to 1km square rectangular grid ensuring comparable spatial resolution between two tessellations. From these two tessellation we set four spatial weight configurations which are; Rook, standard Queen, Distance-Weighted Queen and Hexagonal KNN-6 and then calculated the spatial autocorrelation.

We found Global spatial autocorrelation analysis confirmed significant and positive clustering of A. bungii occurrences across all configurations Moran's I = 0.4525--0.5380, Geary's C = 0.5527--0.5970, Getis-Ord G all p-value = 0.001, demonstrating that the spatial agglomeration characteristic is robust to the choice of tessellation geometry and neighbourhood definition. Local spatial autocorrelation analyses consistently identified two primary hotspot zones at the northern (Kazo--Gyoda) and the southeastern (Soka) parts of Saitama Prefecture across all configurations, providing spatially explicit evidence of persistent infestation core that may serve as priority zones for targeted surveillance and countermeasure deployment.　Across all configurations hotspot areas were consistent, however coldspot delineation was found to be unstable with certain configurations inconsistent coldspot that undermine the accurate identification of management priority zones.

Hexagonal tessellation produced higher counts of similar cells associated with hotspot clusters compared to rectangular configurations due to six- neighbourhood structure equidistant by design, ensuring the spatial relationships are evaluated uniformly in all directions thereby reducing directional bias. In contrast, rectangular configurations evaluate neighbours at unequal distances, introducing directional bias that result in fragmented cluster boundaries and lower hotspot cell counts. The inconsistency of coldspot indicating the spatial clustering outcomes are sensitive to tessellation choice. Practitioners should therefore consider the tessellation configuration carefully when interpreting spatial cluster results for management planning, as different configurations may lead to different prioritisation of monitoring and eradication efforts, specifically at cold spot areas.

The limitation of different geometry structure should be acknowledged. Although, the comparable spatial resolution were used, the exact areal equivalence between rectangular and hexagonal tessellation could not be perfectly achieved due to their difference geometric structure.

Accurate land-cover classification maps are essential geospatial resources that support a wide range of societal applications, including urban development, environmental monitoring, and disaster risk management. With the growing availability of satellite remote sensing data, pixel-wise classification methods, which assign a single land-cover class to each image pixel, have become widely adopted. However, when medium-resolution satellite imagery, such as Sentinel-2, is used, a single pixel often contains more than one land-cover type. This phenomenon, commonly referred to as the mixed-pixel problem, leads to the loss of information about minority classes within a pixel, making it difficult to accurately characterize fine-scale surface conditions and limiting the practical utility of classification results. One promising solution to this issue is compositional classification, which estimates the proportion of each land-cover class present within a pixel rather than assigning a single label. This approach retains information about all classes, including those that occupy only a small fraction of a pixel. Compositional data are subject to two mathematical constraints: all values must be non-negative, and the proportions across all classes must sum to one. Yet, the implementation of such compositional classification is challenging, as machine learning or deep learning model architectures do not consider these characteristics properly. This study proposes a deep learning framework for compositional land-cover estimation at 10 m spatial resolution. We explore combinations of input features, model architectures, and loss functions suited to compositional outputs. Rather than seeking a single definitive solution, we aim to provide practical insights into how these design choices interact and influence estimation quality.

As an experiment, we examined two inputs, three models, and two loss functions. Two types of 10-m level input features were examined. The first was spectral reflectance data from 10 multispectral bands (B2-8, B8A, B11, and B12) in Sentinel-2, which directly capture surface characteristics across visible and infrared wavelengths. The second was Embedding V1, a 64-dimensional feature vector generated by AlphaEarth Foundations, a geospatial embedding model that integrates multiple data streams, including optical, radar, and LiDAR observations from multi-temporal satellite imagery. We used OpenEarthMap (OEM), a publicly available global dataset with high-resolution eight land-cover labels at 0.25–0.5 m spatial resolution, as reference data. To generate training targets at medium resolution, we aggregated all OEM pixels within each 10 m pixel and computed the fractional coverage of each land-cover class. This aggregation produced pixel-level composition vectors that served as ground truth for supervised learning, establishing a direct correspondence between medium-resolution inputs and high-resolution reference proportions. Three deep-learning model architectures were tested: a multilayer perceptron (MLP), a two-dimensional convolutional neural network (2D-CNN), and a three-dimensional convolutional neural network (3D-CNN). All models consisted of two fully connected layers followed by a Softmax activation function to ensure outputs formed valid compositional proportions. The CNN-based models were included to capture spatial context from neighboring pixels, while the MLP was used as a simpler baseline that processes each pixel independently. The dataset was split into 70% for training, 10% for validation, and 20% for testing. We compared two loss functions: Mean Absolute Error (MAE), a widely used metric in regression tasks, and the Aitchison distance, a measure derived from compositional data analysis. The Aitchison distance evaluates the relative differences between components rather than their absolute errors by applying a logarithmic transformation and centering to each proportion. To eliminate redundancy associated with the compositional constraint, we applied a linear transformation using an orthonormal basis, projecting the data into a D−1 dimensional space. This metric enables stable distance computation while reflecting the geometric structure and constraints inherent to compositional data. This combination of architectures and loss functions allowed for a systematic comparison across multiple design dimensions. We evaluated the performance of each combination using MAE and the Aitchison distance to identify the most effective framework for compositional land-cover classification mapping. All processing and training are implemented with open-source Python geospatial/ML tools, and we will release code and experiment configurations to support reproducibility.

The experimental results revealed clear differences based on the combination of inputs and model architecture. When using only Sentinel-2 spectral bands as input, 3D-CNN models outperformed the MLP (MAE of 0.1126 versus 0.1185), indicating that spatial context from convolutional operations benefits models when input features are limited to raw reflectance. In contrast, when using Embedding V1 as input, the MLP achieved the best performance (MAE of 0.0989) among all tested models. This suggests that Embedding V1 already encodes rich spatial and semantic information through its pretraining process, making simpler models better suited to leverage these representations without overfitting to their internal structure. Regarding the choice of loss function, the resulting maps showed a meaningful difference between MAE and Aitchison distance. Models trained with MAE tended to produce smoothed predictions, with estimated proportions concentrated near intermediate values and rarely approaching 0 or 1. This smoothing effect is a known limitation of MAE in compositional contexts, where extreme values are physically meaningful. When the Aitchison distance was used instead, the estimated proportion maps showed sharper contrasts and more realistic spatial patterns that better matched the actual distribution of land-cover classes. These results suggest that using a loss function aligned with the mathematical properties of compositional data leads to more credible and interpretable outputs.

The findings demonstrate that the combination of input features, model architecture, and loss function substantially impacts the quality of fractional land-cover estimation. Under our experimental conditions, the MLP model using Embedding V1 and the Aitchison distance achieved the highest accuracy, providing practical insights into the design of compositional land cover classification frameworks in remote sensing.

Floods are becoming more frequent and severe worldwide due to climate change. These disasters cause significant human and economic losses, making the development of rapid and highly accurate flood mapping technologies essential for clearly identifying damage conditions and supporting subsequent rescue activities.

Remote sensing with openly available Synthetic Aperture Radar (SAR) plays an essential role in disaster monitoring. As an active sensor that transmits its own microwave signals, SAR can acquire observations regardless of weather or time of day, enabling reliable monitoring even in extreme conditions.

However, conventional SAR-based flood detection methods face several challenges. Threshold-based approaches applied to a single SAR image struggle to distinguish permanent water bodies from newly inundated areas. To address this, methods that compare pre-flood and flood-time images have been developed. Yet openly available data such as Sentinel-1 have satellite revisit cycles that create temporal gaps of several days to weeks between compared images. During this interval, seasonal vegetation changes and speckle noise differences are frequently misidentified as flood-induced changes, leading to increased false positives. To address these false positives caused by temporal misalignment, this study proposes a flood detection method that synchronizes the temporal axes of the compared images.

In our experiment, we used the openly available "UrbanSARFloods" dataset. This dataset was created for global-scale flood detection and covers 18 flood events from five continents, including flood and non-flood label data. Each image chip has a spatial resolution of 10 m and a size of 512 × 512 pixels, and data from a selected 628 area of interest were used in this study.

We used Sentinel-1 Ground Range Detected (GRD) data as our primary flood detection source. Sentinel-1 GRD is a multi-looked SAR product projected to ground range using an Earth ellipsoid model. In Interferometric Wide (IW) mode, it provides roughly 10-meter resolution with dual-polarization (VV and VH) across a 250 km swath. Sentinel-1 acquires observations regardless of weather or time of day, with a 6–12 day revisit cycle. We collected SAR time series from Google Earth Engine (GEE), covering one year prior to each flood event up to 40 images per location. These data include two polarization modes: VV and VH. To ensure consistency, only data acquired from the same orbit (ascending or descending) as the flood-time observation were used.

GRD data contain inherent speckle noise that requires filtering. To improve data quality, we reduced speckle noise using a median filter (3 × 3) and performed normalization using Min–Max scaling. For datasets with missing time-series images, we interpolated missing data using three-dimensional spline interpolation.

The revisit gap can cause temporal misalignment between compared images. To obtain an image representing non-flood conditions, we simulated imagery from past SAR GRD sequences at the timing of the flood observation. A three-layer ConvLSTM model captured spatio-temporal dependencies and generated predicted "non-flood" SAR images. The model used a hybrid Mean Absolute Error (MAE)-Structural Similarity (SSIM) loss function to preserve edge details and suppress speckle noise. This function prioritizes structural and statistical similarity over standard MSE.

We created multiple image chip pairs by using six consecutive time steps as one set and sliding the window forward by one time step. Each image chip pair consists of five continuous input images and one subsequent image as reference. The model training was performed by minimizing the error between the predicted image generated from the five input images and the reference image. A dedicated ConvLSTM model was built for each location to learn its unique surface characteristics. After training, we input the five most recent image chips into the model to generate a predicted image for the next time step.

Finally, we performed flood detection using a Siamese Network to compare the simulated images with the observed images acquired during flooding. This CNN-based model uses an encoder–decoder architecture with skip connections, extracting high-dimensional features while preserving spatial boundary information. After computing pixel-wise differences between the feature maps, we applied a sigmoid function to produce a flood probability map. The final flooded areas were detected through thresholding (0.5).

Quantitative evaluation using a pre-split test dataset consisting of 30 samples demonstrated considerable performance improvements. The proposed method achieved an F1 score of 0.605 (Precision: 0.556, Recall: 0.661), representing a 42% improvement over the conventional pre-/post-flood comparison method, which achieved only 0.427 (Precision: 0.328, Recall: 0.610). This improvement was primarily attributed to the substantial reduction in false positives caused by temporal misalignment.

The precision improvement from 0.328 to 0.556 indicates that the proposed temporal synchronization approach effectively suppresses misdetections caused by seasonal vegetation changes and speckle noise variations. Meanwhile, the recall score of 0.661 demonstrates the method's capability to detect actual flood events without sacrificing sensitivity.

Qualitative evaluation further confirmed the effectiveness of the proposed approach. Visual inspection revealed that the method successfully suppresses speckle-like noise patterns that frequently lead to false positives in conventional methods. The approach accurately detected both large-scale inundated areas, such as flooded river plains and urban zones, as well as small, isolated flooded regions that are typically difficult to identify. The generated flood probability maps showed clear spatial boundaries between flooded and non-flooded areas, enabling more reliable damage assessment.

While we demonstrated the usefulness of the flood detection method, several challenges were identified. In regions with complex topography, frequent misdetections occurred due to terrain-induced radar shadow, geometric distortions (layover), and changes in vegetation's dielectric properties. In urban areas, flood detection often misses actual flooding because double-bounce scattering reduces the backscatter intensity differences between images. Future work will consider incorporating Digital Elevation Model (DEM) data and coherence data, which quantify the phase correlation of microwave signals between two temporal images, to account for physical terrain constraints and further suppress misdetections in urban areas, enabling more effective flood detection.

In conclusion, this study proposed a flood detection method that eliminates temporal lag between compared datasets and demonstrated improved performance in detecting flood-induced changes. By enabling near-real-time flood detection under extreme weather conditions, this method is expected to strengthen societal response capabilities to flood disasters.

Forest vertical structural diversity is a key indicator of ecosystem complexity, habitat heterogeneity, and biodiversity potential. Foliage height diversity (FHD), derived from vertical vegetation profiles, is widely used to quantify this structural heterogeneity. Airborne laser scanning (ALS) provides accurate three-dimensional forest structure information, but its limited spatial coverage and high cost hinder large-scale monitoring. The Global Ecosystem Dynamics Investigation (GEDI) mission has enabled global sampling of forest vertical structure using spaceborne LiDAR. However, its footprint-based sampling produces spatially discontinuous observations. As a result, continuous regional-scale mapping of forest vertical structural diversity remains a major challenge. Previous studies have focused primarily on height-related metrics or canopy cover estimation. Regional-scale mapping of foliage height diversity remains limited, especially in cool-temperate and boreal forest ecosystems like northern Japan, where complex terrain, climate gradients, and forest management regimes interact.
This study proposes a large-scale framework for estimating forest vertical structural diversity across Hokkaido, Japan, by integrating GEDI-derived FHD with multi-source satellite remote sensing data and machine learning. The novelty of this work lies in two key contributions: (i) integrating a wide range of complementary satellite data sources including multi-frequency SAR, optical imagery, climate variables, land cover products, topography, and nighttime light data, for FHD estimation, and (ii) explicitly interpreting model behavior using Shapley additive explanations (SHAP) to identify physically meaningful drivers of forest vertical structural diversity.
We used diverse openly available geospatial data. GEDI Level 2B observations acquired in 2023 served as reference data. After quality filtering, GEDI footprints were spatially matched with satellite-derived features, including Sentinel-1 SAR, Sentinel-2 spectral bands, ALOS L-band SAR metrics, a digital elevation model (DEM) from JAXA, TerraClimate climate variables, forest type classes from Copernicus Global Land Service (CGLS), Dynamic World-derived forest probability, and VIIRS nighttime light intensity. These variables were selected to represent complementary aspects of canopy structure, vegetation condition, terrain-driven environmental gradients, climatic constraints, and anthropogenic disturbance. A Light Gradient Boosting Machine (LightGBM) regression model was trained to predict GEDI-derived FHD from the multi-source feature set. Model evaluation was conducted using out-of-fold (OOF) predictions to reduce optimistic bias and to provide a robust estimate of generalization performance. The resulting predictions were used to produce a spatially continuous wall-to-wall map of FHD for the entire region.
The model achieved an RMSE of 0.360 and an R² of 0.306 in predicting GEDI-derived FHD. This indicates that multi-source satellite observations capture part of the spatial variability in forest vertical structural diversity across heterogeneous landscapes.
The predicted wall-to-wall FHD map reveals spatially coherent patterns across Hokkaido. Higher values appear in mountainous forested regions, while lower values occur in flatter or more human-influenced areas. Extremely high-elevation areas show lower FHD values, likely reflecting harsher climatic conditions and simplified forest structures near the treeline.
These spatial patterns align with known ecological gradients in forest composition and management intensity. This suggests the model captures meaningful large-scale structural variability rather than random noise. To interpret the model beyond simple prediction, we used SHAP analysis to understand how individual features contribute to predicted FHD values.
SHAP-based global feature importance showed that tree cover probability (from land cover products) was the most influential predictor. This highlights the fundamental role of forest presence and canopy continuity in determining vertical structural diversity. Optical spectral bands from Sentinel-2, particularly visible and red-edge bands, contributed strongly, reflecting how spectral responses vary with canopy density and vegetation condition. SAR backscatter from ALOS L-band and Sentinel-1 also showed high importance, indicating that longer-wavelength radar signals capture canopy structure and woody biomass information relevant to vertical heterogeneity. Topographic variables (DEM) contributed significantly, suggesting that elevation-related environmental gradients influence forest structure through climate and disturbance patterns. Nighttime light intensity and its distance-based metrics showed measurable but secondary contributions, implying that human pressure is associated with reduced vertical structural diversity in more developed areas. SHAP summary plots further revealed nonlinear and asymmetric relationships between key predictors and FHD. Higher tree cover probability and stronger L-band HV backscatter were associated with positive contributions to predicted FHD, while increasing nighttime light intensity tended to reduce it. Elevation showed both positive and negative contributions depending on context, reflecting complex interactions between topography, forest type, and management practices. These results demonstrate that the machine learning model captures ecologically meaningful relationships rather than purely statistical correlations.
This study's key contribution is demonstrating that foliage height diversity can be estimated at regional scale by combining spaceborne LiDAR data with multi-modal satellite observations. The resulting model can be meaningfully interpreted in ecological terms using SHAP analysis. This interpretability distinguishes our framework from purely predictive approaches and reveals how canopy cover, radar-derived structure, topography, and human influence shape forest vertical structure.
While predictive accuracy remains moderate, reflecting the inherent challenge of inferring three-dimensional forest structure from two-dimensional satellite observations, the results show that broad regional patterns of forest vertical structural diversity can be captured consistently.