我一直在尝试检查DEM栅格上的过滤器以进行模式识别，并且它总是导致缺少最后一行和最后一列（例如.20）。
我尝试了PIL库和图像加载。然后用numpy。输出是相同的。

我想，我的循环出了点问题，当检查数组中的值时（只是用ArcCatalog中的Identification拾取像素），我意识到像素值未加载到数组中。

因此，只需打开，放入数组并保存数组中的图像即可：

a=numpy.array(Image.open(inraster)) #raster is .tif Float32, size 561x253
newIm=Image.new(Im.mode, Im.size)
Image.fromarray(a).save(outraster)


结果是切掉了最后的行和列。抱歉，不能发布图片

任何人都可以帮助您理解为什么？并建议一些解决方案吗？
编辑：我想这是关于numpy数组的限制，所以数组会像这样自动重塑或变形……所以例如：

Traceback (most recent call last):
  File "<pyshell#36>", line 1, in <module>
    ima=numpy.array(inDs.GetRasterBand(1).ReadAsArray())
  File "C:\Python25\lib\site-packages\osgeo\gdal.py", line 835, in ReadAsArray
    buf_xsize, buf_ysize, buf_obj )
  File "C:\Python25\lib\site-packages\osgeo\gdal_array.py", line 140, in BandReadAsArray
    ar = numpy.reshape(ar, [buf_ysize,buf_xsize])
  File "C:\Python25\lib\site-packages\numpy\core\fromnumeric.py", line 108, in reshape
    return reshape(newshape, order=order)
ValueError: total size of new array must be unchanged


重点是我不想要逐块读取数据，因为我需要过滤，请使用不同的过滤器，不同的大小多次。

#1 楼

如果您具有python-gdal绑定：

import numpy as np
from osgeo import gdal
ds = gdal.Open("mypic.tif")
myarray = np.array(ds.GetRasterBand(1).ReadAsArray())


您已完成：

myarray.shape
(2610,4583)
myarray.size
11961630
myarray
array([[        nan,         nan,         nan, ...,  0.38068664,
     0.37952521,  0.14506227],
   [        nan,         nan,         nan, ...,  0.39791253,
            nan,         nan],
   [        nan,         nan,         nan, ...,         nan,
            nan,         nan],
   ..., 
   [ 0.33243281,  0.33221543,  0.33273876, ...,         nan,
            nan,         nan],
   [ 0.33308044,  0.3337177 ,  0.33416209, ...,         nan,
            nan,         nan],
   [ 0.09213851,  0.09242494,  0.09267616, ...,         nan,
            nan,         nan]], dtype=float32)


评论

是的，使用gdal我想我没有问题，但是我正尝试使用较少的库...而且numpy在“搜索时”看来非常受欢迎。确实有任何想法，为什么numpy / PIL停止加载？？？

– najuste
2012-09-10 7:51


我不知道。 PIL应该足够健壮，因此它随python一起提供。但是，imho geotiff不仅仅是图像-例如，它们包含大量元数据-并且PIL并不是（再次恕我直言）正确的工具。

–尼克斯
2012年9月10日上午11:37
打开数据时，我有时有时会讨厌那些引号和斜杠要求。但是将numpy数组写回Raster呢？它与PIL库一起使用，但是使用outputRaster.GetRasterBand（1）.WriteArray（myarray）会生成无效的栅格。

– najuste
2012年9月10日在21:44
不要忘记使用outBand.FlushCache（）将数据刷新到磁盘。您可以在此处找到一些教程：gis.usu.edu/~chrisg/python/2009

–尼克斯
2012-09-10 21:53
检查“ lists.osgeo.org/pipermail/gdal-dev/2010-January/023309.html”-看来您已经用光了内存。

–尼克斯
2012年9月23日在16:40


#2 楼

您可以使用rasterio与NumPy数组接口。要读取数组中的栅格：

import rasterio

with rasterio.open('/path/to/raster.tif', 'r') as ds:
    arr = ds.read()  # read all raster values

print(arr.shape)  # this is a 3D numpy array, with dimensions [band, row, col]


这会将所有内容读取到3D numpy数组arr中，其尺寸为[band, row, col]。


以下是一个高级示例，可以读取，编辑像素，然后将其保存回栅格中：

with rasterio.open('/path/to/raster.tif', 'r+') as ds:
    arr = ds.read()  # read all raster values
    arr[0, 10, 20] = 3  # change a pixel value on band 1, row 11, column 21
    ds.write(arr)


栅格将被写入并在结尾处关闭“ with”语句。

评论

当我写print（arr）时为什么看不到所有值。它使用...，...，？分隔值。

–穆斯塔法·乌萨尔
17年11月22日在14:11
@MustafaUçar这是NumPy打印数组的方式，您可以对其进行修改。或切片阵列的一个窗口进行打印，还有许多其他的Numpy技巧。

– Mike T
17年11月22日在22:29
一个普遍的问题。如果要输出具有多个场景的单个数组，且场景为四个维度（场景，高度，宽度，带），则应如何修改此代码段？

–里卡多·巴罗斯·洛伦索
17年12月19日在20:15
@RicardoBarrosLourenço我猜您的第四个维度（场景？）存储在每个文件中。我先填充一个空的4D numpy数组，然后遍历每个文件（场景）并插入每个文件的3D部分。您可能需要arr.transpose（（1、2、0））从每个文件中获取（高度，宽度，带）。

– Mike T
17/12/19在22:33


@MikeT这个人口会像np.append（）吗？

–里卡多·巴罗斯·洛伦索
17年12月19日在22:37

#3 楼

当然，我正在阅读一个普通的旧png图像，但是使用scipy可以正常工作（尽管imsave使用PIL）：

评论

谢谢，但是我正在尝试使用更少的库。.现在，我可以使用gdal + numpy ...（希望没有PIL）实现它。

– najuste
2012-09-10 21:47
@najuste正在使用什么操作系统？ Mac和大多数Linux版本都带有scipy和numpy。

–乍得·库珀
2012-09-10 22:18
显然...我在Windows上使用Win的各种版本。 ：/

– najuste
2012年9月11日7:31

#4 楼

我使用gdal的解决方案如下所示。我认为它非常可重用。

import gdal
import osgeo.gdalnumeric as gdn

def img_to_array(input_file, dim_ordering="channels_last", dtype='float32'):
    file  = gdal.Open(input_file)
    bands = [file.GetRasterBand(i) for i in range(1, file.RasterCount + 1)]
    arr = np.array([gdn.BandReadAsArray(band) for band in bands]).astype(dtype)
    if dim_ordering=="channels_last":
        arr = np.transpose(arr, [1, 2, 0])  # Reorders dimensions, so that channels are last
    return arr