1

首先，这是脚本：

import numpy as np
import osgeo.gdal
import os

ArbitraryXCoords = np.arange(435531.30622598,440020.30622598,400)
ArbitraryYCoords = np.arange(5634955.28972479,5638945.28972479,400)

os.chdir('/home/foo/GIS_Summer2013')
dataset = osgeo.gdal.Open("Raster_DEM")
gt = dataset.GetGeoTransform()

def XAndYArrays(spacing):
    XPoints = np.arange(gt[0], gt[0] + dataset.RasterXSize * gt[1], spacing)
    YPoints = np.arange(gt[3] + dataset.RasterYSize * gt[5], gt[3], spacing)
    return (XPoints, YPoints)

def RasterPoints(XCoords,YCoords):
    a=[]
    for row in YCoords:
        for col in XCoords:
            rasterx = int((col - gt[0]) / gt[1])
            rastery = int((row - gt[3]) / gt[5])
            band = int(dataset.GetRasterBand(1).ReadAsArray(rasterx,rastery, 1, 1)[0][0])
            a[len(a):] = [band]
    foo = np.asarray(a)
    bar = foo.reshape(YCoords.size,XCoords.size)
    return bar

当我加载上面提供的脚本时，我无法将函数 XAndYArrays 的输出用作函数 RasterPoints 的输入。但我可以使用手动定义的 numpy.ndarray 作为函数 RasterPoints 中的输入。但这还不够好。我需要能够将 XAndYArrays 的输出用作 RasterPoints 中的输入。以下是我在 PyDev 交互式控制台中使用的命令：

>>> Eastings,Northings = XAndYArrays(400)
>>> Eastings
Out[1]: 
array([ 435530.30622598,  435930.30622598,  436330.30622598,
        436730.30622598,  437130.30622598,  437530.30622598,
        437930.30622598,  438330.30622598,  438730.30622598,
        439130.30622598,  439530.30622598,  439930.30622598])
>>> Northings
Out[1]: 
array([ 5634954.28972479,  5635354.28972479,  5635754.28972479,
        5636154.28972479,  5636554.28972479,  5636954.28972479,
        5637354.28972479,  5637754.28972479,  5638154.28972479,
        5638554.28972479,  5638954.28972479])
>>> RasterPoints(Eastings, Northings)
ERROR 5: MergedDEM_EPSG3159_Reduced, band 1: Access window out of range in RasterIO().  Requested (0,246) of size 1x1 on raster of 269x246.
Traceback (most recent call last):
  File "/usr/lib/python2.7/dist-packages/IPython/core/interactiveshell.py", line 2538, in run_code
    exec code_obj in self.user_global_ns, self.user_ns
  File "<ipython-input-1-326be9918188>", line 1, in <module>
    RasterPoints(Eastings, Northings)
  File "/home/foo/GIS_Summer2013/src/22July_StackOverflowQuestion.py", line 23, in RasterPoints
    band = int(dataset.GetRasterBand(1).ReadAsArray(rasterx,rastery, 1, 1)[0][0])
TypeError: 'NoneType' object has no attribute '__getitem__'
>>> RasterPoints(ArbitraryXCoords, ArbitraryYCoords)
Out[1]: 
array([[422, 422, 431, 439, 428, 399, 410, 395, 398, 413, 409, 386],
       [414, 428, 421, 430, 426, 403, 409, 410, 406, 408, 412, 406],
       [420, 428, 427, 424, 408, 406, 428, 420, 408, 410, 409, 420],
       [392, 418, 426, 430, 414, 428, 430, 418, 433, 414, 402, 399],
       [400, 411, 420, 406, 401, 405, 398, 420, 419, 400, 401, 414],
       [408, 421, 418, 428, 399, 398, 405, 412, 421, 406, 395, 397],
       [399, 404, 398, 401, 400, 399, 399, 398, 398, 419, 399, 395],
       [401, 410, 407, 407, 404, 400, 398, 397, 397, 399, 400, 398],
       [400, 410, 418, 405, 401, 400, 397, 398, 400, 398, 397, 396],
       [389, 387, 399, 408, 423, 400, 407, 398, 411, 408, 410, 420]])
>>> print "partial success"
partial success
4

1 回答 1

0

您正在尝试读取一个不存在的像素位置，因为它超出了光栅尺寸。

尝试使用以下方法计算您的 Ypoints：

YPoints = np.arange(gt[3] + (ds.RasterYSize-1) * gt[5], gt[3], abs(gt[5]))

您的RasterPoints功能确实是不好的做法。您当时正在访问所有数组值 1 ，将它们存储在一个列表中，然后创建一个数组。这根本没有意义。

此外，我认为它的良好做法是在打开另一个问题之前关闭/解决以前的问题。

于 2013-07-23T08:42:26.770 回答