GEE：landsat5+modis；投影坐标转地理坐标；将两个栅格大小，坐标系统，像元大小配准；

目录

1、下载Modis白天数据 

下载夜间数据

2、单个 landsat5 热红外数据集下载

批量

3、下载的landsat5和土地利用坐标系统不一致

投影坐标转地理坐标

根据栅格1的坐标系统、像元大小、长和宽，依次对栅格2进行坐标系统转化，重采样，裁剪或填充，保证栅格1和栅格2坐标系统、像元大小、长和宽一样

4、单张，统计每个矢量所覆盖的温度

批量，统计每个矢量所覆盖的温度

5、单张输出每个地物类型的温度平均值

批量输出每个地物类型的温度平均值

---

增加导出影像部分

1、下载Modis白天数据 

// 使用湖北省的边界
var hubei = ee.FeatureCollection(table);

// 将地图中心设置为湖北省
Map.centerObject(hubei, 7);

// 获取一年的温度数据
var modis = ee.ImageCollection('MODIS/006/MOD11A1');
var modisLST = modis.filterBounds(hubei)
                    .filterDate('2021-01-01', '2021-12-31')
                    .select('LST_Day_1km');

// 将温度转换为摄氏度
modisLST = modisLST.map(function(img) {
  var date = img.get('system:time_start');
  return img.multiply(0.02).subtract(273.15).set('system:time_start', date);
});

// 在地图上显示平均温度图层
var meanLST = modisLST.mean().clip(hubei);
Map.addLayer(meanLST, {min: 10, max: 30, palette: ['green', 'yellow', 'red']}, 'LST');

// 导出影像
Export.image.toDrive({
  image: meanLST,
  description: 'MODIS_LST_Hubei_2021',
  scale: 1000,  // 设置导出影像的分辨率，单位为米
  region: hubei.geometry(),  // 设置导出影像的范围
  fileFormat: 'GeoTIFF',  // 设置导出影像的格式
  folder: 'GEE_exports',  // 设置导出到 Google Drive 的文件夹
  maxPixels: 1e9  // 设置导出影像的最大像素数
});

下载夜间数据

// 使用湖北省的边界
var hubei = ee.FeatureCollection(table);

// 将地图中心设置为湖北省
Map.centerObject(hubei, 7);

// 获取一年的温度数据
// 获取一年的温度数据
var modis = modis2;
var modisLST = modis.filterBounds(hubei)
                    .filterDate('2001-01-01', '2001-12-31')
                    .select('LST_Night_1km');  // 使用晚上的温度波段

// 将温度转换为摄氏度
modisLST = modisLST.map(function(img) {
  var date = img.get('system:time_start');
  return img.multiply(0.02).subtract(273.15).set('system:time_start', date);
});

// 在地图上显示平均温度图层
var meanLST = modisLST.mean().clip(hubei);
Map.addLayer(meanLST, {min: 10, max: 30, palette: ['green', 'yellow', 'red']}, 'LST');

// 导出影像
Export.image.toDrive({
  image: meanLST,
  description: 'MODIS_LST_Night_Hubei_2001',
  scale: 1000,  // 设置导出影像的分辨率，单位为米
  region: hubei.geometry(),  // 设置导出影像的范围
  fileFormat: 'GeoTIFF',  // 设置导出影像的格式
  folder: 'GEE_exports',  // 设置导出到 Google Drive 的文件夹
  maxPixels: 1e9  // 设置导出影像的最大像素数
});

2、单个 landsat5 热红外数据集下载

var roi = table //感兴趣的区域信息
var style_set = {color:"red",fillColor:"00000000"}; //设置地图中要素的颜色和填充颜色
Map.addLayer(roi.style(style_set),{},"shape") //使用之前定义的样式集将roi添加到地图中。该地图层默认使用几何形状（例如多边形）来表示区域
Map.centerObject(roi,8) //将地图中心设置为roi对象，并设置缩放级别为10
 
//本示例演示了使用Landsat 8 Collection 2，Level 2的QA_PIXEL波段（CFMask）来屏蔽不需要的像素。
 
