栅格数据用 GPU 计算

前几天做了一下土地覆被二级分类转一级分类，操作就是把二级类的像元值的十位数提取出来，作为一级类的数值，这个可以在QGIS中完成，但是转换后的数据数据量极大，在ChatGPT的帮助下，我做了个Python GPU的计算，计算结果非常理想，计算速度快，数据量很小。

QGIS重分类

QGIS的重分类是我经常使用的工具，可以完成土地覆被重分类工作，但是经实验发现，QGIS转换速度很慢，转换后的文件数据量巨大，全国30米土地覆被一幅高大20GB，转换效果差。

Python设置代码如下：

processing.run("native:reclassifybytable", {'INPUT_RASTER':'G:/Geodata/CASLUCC/1980s.tif','RASTER_BAND':1,'TABLE':['10','12','1','20','24','2','30','33','3','40','46','4','50','53','5','60','67','6','90','99','9'],'NO_DATA':-9999,'RANGE_BOUNDARIES':0,'NODATA_FOR_MISSING':False,'DATA_TYPE':0,'OUTPUT':'I:/LUCC_CAS/L1_1980s.tif'})

Python CPU栅格计算

我用的服务器是至强E5双路CPU，256GB内存，运行下面的程序，居然显示内存不足了，ChatGPT表示并行计算建议将分块设置的小一点，下面代码默认是1024，我估计256可能会好一些吧，计算速度极慢，等了好半天，不想等了，突然想起GPU也是可以计算的嘛，干嘛不试试GPU计算，于是就有了后面第三部分的代码。

import os
import numpy as np
import rasterio
from concurrent.futures import ProcessPoolExecutor
from tqdm import tqdm  # 导入tqdm库

# Define input and output directories
input_dir = r'G:\Geodata\CASLUCC'
output_dir = r'I:\LUCC_CAS'

# Create output directory if it doesn't exist
os.makedirs(output_dir, exist_ok=True)

# Get list of all .tif files in the input directory
tif_files = [f for f in os.listdir(input_dir) if f.endswith('.tif')]

# Function to process each raster file
def process_raster(input_path, output_path):
    with rasterio.open(input_path) as src:
        # Read metadata
        meta = src.meta.copy()

        # Initialize the processed data array
        processed_data = np.zeros((src.height, src.width), dtype=np.uint8)
        
        # Define block size (you can adjust this based on available memory)
        block_size = 1024  # 建议改256，减小块减小内存使用
        
        # Calculate total number of blocks
        num_blocks_x = (src.width + block_size - 1) // block_size  # Horizontal blocks
        num_blocks_y = (src.height + block_size - 1) // block_size  # Vertical blocks

        # Process the raster in blocks to reduce memory usage
        for i in tqdm(range(0, src.height, block_size), desc=f'Processing {os.path.basename(input_path)}', total=num_blocks_y):
            for j in range(0, src.width, block_size):
                window = rasterio.windows.Window(j, i, block_size, block_size)
                data = src.read(1, window=window)
                
                # Create a mask for valid data (values between 11 and 99)
                mask = (data >= 11) & (data <= 99)

                # Process data (compute tens digit)
                processed_data[i:i+block_size, j:j+block_size] = np.where(mask, data // 10, 0)

        # Update metadata for output
        meta.update(dtype=rasterio.uint8, compress='LZW', tiled=True, bigtiff='YES')

        # Write the processed data to the output file with high compression
        with rasterio.open(output_path, 'w', **meta) as dst:
            dst.write(processed_data, 1)

# Use ProcessPoolExecutor to process each .tif file in parallel
def parallel_processing():
    with ProcessPoolExecutor() as executor:
        # Map the process_raster function to each file in the list, and track progress with tqdm
        futures = [
            executor.submit(process_raster, os.path.join(input_dir, tif_file), os.path.join(output_dir, tif_file))
            for tif_file in tif_files
        ]
        
        # Use tqdm to track the overall progress of the futures
        for future in tqdm(futures, desc="Overall progress", total=len(futures)):
            future.result()  # This will raise exceptions if there were any errors during processing

# Run the parallel processing
if