The main workflow for karez mapping is illustrated in Fig. 2. Initially, remote sensing data of the study area were downloaded and pre-processed. Subsequently, a portion of the karez shafts was annotated for training deep learning models, and the karez shafts in the whole area were detected by using the best model. Following this, the detected shafts undergo a screening process employing various methods. After screening, the shafts were clustered to identify those belonging to individual karez systems. Ultimately, the shafts were sorted, and karez lines were generated.
Main workflow of karez mapping in this study.
Study area
In China, the karez is predominantly found in the Xinjiang Uyghur Autonomous Region, particularly in the Turpan Basin, which lies at the southern base of the Tianshan Mountain range (Fig. 3). Turpan is a renowned closed basin characterized by an extremely arid climate. The basin’s terrain exhibits substantial altitude variations, with the surrounding mountains reaching a maximum height of 5445 m at Tianshan Mountain and the basin’s lowest point, Ayding Lake, sitting at 154 m below sea level. These altitude differences significantly influence the basin’s climate and water distribution. Given the large elevation drop, the abundant meltwater from the mountains, and the arid conditions, the excavation of karez became essential for irrigation and other water needs. The historical record of the karez extends back to the Han Dynasty in China18. Notably, the karez in Turpan were recognized as one of the World Heritage Irrigation Structures on September 3, 202443. According to Kobori19, there were 1237 karez systems in Turpan in 1957, which declined to 829 in 1980, with a combined length of approximately 2360 km. During the Third National Cultural Heritage Survey conducted between 2007 and 2011, 1108 karez were archived as cultural heritage in Turpan, which consists 71.95% of total in China44.
a The location of China in the world. b The location of Turpan in China. c 3D scene of the Turpan basin (the terrain was exaggerated vertically by eight times for better visualization). d The patches for downloading images, and the selected ones for postprocessing.
Data collection and preprocessing
Bing Virtual Earth images with a spatial resolution of 0.50 m were downloaded using QGIS software version 3.38.0. This resolution allowed for clear identification of most karez shafts. The selection of the study area was based on the distribution of human activities and the presence of an oasis. Training deep learning models and performing object detection on large image datasets requires substantial computation resources. To streamline downloading, the study area was divided into 65 blocks, each measuring 30 km by 16 km, as shown in Fig. 3d. Subsequently, each of these blocks was further subdivided into nine equal smaller blocks to facilitate the detection of karez shafts using a deep learning model.
The Copernicus DEM data, featuring a spatial resolution of 90 m, was downloaded from the Copernicus website (https://dataspace.copernicus.eu/) to generate slope and aspect angle data. However, some regions were identified as flat areas with the 90 m resolution DEM data. To address this, the DEM data was resampled to 500 m resolution, allowing for the accurate derivation of aspect angle in those previously problematic areas.
Preparation of training data
The blocks numbered 32, 42, and 59 in Fig. 3d were selected for labeling karez shafts due to their representative characteristics. The shafts located in various environments were labeled utilizing the ArcGIS Pro 3.1.5 software, resulting in a total of 5825 karez shafts being marked with circular symbols, as depicted in Fig. 4. The PASCAL Visual Object Classes method within ArcGIS Pro was utilized to generate the training data. The output tile size was set to 256 pixels, with strides of 128 pixels in both the X and Y directions. No rotation angles were applied for image data augmentation. Consequently, a total of 8395 images were created, each with dimensions of 3 * 256 * 256, where three represents the red, green, and blue channels, and 256 denotes the size in both the X and Y directions. Furthermore, one image with 10 km by 5.33 km from the block of 35 was selected to test the models, in which 1607 shafts were labeled.
The green circles represent the labeled karez shafts.
Deep learning models for karez shaft detection
Object detection, a cornerstone task in computer vision, involves recognizing and locating diverse objects within visual data, thereby empowering machines to interpret and comprehend their surroundings45. Various object detection algorithms exist, and for the purpose of karez shaft detection, the deep learning library of ArcGIS Pro 3.1.5 was used. Specifically, the Faster R-CNN, SSD (Single Shot Multibox Detector), YOLOv3 (You Only Look Once version 3), and MMDetection (OpenMMLab Detection Toolbox) were implemented, each equipped with different backbone models. For Faster R-CNN, SSD, and MMDetection, ResNet34 and ResNet50 were used as the backbone models, while DarkNet53 was employed for YOLOv3. The parameters of these deep learning models were detailed in Table 1.
Faster R-CNN is an object detection model that builds upon Fast R-CNN by incorporating a Region Proposal Network (RPN) in conjunction with the CNN model46. The RPN shares convolutional features with the detection network across the entire image, thereby facilitating region proposals at minimal computational cost. As a fully convolutional network, the RPN simultaneously predicts object boundaries and foreground probability scores at each spatial position. Trained end-to-end, the RPN generates high-quality region proposals that Fast R-CNN subsequently utilizes for detection. Notably, Faster R-CNN and its RPN component were instrumental in achieving 1st-place wins in several tracks of the ILSVRC and COCO 2015 competitions.
SSD is an advanced object detection method that employs a single-stage detection network to efficiently and accurately detect objects in images. This network integrates predictions from multiscale features, enabling it to capture diverse spatial resolutions and contextual information47. Unlike two-stage detectors such as Faster R-CNN, which involve a separate region proposal stage followed by classification and bounding box refinement, SSD performs detection in a single pass, making it computationally more efficient and faster.
The YOLOv3 algorithm, introduced by Redmon and Farhadi in 2018, initially divides an image into a grid48. Each grid cell is tasked with predicting a certain number of bounding boxes (also known as anchor boxes) around objects that closely match predefined classes, assigning a high confidence score to these predictions. Each bounding box has a respective confidence score of how accurately it assumes that prediction should be. Importantly, each bounding box is designated for a single object. These bounding boxes are generated through clustering the dimensions of the ground truth boxes from the training dataset, allowing YOLOv3 to identify the most common object shapes and sizes. Unlike systems such as R-CNN and Fast R-CNN, YOLOv3 can simultaneously perform classification and bounding box regression. This dual capability ensures that YOLOv3 provides both efficient and precise predictions for the detected bounding boxes.
MMDetection is an object detection toolbox that encompasses a comprehensive collection of object detection and instance segmentation algorithms, along with pertinent components and modules. It distinguished itself by winning the detection track of the COCO Challenge 201849, showcasing its effectiveness and robustness in the field of computer vision.
After training, the model with the highest precision was selected to detect karez shafts across the entire study area. Table 2 presents the parameters employed for detecting karez shafts using the optimal model. The training of the data and the detection of karez shafts were performed on a graphic workstation equipped with a 13th Gen Intel® Core™ i9-13900KF processor clocked at 3.00 GHz, 128 GB of RAM, and an NVIDIA GeForce RTX 4080 GPU.
Precision (P), recall (R), f1-score, mean average precision at IoU = 0.50 (mAP50), and mean average precision at IoU = [0.50, 0.55, …, 0.95] (mAP50-95) were used to evaluate the models on test data using the “compute accuracy for object detection” function of ArcGIS Pro.
Post-processing of detected shafts
The initially detected karez shafts were square-shaped. For further processing, these polygons were converted into point features, ensuring that the center of each square was recorded as a single shaft. Given the presence of some incorrectly detected or missing shafts, it is essential to employ screening methods to achieve more accurate data. In this study, five steps of screening were used to eliminate incorrectly detected shafts, as outlined below:
The first step is the deletion of points too close together (<5 m) to avoid duplication. The downloaded image patches were overlapped to a certain extent to avoid missing shafts, and generally, the shafts should not be positioned too close to each other. Therefore, any shafts found within a distance of less than 5 m from another shaft were deleted to eliminate duplication. This procedure was performed using the ArcGIS Pro software.
Next step is the deletion of points on larger slopes. Based on experience, shafts located on slopes steeper than 20° were removed. The slopes were calculated using the Copernicus DEM with a spatial resolution of 90 m
In the third step, the directed fan-shaped buffer method was conducted. Typically, the karez shafts distribute linearly and extend predominantly from higher terrain to lower terrain. To accommodate this, aspect angles of the terrain were derived using the Copernicus DEM with a spatial resolution of 500 m. For each karez shaft, fan-shaped buffer zones were created, encompassing a fan angle range of aspect ±45° (Fig. 5). Recognizing that some karez lines do not strictly adhere to the aspect angle, slight adjustments were made to the angle for certain areas (this adjustment was controlled in ±20°). The largest shaft distance was selected as the buffer distance in different patches. Subsequently, only the intersected buffer zones were selected, and shafts located outside these zones were removed. This step can effectively screen out the isolated points that did not align with the direction of the karez lines and out of the largest shaft distance. This process was executed within the Python environment of ArcGIS Pro.
a Generated buffer according to aspect angle. b Selected buffers according to intersection. The colored buffer zones were the selected ones, so the non-selected ones and shaft points will be deleted.
Subsequently, small clusters were screened out based on a larger distance threshold. Utilizing the single-linkage clustering method, clusters were formed with a distance threshold set at 0.001°. Then, clusters containing three or fewer points were removed. This process was facilitated using the SciPy open-source Python library.
Finally, manual adjustments were applied. Despite the preceding steps, some erroneous points remained. Therefore, regions of interest were delineated based on expert experience, specifically excluding mountainous areas and other regions where karez are unlikely to be found. Following this, a manual screening process was conducted to eliminate points that were not linear or highly irregular, combined with satellite images. This step was executed using ArcGIS Pro software.
Line generation for karez segments
Data clustering organizes objects into groups based on their similarity. Among the hierarchical clustering methods, the method stands out by offering multiple, layered partitions. It employs divisive and agglomerative approaches to split or merge clusters, respectively, based on measures of dissimilarity or similarity that are tailored to accommodate diverse data sets and applications50. Hierarchical clustering is widely utilized in cluster analysis for data mining due to its ability to facilitate the visualization of similar objects grouped into clusters. In each iteration of hierarchical clustering, specific objects or clusters are linked together based on a criterion referred to as the “Linkage Criteria”. Single-linkage clustering is particularly advantageous when aiming to gain insights into the fine granularity of the data. It performs effectively when clusters naturally exhibit a chain-like or non-compact structure51,52. The fundamental principles of this method are outlined in the paper by Manning and Schutze51.
To identify the shafts belonging to the same karez system, the single-linkage clustering method of the hierarchical clustering approach was employed. The clustering analysis was executed using the SciPy open-source Python library. Due to variations in the distances between shafts across different regions, different maximum distance thresholds were used for clustering shafts in various blocks. However, some karez shafts were undetected either because they had been destroyed or because they were not sufficiently distinctive compared to their surroundings, resulting in shafts belonging to the same karez segment being mistakenly clustered into several separate parts. Conversely, some shafts belonging to different karez were too close to each other and were inadvertently clustered as part of the same karez. In such instances, the automatic clustering method alone was insufficient to yield accurate clusters. Therefore, manual adjustments were necessary to correct these clusters. Ultimately, the clusters from different patches were merged and assigned the same cluster IDs.
After clustering, the shafts are sorted utilizing the “sort” tool in ArcGIS Pro. The X coordinates are generated specifically for shafts where the karez lines are oriented in an east-west direction and are subsequently employed for the sorting process. Ultimately, the “Points to Line” tool was utilized to generate the lines for karez segments.
Field verification
A comprehensive field verification campaign was conducted from 28th October 2024 to 10th April 2025 as an integral component of the fourth National Cultural Heritage survey. The investigation encompassed all identified karez shafts within the study area through systematic ground-truthing procedures. Four specialized survey teams were deployed, with each team consisting of four qualified members. All members were equipped with mobile devices running the Omap (Ovi Interactive Map) application (Beijing Yuansheng Huawang Software Co., Beijing, China), enabling real-time visualization of detected shaft locations and navigation to target sites via satellite imagery. To ensure complete documentation of the extensive karez networks, each team was equipped with a DJI Mavic 3 drone system (DJI Innovations Science and Technology Co., Ltd, Shenzhen, China) for aerial photography, particularly crucial for shafts extending several kilometers in length. Every identified shaft underwent rigorous field verification, with precise location data collected and detailed distribution maps produced based on direct field observations.