//定义函数maskL8sr，接受一个图像作为输入，并对图像进行处理
function maskL5sr(image) {
  // Bit 0 - Fill
  // Bit 1 - Dilated Cloud
  // Bit 2 - Cirrus
  // Bit 3 - Cloud
  // Bit 4 - Cloud Shadow
  var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0);
  //从输入图像中选择QA_PIXEL波段，使用位运算和掩码来识别填充、云、云影等像素
  var saturationMask = image.select('QA_RADSAT').eq(0); //从输入图像中选择QA_RADSAT波段，并识别未饱和的像素。
  
  // 将缩放因子应用于适当的频带
  // var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  // 从输入图像中选择光学波段，并应用归一化处理。
  // var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0)
  // 从输入图像中选择热红外波段，并应用归一化处理。
  // addBands(opticalBands, null, true)
  // 用缩放的带替换原始带并应用掩码。
  return image.addBands(thermalBands, null, true)
      .updateMask(qaMask)
      .updateMask(saturationMask);
}

 
 
// 将函数映射到一年的数据上
var collection = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
                     .filterDate('1990-1-1','1990-2-1')
                     .map(maskL5sr)
                     .median() //中值合成
                     .clip(roi); //裁剪
 
// Display the results.
// Map.setCenter(-4.52, 40.29, 7);  // Iberian Peninsula
 
// print(dataset)
var img = collection.select("ST_B6") //从处理后的图像集合中选择热红外波段（'ST_B10'）
var lst = img.expression(
    'B1-273.15',
    {
        B1:img.select('ST_B6'), 
    }); //对选择的热红外波段进行计算表达式操作
    
print("LST直方图", ui.Chart.image.histogram(lst, roi, 100, 258)) //打印直方图，显示热红外波段处理后的数据分布情况
print(lst) //打印热红外波段处理后的数据
 
Map.addLayer(lst, {'min':0,'max':40,'palette':["eff3ff","c6dbef","9ecae1","6baed6","4292c6","2171b5","084594",
"fff5f0","fee0d2","fcbba1","fc9272","fb6a4a","ef3b2c","cb181d","99000d"]}, 'lst')
// 将处理后的热红外波段数据添加到地图上，并设定显示范围和颜色映射
 
function exportImage(image, roi, fileName) {  
    Export.image.toDrive({  
       image: image,  
       description: "Landsat5"+fileName,  
       fileNamePrefix: fileName,  //文件命名
       folder: "Landsat 5",  //保存的文件夹
       scale: 30,  //分辨率
       region: roi,  //研究区
       maxPixels: 1e13,  //最大像元素
       crs: "EPSG:4326"  //设置投影
   });  
 }
exportImage(lst,roi,"lst5")
 

批量

// 感兴趣的区域信息
var roi = table;
 
// 设置地图中要素的颜色和填充颜色
var style_set = { color: "red", fillColor: "00000000" };
 
// // 将roi添加到地图中
// Map.addLayer(roi.style(style_set), {}, "shape");
 
// // 将地图中心设置为roi对象，并设置缩放级别为10
// Map.centerObject(roi, 8);
 
// 定义函数maskL8sr，接受一个图像作为输入，并对图像进行处理
function maskL8sr(image) {
  // Bit 0 - Fill
  // Bit 1 - Dilated Cloud
  // Bit 2 - Cirrus
  // Bit 3 - Cloud
  // Bit 4 - Cloud Shadow
  var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0);
 
  // 从输入图像中选择QA_PIXEL波段，使用位运算和掩码来识别填充、云、云影等像素
  var saturationMask = image.select('QA_RADSAT').eq(0);
 
  // 用缩放因子应用于适当的频带
  // var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
 
  // 用缩放的带替换原始带并应用掩码
  return image.addBands(thermalBands, null, true)
      .updateMask(qaMask)
      .updateMask(saturationMask);
}
 
// 创建一个空数组来存储每个时间段的影像
var finalImages = [];
 
// 循环遍历所需的时间段
for (var i = 0; i < 20; i++) {
  var startDate = ee.Date('1990-3-01').advance(i * 3, 'month');
  var endDate = ee.Date('1995-3-01').advance((i + 1) * 3, 'month');
 
  // 为当前时间段筛选图像集合
  var collection = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterDate(startDate, endDate)
    .map(maskL8sr)
    .median()
    .clip(roi);
 
  // 选择所需的波段（例如，'ST_B6'）
  var img = collection.select('ST_B6');
 
  // 计算LST
  var lst = img.expression(
    'B1 - 273.15',
    {
      B1: img.select('ST_B6')
    }
  );
 
  // 打印LST直方图
  // print('LST直方图', ui.Chart.image.histogram(lst, roi, 100, 258));
 
  // 将当前时间段的影像添加到 finalImages 数组中
  finalImages.push(lst);
 
  // 为当前时间段导出影像
  exportImage(lst, roi, 'lst5_' + i);
}
 
