Menulis array numpy ke file raster

30

Saya baru mengenal GIS.

Saya memiliki beberapa kode yang mengubah gambar inframerah Mars menjadi peta inersia termal, yang kemudian disimpan sebagai array numpy 2D. Saya telah menyimpan peta-peta ini sebagai file hdf5 tetapi saya benar-benar ingin menyimpannya sebagai gambar raster sehingga saya dapat memprosesnya di QGIS. Saya telah melalui beberapa pencarian untuk menemukan cara melakukan ini tetapi tidak berhasil. Saya sudah mencoba mengikuti petunjuk dalam tutorial di http://www.gis.usu.edu/~chrisg/python/ tetapi file yang saya hasilkan menggunakan kode contohnya terbuka sebagai kotak abu-abu ketika saya mengimpornya ke QGIS. Saya merasa jika seseorang dapat menyarankan prosedur yang paling sederhana untuk contoh sederhana dari apa yang ingin saya lakukan maka saya mungkin dapat membuat beberapa kemajuan. Saya memiliki QGIS dan GDAL, saya akan sangat senang menginstal kerangka kerja lain yang bisa direkomendasikan oleh siapa pun. Saya menggunakan Mac OS 10.7.

Jadi jika misalnya saya memiliki array inertia termal yang numpy yang terlihat seperti:

TI = ( (0.1, 0.2, 0.3, 0.4),
       (0.2, 0.3, 0.4, 0.5),
       (0.3, 0.4, 0.5, 0.6),
       (0.4, 0.5, 0.6, 0.7) )

Dan untuk setiap piksel saya memiliki garis lintang dan bujur:

lat = ( (10.0, 10.0, 10.0, 10.0),
        ( 9.5,  9.5,  9.5,  9.5),
        ( 9.0,  9.0,  9.0,  9.0),
        ( 8.5,  8.5,  8.5,  8.5) )
lon = ( (20.0, 20.5, 21.0, 21.5),
        (20.0, 20.5, 21.0, 21.5),
        (20.0, 20.5, 21.0, 21.5),
        (20.0, 20.5, 21.0, 21.5) ) 

Prosedur mana yang akan direkomendasikan orang untuk mengubah data ini menjadi file raster yang dapat saya buka di QGIS?

EddyTheB
sumber
Slide mana pada tutorial yang Anda maksud?
RK

Jawaban:

23

Salah satu solusi yang mungkin untuk masalah Anda: Konversikan menjadi ASCII Raster, dokumentasi yang ada di sini . Ini seharusnya cukup mudah dilakukan dengan python.

Jadi dengan contoh data Anda di atas, Anda akan berakhir dengan yang berikut dalam file .asc:

ncols 4
nrows 4
xllcorner 20
yllcorner 8.5
cellsize 0.5
nodata_value -9999
0.1 0.2 0.3 0.4
0.2 0.3 0.4 0.5
0.3 0.4 0.5 0.6
0.4 0.5 0.6 0.7

Ini berhasil menambah QGIS dan ArcGIS, dan bergaya di ArcGIS seperti ini: versi raster di atas

Tambahan: Meskipun Anda dapat menambahkannya ke QGIS seperti yang disebutkan, jika Anda mencoba dan masuk ke properti untuk itu (untuk menyesuaikan dgn mode itu), QGIS 1.8.0 hang. Saya akan melaporkan itu sebagai bug. Jika ini terjadi pada Anda juga, maka ada banyak SIG gratis lainnya di luar sana.

GIS-Jonathan
sumber
Fantastis, terima kasih. Dan saya membayangkan bahwa setelah menulis array saya sebagai file ascii saya dapat mengubahnya menjadi format biner menggunakan fungsi konversi pra-tertulis.
EddyTheB
FYI, saya tidak memiliki masalah menggantung dengan QGIS, saya juga memiliki versi 1.8.0.
EddyTheB
31

Di bawah ini adalah contoh yang saya tulis untuk lokakarya yang memanfaatkan modul Python numpy dan gdal. Bunyinya data dari satu file .tif ke array numpy, melakukan klasifikasi ulang dari nilai-nilai dalam array dan kemudian menulisnya kembali ke .tif.

