用gdal和python计算光栅参数

2024-03-28 19:56:52 发布

您现在位置:Python中文网/ 问答频道 /正文

我试图创建一个脚本,它读取一个光栅数组并计算一个参数,但是我得到了一个错误:操作数不能与形状一起广播(27184310)(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()

谢谢你的帮助


Tags: columnsimageimportnewoutputosdirpaths
1条回答
网友
1楼 · 发布于 2024-03-28 19:56:52

据我所知,你要做以下事情:

您有三个不同的.tif文件,每个文件包含一个波段。将这些图像的内容读入(numpy)数组,并将它们存储在一个数组列表中,称为bands。所以band[0]是etp的数组,band[1]是雨的数组,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

上面的例子处理的是1D数组,但是2D数组也是如此。现在,如果python告诉您它不能将一个数组“广播”到另一个数组中,那么这意味着逐项计算失败,因为数组的形状不同。在上面的示例中,如果x的长度为5,而y的长度为4,则numpy/python无法知道如何处理丢失的项并抛出一个错误。在

你可以通过查看你的带子项目的形状来检查(确保你import numpy as np

^{pr2}$

可能他们会有所不同,这会导致你的错误。因此,在你执行像素代数之前,首先你必须确保你所有的图像都有相同的大小。在

相关问题 更多 >