Raster vollständig in ein Numpy-Array laden?

26

Ich habe versucht, meine Filter im DEM-Raster auf Mustererkennung zu überprüfen. Dabei fehlen immer die letzten Zeilen und Spalten (wie z . B. 20) . Ich habe mit PIL-Bibliothek versucht, Bilder zu laden. Dann mit Numpy. Die Ausgabe ist die gleiche.

Ich dachte, mit meinen Schleifen stimmt etwas nicht, als ich die Werte im Array überprüfte (indem ich nur die Pixel mit "Identifikation" in ArcCatalog auswählte), stellte ich fest, dass die Pixelwerte nicht in ein Array geladen wurden.

Öffnen Sie einfach das Array und speichern Sie das Bild aus dem Array:

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

Entfernt die letzten Zeilen und Spalten. Das Bild kann leider nicht # gepostet werden

Kann jemand helfen zu verstehen, warum? Und eine Lösung empfehlen?

BEARBEITEN:

Ich habe es also geschafft, mit Hilfe von Leuten kleine Raster in Numpy-Arrays zu laden, aber wenn ich ein größeres Bild habe, bekomme ich Fehler. Ich nehme an, es geht um die Grenzen des Numpy-Arrays, und so wird das Array automatisch umgeformt oder ähnliches ... Also Bsp .:

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

Der Punkt ist, dass ich Block für Block nicht lesen möchte, da ich mehrmals mit verschiedenen Filtern und verschiedenen Größen filtern muss

najuste
quelle

Antworten:

42

Wenn Sie Python-GDAL-Bindungen haben:

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

Und du bist fertig:

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)
nickves
quelle
Ja, mit gdal hatte ich wohl kein Problem, aber ich versuche, so wenig Bibliotheken wie möglich zu verwenden. Irgendeine Idee, in der Tat, warum das Laden von numpy / PIL aufhört ???
Najuste
Ich weiß es nicht. PIL sollte robust genug sein, damit es mit Python ausgeliefert wird. Aber imho Geotiff ist mehr als nur Bilder - sie enthalten zum Beispiel viele Metadaten - und PIL ist (wieder imho) nicht das richtige Werkzeug.
nickves
Ich hasse nur manchmal die Anforderungen für diff-Anführungszeichen und Schrägstriche beim Öffnen von Daten. Aber was ist mit dem Zurückschreiben von numpy-Arrays in Raster? Es arbeitet mit PIL - Bibliothek, aber unter Verwendung von outputRaster.GetRasterBand (1) .WriteArray (myarray) erzeugt ungültige Raster ..
najuste
Vergessen Sie nicht, die Daten mit outBand.FlushCache () auf die Festplatte zu kopieren. Einige Tutorials finden Sie hier: gis.usu.edu/~chrisg/python/2009
nickves
21

Mit rasterio können Sie eine Schnittstelle zu NumPy-Arrays herstellen. So lesen Sie ein Raster in ein Array:

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]

Dadurch wird alles in ein 3D-Zahlenfeld arrmit Abmessungen eingelesen [band, row, col].


Hier ist ein erweitertes Beispiel zum Lesen, Bearbeiten und Speichern eines Pixels im Raster:

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)

Das Raster wird am Ende der "with" -Anweisung geschrieben und geschlossen .

Mike T
quelle
Warum können wir nicht alle Werte sehen, wenn ich print (arr) schreibe? Es trennt Werte mit diesem ..., ...,?
Mustafa Uçar
@ MustafaUçar So druckt NumPy Arrays, die Sie ändern können . Oder schneiden Sie unter vielen anderen Numpy-Tricks ein Fenster des zu druckenden Arrays.
Mike T
Eine allgemeine Frage. Wenn ich ein einzelnes Array mit mehreren Szenen als vier Dimensionen (Szene, Höhe, Breite, Bänder) ausgeben möchte, wie soll ich dieses Snippet ändern?
Ricardo Barros Lourenço
@ RicardoBarrosLourenço Ich vermute, deine vierte Dimension (Szene?) Ist in jeder Datei gespeichert. Ich würde zuerst ein leeres 4D Numpy Array füllen, dann jede Datei (Szene) durchlaufen und den 3D-Teil von jeder einfügen. Möglicherweise müssen arr.transpose((1, 2, 0))Sie (Höhe, Breite, Bänder) aus jeder Datei abrufen.
Mike T
@MikeT diese Bevölkerung wäre wie np.append()?
Ricardo Barros Lourenço
3