// 导出影像的函数
function exportImage(image, roi, fileName) {
  Export.image.toDrive({
    image: image,
    description: fileName,
    fileNamePrefix: fileName,
    folder: "Landsat 5",
    scale: 30,
    region: roi,
    maxPixels: 1e13,
    crs: "EPSG:4326"
  });
}

3、下载的landsat5和土地利用坐标系统不一致

武汉大学遥感院黄昕&李家艺团队 (whu.edu.cn)

from osgeo import gdal

def get_raster_extent(raster_path):
    # 打开栅格文件
    raster = gdal.Open(raster_path)
    
    # 获取地理变换信息
    # 这是一个六元组，包含以下元素：
    # (左上角x坐标, 东西方向像素分辨率, 旋转(通常为0), 左上角y坐标, 旋转(通常为0), 南北方向像素分辨率(通常为负值))
    geo_transform = raster.GetGeoTransform()
    
    # 获取栅格的列数和行数
    cols = raster.RasterXSize
    rows = raster.RasterYSize
    
    # 计算栅格的四个边界坐标
    min_x = geo_transform[0]
    max_x = geo_transform[0] + geo_transform[1] * cols
    max_y = geo_transform[3]
    min_y = geo_transform[3] + geo_transform[5] * rows
    
    return (min_x, min_y, max_x, max_y)

# 栅格文件的路径
raster_path_1 = 'wendu/90-1.tif'
raster_path_2 = 'data/landuse/CLCD_v01_1990_albert.tif'

# 获取栅格范围
extent_1 = get_raster_extent(raster_path_1)
extent_2 = get_raster_extent(raster_path_2)

print('栅格1的范围:', extent_1)
print('栅格2的范围:', extent_2)

投影坐标转地理坐标

from osgeo import gdal, osr
 
def reproject_to_geographic(raster_file, output_file):
    """如果为投影坐标系统，则将栅格影像转换为地理坐标系统"""
    # 打开原始影像
    src_ds = gdal.Open(raster_file)
    src_proj = src_ds.GetProjection()
 
    # 创建空间参考对象
    src_srs = osr.SpatialReference()
    src_srs.ImportFromWkt(src_proj)
 
    # 检查是否为投影坐标系统
    if src_srs.IsProjected():
        # 定义目标坐标系统为WGS84
        dst_srs = osr.SpatialReference()
        dst_srs.ImportFromEPSG(4326)
 
        # 创建转换对象
        transform = osr.CoordinateTransformation(src_srs, dst_srs)
 
        # 使用gdal.Warp重新投影影像
        gdal.Warp(output_file, src_ds, dstSRS=dst_srs.ExportToWkt())
 
        print(f"影像已从投影坐标系统转换为地理坐标系统，输出文件：{output_file}")
    else:
        print("影像不是投影坐标系统，无需转换。")
 
# 输入和输出文件路径
input_raster = 'data/landuse/CLCD_v01_1990_albert.tif'
output_raster = 'data/landuse/CLCD_v01_1990_albert_1.tif'
 
# 调用函数
reproject_to_geographic(input_raster, output_raster)
 

根据栅格1的坐标系统、像元大小、长和宽，依次对栅格2进行坐标系统转化，重采样，裁剪或填充，保证栅格1和栅格2坐标系统、像元大小、长和宽一样

from osgeo import gdal