Dari penjelasan Anda, sepertinya Anda mungkin berhasil menulis file yang valid, tetapi Anda hanya perlu menyimbolkannya di QGIS. Jika saya ingat benar, ketika Anda pertama kali menambahkan raster, sering muncul semua satu warna jika Anda tidak memiliki peta warna yang sudah ada sebelumnya.

import numpy, sys
from osgeo import gdal
from osgeo.gdalconst import *


# register all of the GDAL drivers
gdal.AllRegister()

# open the image
inDs = gdal.Open("c:/workshop/examples/raster_reclass/data/cropland_40.tif")
if inDs is None:
  print 'Could not open image file'
  sys.exit(1)

# read in the crop data and get info about it
band1 = inDs.GetRasterBand(1)
rows = inDs.RasterYSize
cols = inDs.RasterXSize

cropData = band1.ReadAsArray(0,0,cols,rows)

listAg = [1,5,6,22,23,24,41,42,28,37]
listNotAg = [111,195,141,181,121,122,190,62]

# create the output image
driver = inDs.GetDriver()
#print driver
outDs = driver.Create("c:/workshop/examples/raster_reclass/output/reclass_40.tif", cols, rows, 1, GDT_Int32)
if outDs is None:
    print 'Could not create reclass_40.tif'
    sys.exit(1)

outBand = outDs.GetRasterBand(1)
outData = numpy.zeros((rows,cols), numpy.int16)


for i in range(0, rows):
    for j in range(0, cols):

    if cropData[i,j] in listAg:
        outData[i,j] = 100
    elif cropData[i,j] in listNotAg:
        outData[i,j] = -100
    else:
        outData[i,j] = 0


# write the data
outBand.WriteArray(outData, 0, 0)

# flush data to disk, set the NoData value and calculate stats
outBand.FlushCache()
outBand.SetNoDataValue(-99)

# georeference the image and set the projection
outDs.SetGeoTransform(inDs.GetGeoTransform())
outDs.SetProjection(inDs.GetProjection())

del outData
DavidF
sumber
1
+1 untuk disiram - membenturkan kepala ke dinding mencoba mencari cara untuk 'menyelamatkan' benda itu!
badgley
Saya harus menambahkan outDs = Noneuntuk menyelamatkannya
JaakL
23

Saya akhirnya menemukan solusi ini, yang saya dapatkan dari diskusi ini ( http://osgeo-org.1560.n6.nabble.com/gdal-dev-numpy-array-to-raster-td4354924.html ). Saya suka karena saya bisa langsung dari array numpy ke file raster tif. Saya akan sangat berterima kasih atas komentar yang dapat meningkatkan solusi. Saya akan memposting di sini kalau-kalau ada orang lain mencari jawaban yang sama.

import numpy as np
from osgeo import gdal
from osgeo import gdal_array
from osgeo import osr
import matplotlib.pylab as plt

array = np.array(( (0.1, 0.2, 0.3, 0.4),
                   (0.2, 0.3, 0.4, 0.5),
                   (0.3, 0.4, 0.5, 0.6),
                   (0.4, 0.5, 0.6, 0.7),
                   (0.5, 0.6, 0.7, 0.8) ))
# My image array      
lat = np.array(( (10.0, 10.0, 10.0, 10.0),
                 ( 9.5,  9.5,  9.5,  9.5),
                 ( 9.0,  9.0,  9.0,  9.0),
                 ( 8.5,  8.5,  8.5,  8.5),
                 ( 8.0,  8.0,  8.0,  8.0) ))
lon = np.array(( (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5) ))
# For each pixel I know it's latitude and longitude.
# As you'll see below you only really need the coordinates of
# one corner, and the resolution of the file.

xmin,ymin,xmax,ymax = [lon.min(),lat.min(),lon.max(),lat.max()]
nrows,ncols = np.shape(array)
xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)
geotransform=(xmin,xres,0,ymax,0, -yres)   
# That's (top left x, w-e pixel resolution, rotation (0 if North is up), 
#         top left y, rotation (0 if North is up), n-s pixel resolution)
# I don't know why rotation is in twice???

