Sepenuhnya memuat raster ke array numpy?

26

Saya telah mencoba untuk memeriksa filter saya pada DEM raster untuk pengenalan pola dan selalu menghasilkan baris dan kolom terakhir yang hilang (seperti..20) . Saya sudah mencoba dengan PIL library, memuat gambar. Lalu dengan numpy. Outputnya sama.

Saya pikir, ada yang salah dengan loop saya, ketika memeriksa nilai dalam array (hanya memilih piksel dengan Identifikasi di ArcCatalog) saya menyadari bahwa nilai piksel tidak dimuat ke dalam array.

Jadi, cukup buka, masukkan ke dalam array dan simpan gambar dari 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)

Hasil dalam memotong baris dan kolom terakhir. Maaf, tidak dapat memposting gambar

Adakah yang bisa membantu untuk memahami mengapa? Dan menyarankan beberapa solusi?

EDIT:

Jadi, saya berhasil memuat raster kecil ke array numpy dengan bantuan teman-teman, tetapi ketika memiliki gambar yang lebih besar saya mulai mendapatkan kesalahan. Saya kira ini tentang batas array numpy, dan array secara otomatis dibentuk ulang atau bertaut seperti itu ... Jadi ex:

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

Intinya adalah saya tidak ingin membaca blok demi blok karena saya perlu penyaringan, beberapa kali dengan filter berbeda, ukuran berbeda .. Apakah ada pekerjaan di sekitar atau saya harus belajar rading dengan blok: O

nafsu makan
sumber

Jawaban:

42

jika Anda memiliki binding python-gdal:

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

Dan Anda selesai:

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)
Nickes
sumber
Ya, dengan gdal kurasa aku tidak punya masalah, tapi aku mencoba menggunakan perpustakaan yang lebih sedikit ... Dan numpy sepertinya sangat populer untuk itu 'sambil googling'. Tahu, memang, mengapa numpy / PIL berhenti memuat ???
najuste
Saya tidak tahu PIL harus cukup kuat sehingga dikirimkan dengan python. Tapi imho geotiff lebih dari gambar - mereka membawa banyak metadata misalnya - dan PIL bukan (lagi-lagi imho) alat yang tepat.
Nickes
Saya kadang-kadang benci persyaratan kutip dan slash yang berbeda itu, ketika membuka data .. Tapi bagaimana dengan menulis array numpy kembali ke Raster? Ia bekerja dengan PIL library, tetapi menggunakan outputRaster.GetRasterBand (1) .WriteArray (myarray) menghasilkan raster yang tidak valid ..
najuste
jangan lupa untuk menyiram data ke disk, dengan outBand.FlushCache (). Anda dapat menemukan beberapa tutorial di sini: gis.usu.edu/~chrisg/python/2009
nickves
1
Periksa " lists.osgeo.org/pipermail/gdal-dev/2010-January/023309.html " - sepertinya Anda kehabisan atau ram.
nickes
21

Anda dapat menggunakan rasterio untuk berinteraksi dengan array NumPy. Untuk membaca raster ke 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]

Ini akan membaca semuanya menjadi array numpy 3D arr, dengan dimensi [band, row, col].


Berikut adalah contoh lanjutan untuk membaca, mengedit piksel, lalu menyimpannya kembali ke 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)

Raster akan ditulis dan ditutup pada akhir pernyataan "dengan" .

Mike T
sumber
Mengapa kita tidak bisa melihat semua nilai ketika saya menulis print (arr). Ini memisahkan nilai dengan ini ..., ...,?
Mustafa Uçar
@ MustafaUçar, inilah cara NumPy mencetak array, yang dapat Anda modifikasi . Atau potong jendela array untuk dicetak, di antara banyak trik Numpy lainnya.
Mike T
Pertanyaan umum. Jika saya ingin menampilkan array tunggal dengan beberapa adegan, menjadi empat dimensi seperti (adegan, tinggi, lebar, pita), bagaimana saya harus memodifikasi cuplikan ini?
Ricardo Barros Lourenço
@ RicardoBarrosLourenço Saya menduga dimensi keempat Anda (adegan?) Disimpan di setiap file. Saya pertama-tama akan mengisi array numpy 4D kosong, kemudian loop melalui setiap file (adegan) dan masukkan bagian 3D masing-masing. Anda mungkin perlu arr.transpose((1, 2, 0))mendapatkan (tinggi, lebar, pita) dari setiap file.
Mike T
@ MikeT populasi ini akan seperti np.append()?
Ricardo Barros Lourenço
3

Memang saya membaca gambar png tua biasa, tetapi ini bekerja menggunakan scipy ( imsavemenggunakan 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)

Png saya yang dihasilkan juga 81 x 90 piksel.

Chad Cooper
sumber
Terima kasih, tapi saya mencoba menggunakan lebih sedikit perpustakaan .. Dan untuk saat ini saya dapat membuatnya dengan gdal + numpy ... (semoga tanpa PIL).
najuste
1
@najuste OS apa yang aktif? Mac dan kebanyakan rasa Linux datang dengan scipydan numpy.
Chad Cooper
Rupanya ... saya di Windows, berbagai versi Win. : /
najuste
2

Solusi saya menggunakan gdal terlihat seperti ini. Saya pikir ini sangat dapat digunakan kembali.

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
sumber
0

saya menggunakan gambar hyperspectral dengan 158 band. Saya ingin menghitung raster. tapi saya mengerti

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)

untuk print(data1)saya hanya mendapat beberapa "1", tetapi nilai sebenarnya adalah beberapa mengapung

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

Nilai piksel 0,139200

Tolong bantu untuk menemukan kesalahan

Kais Tounsi
sumber