Zapisywanie tablicy numpy do pliku rastrowego

30

Jestem nowy w GIS.

Mam kod, który konwertuje obrazy Marsa w podczerwieni na mapy bezwładności cieplnej, które są następnie zapisywane jako tablice dwuwymiarowe. Zapisuję te mapy jako pliki hdf5, ale naprawdę chciałbym je zapisać jako obrazy rastrowe, aby móc je przetwarzać w QGIS. Przeszukałem wiele wyszukiwań, aby dowiedzieć się, jak to zrobić, ale bez powodzenia. Próbowałem postępować zgodnie z instrukcjami w samouczku na stronie http://www.gis.usu.edu/~chrisg/python/, ale pliki, które tworzę przy użyciu jego przykładowego kodu, otwierają się jako zwykłe szare pola podczas importowania ich do QGIS. Wydaje mi się, że gdyby ktoś mógł zasugerować najprostszą możliwą procedurę uproszczonemu przykładowi tego, co chciałbym zrobić, być może będę w stanie zrobić pewne postępy. Mam QGIS i GDAL, chętnie zainstaluję inne frameworki, które każdy może polecić. Używam Mac OS 10.7.

Więc jeśli na przykład mam tablicę liczb bezwładności, która wygląda następująco:

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) )

I dla każdego piksela mam szerokość i długość geograficzną:

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) ) 

Jaką procedurę poleciliby ludzie przekonwertować te dane na plik rastrowy, który mogę otworzyć w QGIS?

EddyTheB
źródło
Do którego slajdu w samouczku chodzi?
RK

Odpowiedzi:

23

Jedno z możliwych rozwiązań Twojego problemu: Przekształć go w raster ASCII, którego dokumentacja jest tutaj . Powinno to być dość łatwe w przypadku Pythona.

Tak więc z powyższymi przykładowymi danymi, w pliku .asc uzyskasz następujące elementy:

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

Z powodzeniem dodaje to zarówno do QGIS, jak i ArcGIS, a stylizowany w ArcGIS wygląda następująco: wersja rastrowa powyżej

Dodatek: Chociaż możesz dodać go do QGIS, jak wspomniano, jeśli spróbujesz przejść do jego właściwości (aby go stylizować), QGIS 1.8.0 zawiesza się. Mam zamiar zgłosić to jako błąd. Jeśli ci się to przydarzy, istnieje wiele innych bezpłatnych GIS.

GIS-Jonathan
źródło
To fantastycznie, dzięki. Wyobrażam sobie, że po zapisaniu mojej tablicy jako pliku ascii mógłbym ją przekonwertować na format binarny za pomocą wcześniej napisanej funkcji konwersji.
EddyTheB,
Do twojej wiadomości, nie miałem problemu z zawieszeniem się w QGIS, mam również wersję 1.8.0.
EddyTheB,
31

Poniżej znajduje się przykład, który napisałem na warsztaty z wykorzystaniem modułów numpy i gdal Python. Odczytuje dane z jednego pliku .tif do tablicy numpy, dokonuje przeklasyfikowania wartości w tablicy, a następnie zapisuje je z powrotem do pliku .tif.

Z twojego wyjaśnienia wygląda na to, że udało ci się napisać prawidłowy plik, ale wystarczy go symbolizować w QGIS. Jeśli dobrze pamiętam, po pierwszym dodaniu rastra często wyświetla się jeden kolor, jeśli nie masz wcześniej istniejącej mapy kolorów.

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
źródło
1
+1 za spłukiwanie - waliłem głową w ścianę, próbując wymyślić, jak „uratować” rzecz!
badgley,
Musiałem dodać, outDs = Noneaby go uratować
JaakL
23

W końcu wpadłem na to rozwiązanie, które zyskałem podczas tej dyskusji ( http://osgeo-org.1560.n6.nabble.com/gdal-dev-numpy-array-to-raster-td4354924.html ). Podoba mi się, ponieważ mogę przejść prosto z tablicy numpy do pliku rastrowego tif. Byłbym bardzo wdzięczny za komentarze, które mogłyby ulepszyć rozwiązanie. Zamieszczę go tutaj, na wypadek gdyby ktoś szukał podobnej odpowiedzi.

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
źródło
3
„Obrót jest dwa razy”, aby uwzględnić wpływ obróconego bitu y na x i obróconego bitu x na y. Widzieć lists.osgeo.org/pipermail/gdal-dev/2011-July/029449.html, który próbuje wyjaśnić wzajemne powiązania między parametrami „rotacji”.
Dave X
Ten post jest naprawdę przydatny, dzięki. Jednak w moim przypadku otrzymuję plik tif, który jest całkowicie czarny, gdy otwieram go jako obraz poza ArcGIS. Moje odniesienie przestrzenne to British National Grid (EPSG = 27700), a jednostki to metry.
FaCoffee,
Wysłałem
FaCoffee
Czy wiesz, jak ustawić kodowanie Mars IAU2000: 49900?
K.-Michael Aye
4

Istnieje również dobre rozwiązanie w oficjalnej książce kucharskiej GDAL / OGR dla Pythona.

Ten przepis tworzy raster z tablicy

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
źródło
Ten przepis jest dobry, ale występuje problem z ostatecznym plikiem tiff. Wartości długości lat w pikselach są nieprawidłowe.
Shubham_geo
Możesz zobaczyć dziwne niezgodności między ESRI WKT a OGC WKT: gis.stackexchange.com/questions/129764/…
Adam Erickson
Jedną rzeczą, na którą się natknąłem, jest sposób, o którym wspominasz, z łatwością zmieni tablicę na raster. Ale musimy georeferencję tego rastra ze skoordynowanym górnym lewym i dolnym prawym przy użyciu gdal_translate. Jednym ze sposobów na to jest wykonanie dwóch kroków: 1) Najpierw znajdź górne lewe i prawe dolne lew-lon poprzez gdalinfo 2) Następnie, poprzez gdal_translate wykorzystaj geotiff (wygenerowany z wyżej wspomnianego podejścia konwersji tablicy do rastra) do georeferencji za pomocą współrzędnych lewy górny lewy i dolny prawy dolny.
Shubham_geo
0

Alternatywą dla podejścia sugerowanego w innych odpowiedziach jest użycie rasteriopakietu. Miałem problemy z ich generowaniem gdali uznałem tę stronę za przydatną.

Zakładając, że masz inny plik tif ( other_file.tif) i tablicę numpy ( numpy_array), która ma taką samą rozdzielczość i zasięg jak ten plik, to podejście działało dla mnie:

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
źródło