def align_raster(reference_raster_path, target_raster_path, output_raster_path):
    """根据参考栅格的坐标系统、像元大小、长和宽对目标栅格进行处理"""

    # 打开参考栅格和目标栅格
    ref_ds = gdal.Open(reference_raster_path)
    target_ds = gdal.Open(target_raster_path)

    # 获取参考栅格的地理变换参数、投影信息、大小
    ref_geotransform = ref_ds.GetGeoTransform()
    ref_projection = ref_ds.GetProjection()
    ref_x_size = ref_ds.RasterXSize
    ref_y_size = ref_ds.RasterYSize

    # 使用gdal.Warp对目标栅格进行处理
    gdal.Warp(output_raster_path, target_ds,
              format='GTiff',  # 输出格式
              dstSRS=ref_projection,  # 设置目标坐标系统
              xRes=ref_geotransform[1],  # 设置目标像元宽度
              yRes=abs(ref_geotransform[5]),  # 设置目标像元高度
              width=ref_x_size,  # 设置输出栅格的列数
              height=ref_y_size,  # 设置输出栅格的行数
              resampleAlg=gdal.GRA_NearestNeighbour,  # 设置重采样算法
              outputBounds=(ref_geotransform[0],  # 设置输出栅格的左上角X坐标
                            ref_geotransform[3] + ref_geotransform[5] * ref_y_size,  # 设置输出栅格的左上角Y坐标
                            ref_geotransform[0] + ref_geotransform[1] * ref_x_size,  # 设置输出栅格的右下角X坐标
                            ref_geotransform[3]))  # 设置输出栅格的右下角Y坐标

# 栅格文件路径
reference_raster_path = 'wendumodis/MODIS_LST_Hubei_2021.tif'
target_raster_path = 'data/landuse/CLCD_v01_2001_albert_hubei_Huanggang.tif'
output_raster_path = 'data/landuse/CLCD_v01_2001_albert_hubei_Huanggang_1km.tif'

# 调用函数
align_raster(reference_raster_path, target_raster_path, output_raster_path)

批量

import os
from osgeo import gdal

def align_raster(reference_raster_path, target_raster_path, output_raster_path):
    """根据参考栅格的坐标系统、像元大小、长和宽对目标栅格进行处理"""

    # 打开参考栅格和目标栅格
    ref_ds = gdal.Open(reference_raster_path)
    target_ds = gdal.Open(target_raster_path)

    # 获取参考栅格的地理变换参数、投影信息、大小
    ref_geotransform = ref_ds.GetGeoTransform()
    ref_projection = ref_ds.GetProjection()
    ref_x_size = ref_ds.RasterXSize
    ref_y_size = ref_ds.RasterYSize

    # 使用gdal.Warp对目标栅格进行处理
    gdal.Warp(output_raster_path, target_ds,
              format='GTiff',  # 输出格式
              dstSRS=ref_projection,  # 设置目标坐标系统
              xRes=ref_geotransform[1],  # 设置目标像元宽度
              yRes=abs(ref_geotransform[5]),  # 设置目标像元高度
              width=ref_x_size,  # 设置输出栅格的列数
              height=ref_y_size,  # 设置输出栅格的行数
              resampleAlg=gdal.GRA_NearestNeighbour,  # 设置重采样算法
              outputBounds=(ref_geotransform[0],  # 设置输出栅格的左上角X坐标
                            ref_geotransform[3] + ref_geotransform[5] * ref_y_size,  # 设置输出栅格的左上角Y坐标
                            ref_geotransform[0] + ref_geotransform[1] * ref_x_size,  # 设置输出栅格的右下角X坐标
                            ref_geotransform[3]))  # 设置输出栅格的右下角Y坐标

def process_folder(reference_raster_path, target_folder_path, output_folder_path):
    """遍历目标文件夹中的所有栅格文件，批量处理"""

    # 确保输出文件夹存在
    if not os.path.exists(output_folder_path):
        os.makedirs(output_folder_path)

    # 遍历文件夹中的每一个文件
    for filename in os.listdir(target_folder_path):
        if filename.endswith('.tif'):  # 仅处理 .tif 文件
            target_raster_path = os.path.join(target_folder_path, filename)
            output_raster_path = os.path.join(output_folder_path, f'aligned_{filename}')

            # 调用 align_raster 函数处理每个文件
            align_raster(reference_raster_path, target_raster_path, output_raster_path)
            print(f'Processed {filename} and saved as aligned_{filename}')

# 栅格文件路径
reference_raster_path = 'D:/project/地热研究/咸宁/d到断层距离/d到断层距离.tif'
target_folder_path = 'D:/xianningdata1'  # 目标文件夹路径
output_folder_path = 'D:/xianningdata_peizhun'  # 输出文件夹路径

# 批量处理文件夹中的所有栅格文件
process_folder(reference_raster_path, target_folder_path, output_folder_path)

 

 4、单张，统计每个矢量所覆盖的温度

from osgeo import gdal, ogr
import numpy as np