Zugegeben, ich lese ein einfaches altes PNG-Bild, aber das funktioniert mit scipy ( imsaveverwendet jedoch PIL):

>>> import scipy
>>> import numpy
>>> img = scipy.misc.imread("/home/chad/logo.png")
>>> img.shape
(81, 90, 4)
>>> array = numpy.array(img)
>>> len(array)
81
>>> scipy.misc.imsave('/home/chad/logo.png', array)

Meine resultierende PNG ist auch 81 x 90 Pixel.

Chad Cooper
quelle
Danke, aber ich versuche, so wenig Bibliotheken wie möglich zu verwenden. Und jetzt kann ich es mit gdal + numpy schaffen (hoffentlich ohne PIL).
Najuste
1
@najuste Welches Betriebssystem ist aktiviert? Mac und die meisten Linux-Versionen kommen mit scipyund numpy.
Chad Cooper
Anscheinend ... habe ich unter Windows verschiedene Versionen von Win. : /
najuste
2

Meine Lösung mit gdal sieht so aus. Ich denke, dass es sehr wiederverwendbar ist.

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
MonsterMax
quelle
0

Ich verwende ein hyperspektrales Bild mit 158 ​​Bändern. Ich möchte das Raster berechnen. aber ich verstehe

import gdal # Import GDAL library bindings
from osgeo.gdalnumeric import *
from osgeo.gdalconst import *
import pylab as plt
import numpy as np
import xlrd
# The file that we shall be using
# Needs to be on current directory
filename = ('C:/Users/KIFF/Desktop/These/data/Hyperion/10th_bandmathref')
outFile = ('C:/Users/KIFF/Desktop/These/data/Hyperion/Math')
XLS=('C:/Users/KIFF/Desktop/These/data/Coef/bcoef.xlsx')
wb = xlrd.open_workbook(XLS)
sheet = wb.sheet_by_index(0)
sheet.cell_value(0, 0)


g = gdal.Open(filename, GA_ReadOnly)

# g should now be a GDAL dataset, but if the file isn't found
# g will be none. Let's test this:
if g is None:
    print ("Problem opening file %s!" % filename)
else:
    print ("File %s opened fine" % filename )

#band_array = g.ReadAsArray()
#print(band_array)
print ("[ RASTER BAND COUNT ]: ", g.RasterCount)

for band in range( g.RasterCount ):
    print (band)
    band += 1
    outFile = ('C:/Users/KIFF/Desktop/These/data/Results/Temp/Math_1_sur_value'+str(band)+'.tiff')
    #print ("[ GETTING BAND ]: ", band )
    srcband = g.GetRasterBand(band)
    if srcband is None:
        continue
    data1 = BandReadAsArray(srcband).astype(np.float)
    print(data1)
   # for i in range(3,sheet.nrows):
    b=sheet.cell_value(band+2,1)
    #print(b)
    dataOut = (1/data1)
    driver = gdal.GetDriverByName("ENVI")
    dsOut = driver.Create(outFile, g.RasterXSize, g.RasterYSize, 1)
    CopyDatasetInfo(g,dsOut)
    bandOut=dsOut.GetRasterBand(1)
    BandWriteArray(bandOut, dataOut)

Für die habe print(data1)ich nur eine "1", aber die realen Werte sind einige Schwimmer

0
[[1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 ...
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]]
1
[[1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 ...
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]]
2

Pixelwert 0,139200

PLZ helfen, den Fehler zu finden

Kais Tounsi
quelle