python遥感影像强制重采样，gdal不同坐标系的影像如何强行校正，ArcGIS投影栅格失败如何进行栅格转投影

问题1：

我们在进行遥感影像处理的时候，有的时候会发现两个影像的行列数不一致，从而导致无法进行运算（像计算NDVI这种需要用到矩阵运算的，都必须保持相同的行列数），这时候就需要进行重采样将两景影像的行列数保持一致。

问题2：

有的时候，也会遇到这样的问题，同一区域的两景影像，坐标系不一致，比如一个坐标系是WGS84的，另一个坐标系是CGCS2000的，而在arcgis中没办法实现二者的转换。因为CGCS2000坐标系的很多参数是保密的，ArcGIS中“project raster”工具并没有内置这些参数，因此在对栅格进行校正的时候，会报错，并提示需要额外的参数，就是下图中红色框中的“Geographical Transformation”中需要我们提供转投影的参数才能转换成功。很显然我们拿不到这些参数，ArcGIS投影栅格必然失败。因此这种情况下就需要另一种思路进行重投影。

下面给大家分享一个方法，

可以将目标影像，按照参考影像的行列数，强制校正与参考影像的行列数保持一致。

同时也可以对目标影像进行重投影，将目标影像的投影与参考影像保持一致。

但是前提是，两景影像是同一区域的影像，且具有相同的外接矩形。

代码如下：

import numpy as np
import os
import netCDF4 as nc
from osgeo import gdal, osr
import matplotlib.image as mpimg # mpimg 用于读取图片

def readTifGeoTransform(fileName):
    dataset = gdal.Open(fileName)
    if dataset == None:
        print(fileName + "文件无法打开")
        return
    im_width = dataset.RasterXSize  # 栅格矩阵的列数
    im_height = dataset.RasterYSize  # 栅格矩阵的行数
    im_geotrans = dataset.GetGeoTransform()  # 获取仿射矩阵信息
    im_proj = dataset.GetProjection()  # 获取投影信息
    return im_width,im_height,im_geotrans,im_proj

def reprojection(projected_dataname,unprojected_dataname,output_dataname):
    im_width, im_height, im_geotrans, im_proj=readTifGeoTransform(projected_dataname)
    # 仅针对普通的图像的读取
    ds=mpimg.imread(unprojected_dataname)

    driver = gdal.GetDriverByName('GTiff')
    out_tif_name = output_dataname
    out_tif = driver.Create(out_tif_name, im_width, im_height, 1, gdal.GDT_Float32)  # 创建框架
    out_tif.SetGeoTransform(im_geotrans)
    out_tif.SetProjection(im_proj)
    out_tif.GetRasterBand(1).WriteArray(ds)  # 将数据写入内存，此时没有写入硬盘
    out_tif.FlushCache()  # 将数据写入硬盘
    out_tif = None  # 注意必须关闭tif文件

if __name__ == '__main__':
    reprojection(r'E:\test1\TB.tif',r'E:\test1\SZ.tif',r'E:\test1\SZ_reproject.tif')