# 读取矢量文件和温度栅格文件
vector_file = 'dataset/real_value_geo.shp'
raster_file = 'wendu/drive-download-20240308T030039Z-001/MODIS_LST_Day_Hubei_2001_1.tif'

# 打开矢量文件
vector_ds = ogr.Open(vector_file)
layer = vector_ds.GetLayer()

# 打开栅格文件
raster_ds = gdal.Open(raster_file)
band = raster_ds.GetRasterBand(1)
geo_transform = raster_ds.GetGeoTransform()
inv_geo_transform = gdal.InvGeoTransform(geo_transform)
nodata_value = band.GetNoDataValue()

# 遍历矢量数据中的每个多边形
for feature in layer:
    # 获取多边形的几何对象
    geom = feature.GetGeometryRef()
    polygon_id = feature.GetFID()

    # 计算多边形对应的栅格坐标范围
    min_x, max_x, min_y, max_y = geom.GetEnvelope()
    x_off, y_off = gdal.ApplyGeoTransform(inv_geo_transform, min_x, max_y)
    x_size = int((max_x - min_x) / geo_transform[1]) + 1
    y_size = int((max_y - min_y) / -geo_transform[5]) + 1

    # 裁剪温度栅格
    temp_array = band.ReadAsArray(x_off, y_off, x_size, y_size)

    # 计算裁剪区域内的平均温度，忽略空值
    valid_temps = temp_array[temp_array != nodata_value]
    if valid_temps.size > 0:
        avg_temp = np.mean(valid_temps)
        print(f'Polygon ID: {polygon_id},Average Temperature: {avg_temp}')
    else:
        print('No valid temperature data')

# 关闭数据源
vector_ds = None
raster_ds = None

批量，统计每个矢量所覆盖的温度

from osgeo import gdal, ogr
import numpy as np
import os

# 读取矢量文件和温度栅格文件
vector_file = 'dataset/real_value_geo.shp'
folder_path = 'wendu/drive-download-20240308T030039Z-001'

# 打开矢量文件
vector_ds = ogr.Open(vector_file)
layer = vector_ds.GetLayer()

# 获取文件夹中所有的 TIFF 文件，并按文件名排序
raster_files = [f for f in os.listdir(folder_path) if f.endswith('.tif')]
raster_files.sort()

for filename in raster_files:
    if filename.endswith('.tif'):  # 假设所有栅格文件都是 TIFF 格式
        file_path = os.path.join(folder_path, filename)
        print(f'Processing {filename}...')

    # 打开栅格文件
    raster_ds = gdal.Open(file_path)
    band = raster_ds.GetRasterBand(1)
    geo_transform = raster_ds.GetGeoTransform()
    inv_geo_transform = gdal.InvGeoTransform(geo_transform)
    nodata_value = band.GetNoDataValue()

    # 遍历矢量数据中的每个多边形
    for feature in layer:
        # 获取多边形的几何对象
        geom = feature.GetGeometryRef()
        # polygon_id = feature.GetFID()

        # 计算多边形对应的栅格坐标范围
        min_x, max_x, min_y, max_y = geom.GetEnvelope()
        x_off, y_off = gdal.ApplyGeoTransform(inv_geo_transform, min_x, max_y)
        x_size = int((max_x - min_x) / geo_transform[1]) + 1
        y_size = int((max_y - min_y) / -geo_transform[5]) + 1

        # 裁剪温度栅格
        temp_array = band.ReadAsArray(x_off, y_off, x_size, y_size)

        # 计算裁剪区域内的平均温度，忽略空值
        valid_temps = temp_array[temp_array != nodata_value]
        if valid_temps.size > 0:
            avg_temp = np.mean(valid_temps)
            # print(f'Polygon ID: {polygon_id},Average Temperature: {avg_temp}')
            print(avg_temp)
        else:
            print('No valid temperature data')

# 关闭数据源
vector_ds = None
raster_ds = None

5、单张输出每个地物类型的温度平均值

import numpy as np
from osgeo import gdal

# 读取土地覆盖类型栅格和温度栅格
land_cover_file = 'data/landuse/CLCD_v01_2001_albert_hubei_Huanggang_1km.tif'
temperature_file = 'wendu/drive-download-20240308T030039Z-001/MODIS_LST_Day_Hubei_2001_5.tif'

# 打开土地覆盖类型栅格
land_cover_ds = gdal.Open(land_cover_file)
land_cover_band = land_cover_ds.GetRasterBand(1)
land_cover_array = land_cover_band.ReadAsArray()