output_raster = gdal.GetDriverByName('GTiff').Create('myraster.tif',ncols, nrows, 1 ,gdal.GDT_Float32)  # Open the file
output_raster.SetGeoTransform(geotransform)  # Specify its coordinates
srs = osr.SpatialReference()                 # Establish its coordinate encoding
srs.ImportFromEPSG(4326)                     # This one specifies WGS84 lat long.
                                             # Anyone know how to specify the 
                                             # IAU2000:49900 Mars encoding?
output_raster.SetProjection( srs.ExportToWkt() )   # Exports the coordinate system 
                                                   # to the file
output_raster.GetRasterBand(1).WriteArray(array)   # Writes my array to the raster

output_raster.FlushCache()
EddyTheB
sumber
3
"Rotasi berada dalam dua kali" untuk menjelaskan efek bit yang diputar dari y pada x dan bit yang diputar dari x pada y. Lihat lists.osgeo.org/pipermail/gdal-dev/2011-July/029449.html yang mencoba menjelaskan keterkaitan antara parameter "rotasi".
Dave X
Posting ini sangat berguna, terima kasih. Namun, dalam kasus saya, saya mendapatkan file tif yang benar-benar hitam ketika saya membukanya sebagai gambar di luar ArcGIS. Referensi spasial saya adalah British National Grid (EPSG = 27700), dan unitnya adalah meter.
FaCoffee
Saya telah mengirim pertanyaan di sini: gis.stackexchange.com/questions/232301/…
FaCoffee
Apakah Anda mengetahui cara mengatur pengkodean IAU2000: 49900 Mars?
K.-Michael Aye
4

Ada juga solusi bagus di Cookbook GDAL / OGR resmi untuk Python.

Resep ini membuat raster dari array

import gdal, ogr, os, osr
import numpy as np


def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):

    cols = array.shape[1]
    rows = array.shape[0]
    originX = rasterOrigin[0]
    originY = rasterOrigin[1]

    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Byte)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()


def main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):
    reversed_arr = array[::-1] # reverse array so the tif looks like the array
    array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,reversed_arr) # convert array to raster


if __name__ == "__main__":
    rasterOrigin = (-123.25745,45.43013)
    pixelWidth = 10
    pixelHeight = 10
    newRasterfn = 'test.tif'
    array = np.array([[ 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, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1],
                      [ 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1],
                      [ 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1],
                      [ 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1],
                      [ 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 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]])


    main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array)
Adam Erickson
sumber
Resep ini bagus tetapi ada masalah dengan file tiff akhir. Nilai lat-lon piksel tidak tepat.
Shubham_geo
Anda mungkin melihat ketidakcocokan aneh antara ESRI WKT dan OGC WKT: gis.stackexchange.com/questions/129764/…
Adam Erickson
Satu hal yang saya temui adalah, cara yang Anda sebutkan pasti akan mengubah array menjadi raster dengan mudah. Tetapi kita perlu melakukan georeferensi raster ini dengan kiri atas dan kanan bawah terkoordinasi menggunakan gdal_translate. Salah satu cara untuk melakukannya adalah mengikuti dua langkah: 1) Pertama menemukan nilai lat-lon atas-kiri dan kanan bawah melalui gdalinfo 2) Kemudian, melalui gdal_translate gunakan geotiff (dihasilkan dengan pendekatan yang disebutkan di atas untuk mengubah array menjadi raster) untuk melakukan georeferensi dengan koordinat lat-lon kanan atas dan kiri.
Shubham_geo
0

Alternatif untuk pendekatan yang disarankan dalam jawaban lain adalah dengan menggunakan rasteriopaket. Saya mengalami masalah dalam menghasilkan penggunaan ini gdaldan menganggap situs ini bermanfaat.

Dengan asumsi Anda memiliki file tif lain ( other_file.tif) dan array numpy ( numpy_array) yang memiliki resolusi dan luas yang sama dengan file ini, ini adalah pendekatan yang bekerja untuk saya:

import rasterio as rio    

with rio.open('other_file.tif') as src:
    ras_data = src.read()
    ras_meta = src.profile

# make any necessary changes to raster properties, e.g.:
ras_meta['dtype'] = "int32"
ras_meta['nodata'] = -99

with rio.open('outname.tif', 'w', **ras_meta) as dst:
    dst.write(numpy_array, 1)
Tim Williams
sumber