使用 gdal 和 python 计算栅格参数

Question

我正在尝试创建一个读取栅格数组并计算参数的脚本，但出现此错误：

operands could not be broadcast together with shapes (2718,4310) (70,86).

# -*- coding: cp1252 -*-
from osgeo import gdal
import numpy as np
import sys ,os
from osgeo.gdalconst import *
output_dir='C:\wamp\www\Donnees_CRTS\irrisat_dessimination-smiej'
etp = os.path.join(output_dir,"L8_L8_ETpot_24_30m_2016_223.tif")
rain = os.path.join(output_dir,"GLDAS_NOAH025SUBP_3H.A2016223.sumday.rain.tif")
eff = os.path.join(output_dir,"eff_v11.tif")
driver=gdal.GetDriverByName('GTIff')
driver.Register()
paths = []
paths.append(etp)
paths.append(rain)
paths.append(eff)
raster_px = []
rasters_px = []
bands = []
def open_raster(raster):       
    for i in range(len(paths)) :
       raster = gdal.Open(paths[i])
       columns = raster.RasterXSize
       lines = raster.RasterYSize
       band = raster.GetRasterBand(1).ReadAsArray(0,0,columns,lines)
       bands.append(band)
       gt=raster.GetGeoTransform()
       raster_px.append(gt)
    
       band=None  
   if raster is None :
    print ("Erreur : Impossible d'ouvrir le raster: ")
   try :
    smallest = min(raster_px)
    raster.SetGeoTransform(smallest)

    
   except:
    print ("Erreur : Impossible d'extraire la bande")

  return smallest,  columns, lines



def parameters(parameter) :
    diff = bands[0]-bands[1]
   if parameter == "IWR" :
     iwr=diff/bands[2]
return iwr



def createImage(new_image,columns,lines,smallest):

    new_image=driver.Create("iwr9.tif", columns,lines, 1, gdal.GDT_Byte)
    new_image.SetProjection("EPSG:4326")
    new_image.SetGeoTransform(smallest)
    new_image.GetRasterBand(1).WriteRaster( 0, 0, columns, lines, output_dir 
  )
  new_image=None
  return columns, lines

def main() :

smallest, columns, lines = open_raster(paths)
p = "IWR"

parameters(p)
createImage("iwr9.tif",columns, lines,smallest)
main()

Answer 1

据我了解，您执行以下操作：

您有三个不同的 .tif 文件，每个文件包含一个波段。您将这些图像的内容读入 (numpy) 数组并将它们存储在数组列表中，称为 bands。所以band[0]是etp的数组，band[1]是rain的数组，band[2]是eff的数组。这三个数组具有 x&y 维度（即图像的维度/"picture"）。

bands[1] - bands[0] 是一个数组除法，这意味着您从另一个数组的项目中减去每个数组项目的值。

例如

import numpy as np
x = np.array([10, 20, 30, 40, 50])
y = np.array([1, 2, 3, 4, 5])
z = x - y # print(z): 9 18 27 36 45

上面的例子处理的是一维数组，但同样适用于二维数组。现在如果 python 告诉你它不能 "broadcast" 一个数组到另一个数组，那么这意味着逐项计算失败，因为数组的 形状不同 .在上面的示例中，如果 x 的长度为 5 但 y 的长度为 4，则 numpy/python 无法知道如何处理丢失的项目并抛出错误。

您可以通过查看手环项目的形状来检查（确保 import numpy as np）

np.shape(bands[0])
np.shape(bands[1])
np.shape(bands[2])

它们可能会有所不同，从而导致您的错误。因此，在执行逐像素代数之前，首先您必须确保所有图像具有相同的大小。

使用 gdal 和 python 计算栅格参数

Calculating parameter on raster using gdal and python

python

gdal