# 打开温度栅格
temperature_ds = gdal.Open(temperature_file)
temperature_band = temperature_ds.GetRasterBand(1)
temperature_array = temperature_band.ReadAsArray()

# 初始化用于存储结果的字典
average_temperatures = {}

# 对于每个土地覆盖类型 (1 到 9)，计算对应的温度栅格的平均值
for land_cover_type in range(1, 10):
    mask = land_cover_array == land_cover_type
    masked_temperatures = temperature_array[mask]
    # 使用 np.nanmean 来忽略 NaN 值
    if np.count_nonzero(~np.isnan(masked_temperatures)) > 0:
        average_temperature = np.nanmean(masked_temperatures)
        average_temperatures[land_cover_type] = average_temperature
    else:
        average_temperatures[land_cover_type] = 'No data'

# 输出结果
for land_cover_type, avg_temp in average_temperatures.items():
    print(f'Type {land_cover_type}: {avg_temp}')

# 关闭数据源
land_cover_ds = None
temperature_ds = None

批量输出每个地物类型的温度平均值

import os
import numpy as np
from osgeo import gdal

# 读取土地覆盖类型栅格
land_cover_file = 'data/landuse/CLCD_v01_2001_albert_hubei_Huanggang_1km.tif'

# 定义温度文件夹的路径
temperature_folder = 'wendu/drive-download-20240308T030039Z-001/'

# 打开土地覆盖类型栅格
land_cover_ds = gdal.Open(land_cover_file)
land_cover_band = land_cover_ds.GetRasterBand(1)
land_cover_array = land_cover_band.ReadAsArray()

temperature_files = [f for f in os.listdir(temperature_folder) if f.endswith('.tif')]
temperature_files.sort()

# 遍历排序后的温度文件列表
for filename in temperature_files:
    temperature_file = os.path.join(temperature_folder, filename)    
    # 打开温度栅格
    temperature_ds = gdal.Open(temperature_file)
    temperature_band = temperature_ds.GetRasterBand(1)
    temperature_array = temperature_band.ReadAsArray()
    
    # 初始化用于存储结果的字典
    average_temperatures = {}
    
    # 对于每个土地覆盖类型 (1 到 9)，计算对应的温度栅格的平均值
    for land_cover_type in range(1, 10):
        mask = land_cover_array == land_cover_type
        masked_temperatures = temperature_array[mask]
        if np.count_nonzero(~np.isnan(masked_temperatures)) > 0:
            average_temperature = np.nanmean(masked_temperatures)
            average_temperatures[land_cover_type] = average_temperature
        else:
            average_temperatures[land_cover_type] = 'No data'
    
    # 输出结果
    print(f'Average Temperatures for {filename}:')
    for land_cover_type, avg_temp in average_temperatures.items():
        print(avg_temp)
    
    # 关闭温度数据源
    temperature_ds = None

# 关闭土地覆盖数据源
land_cover_ds = None

6、将两个栅格数据相减，并将结果保存为一个新的栅格文件

from osgeo import gdal

# 读取栅格1和栅格2
raster1_file = 'path/to/your/raster1.tif'
raster2_file = 'path/to/your/raster2.tif'

# 打开栅格文件
raster1_ds = gdal.Open(raster1_file)
raster2_ds = gdal.Open(raster2_file)

# 读取栅格数据
raster1_band = raster1_ds.GetRasterBand(1)
raster2_band = raster2_ds.GetRasterBand(1)

raster1_array = raster1_band.ReadAsArray()
raster2_array = raster2_band.ReadAsArray()

# 计算栅格差异
difference_array = raster1_array - raster2_array

# 创建新栅格文件来保存结果
driver = gdal.GetDriverByName('GTiff')
out_file = 'path/to/your/difference.tif'
out_ds = driver.Create(out_file, raster1_ds.RasterXSize, raster1_ds.RasterYSize, 1, raster1_band.DataType)
out_band = out_ds.GetRasterBand(1)

# 写入差异数据
out_band.WriteArray(difference_array)

# 设置地理变换和投影信息
out_ds.SetGeoTransform(raster1_ds.GetGeoTransform())
out_ds.SetProjection(raster1_ds.GetProjection())

# 关闭数据源
raster1_ds = None
raster2_ds = None
out_ds = None