Запись массива в растровый файл

30

Я новичок в ГИС.

У меня есть код, который преобразует инфракрасные изображения Марса в карты тепловой инерции, которые затем сохраняются в виде двумерных массивов. Я сохранял эти карты в виде файлов hdf5, но мне бы очень хотелось сохранить их как растровые изображения, чтобы я мог обработать их в QGIS. Я прошел несколько поисков, чтобы найти, как это сделать, но безуспешно. Я попытался следовать инструкциям в руководстве по адресу http://www.gis.usu.edu/~chrisg/python/, но файлы, которые я создаю с использованием его примера, открываются в виде простых серых рамок, когда я импортирую их в QGIS. Я чувствую, что если бы кто-то мог предложить простейшую возможную процедуру для упрощенного примера того, что я хотел бы сделать, тогда я мог бы добиться некоторого прогресса. У меня есть QGIS и GDAL, я был бы очень рад установить другие фреймворки, которые кто-нибудь мог бы порекомендовать. Я использую Mac OS 10.7.

Так что, если, например, у меня есть массив тепловых инерций, который выглядит следующим образом:

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

И для каждого пикселя у меня есть широта и долгота:

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

Какую процедуру люди рекомендуют преобразовать эти данные в растровый файл, который я могу открыть в QGIS?

EddyTheB
источник
На какой слайд учебника вы ссылаетесь?
РК

Ответы:

23

Одно из возможных решений вашей проблемы: преобразовать его в ASCII Raster, документация по которому находится здесь . Это должно быть довольно легко сделать с Python.

Итак, с вашими данными примера выше вы получите в файле .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

Это успешно добавляет как в QGIS, так и в ArcGIS, а стилизовано в ArcGIS и выглядит так: растровая версия выше

Приложение: Хотя вы можете добавить его в QGIS, как уже отмечалось, если вы попытаетесь зайти в свойства для него (чтобы стилизовать его), QGIS 1.8.0 зависает. Я собираюсь сообщить об этом как об ошибке. Если это случится и с вами, есть множество других бесплатных ГИС.

ГИС-Jonathan
источник
Это фантастика, спасибо. И я представляю, что, написав мой массив как файл ascii, я мог бы преобразовать его в двоичный формат, используя предварительно написанную функцию преобразования.
EddyTheB
К вашему сведению, у меня не было проблемы с зависанием QGIS, у меня тоже версия 1.8.0.
EddyTheB
31

Ниже приведен пример, который я написал для семинара, в котором используются модули Python numpy и gdal. Он считывает данные из одного файла .tif в пустой массив, выполняет переклассификацию значений в массиве и затем записывает их обратно в .tif.

Из вашего объяснения может показаться, что вам удалось записать действительный файл, но вам просто нужно обозначить его в QGIS. Если я правильно помню, когда вы впервые добавляете растр, он часто показывает все один цвет, если у вас нет уже существующей карты цветов.

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
источник
1
+1 за смывание - ударился головой о стену, пытаясь понять, как «спасти» вещь!
Badgley
Я должен был добавить, outDs = Noneчтобы сохранить это
JaakL
23

Я наконец-то нашел это решение, которое я получил из этого обсуждения ( http://osgeo-org.1560.n6.nabble.com/gdal-dev-numpy-array-to-raster-td4354924.html ). Мне это нравится, потому что я могу перейти от простого массива к растровому файлу TIF. Буду очень признателен за комментарии, которые могут улучшить решение. Я опубликую это здесь в случае, если кто-то еще ищет подобный ответ.

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
источник
3
«Вращение вдвое больше», чтобы учесть влияние повернутого бита у на х и повернутого бита х на у. См. Lists.osgeo.org/pipermail/gdal-dev/2011-July/029449.html, который пытается объяснить взаимосвязь между параметрами «вращения».
Дэйв Икс
Этот пост действительно полезен, спасибо. В моем случае, однако, я получаю TIF-файл, который полностью черный, когда я открываю его как изображение вне ArcGIS. Моя пространственная привязка - Британская национальная сеть (EPSG = 27700), а единицы измерения - метры.
FaCoffee
Я разместил вопрос здесь: gis.stackexchange.com/questions/232301/…
FaCoffee
Вы узнали, как установить кодировку IAU2000: 49900 Mars?
К.-Майкл Ай
4

Есть также хорошее решение в официальной GDAL / OGR Cookbook для Python.

Этот рецепт создает растр из массива

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)
Адам Эриксон
источник
Этот рецепт хорош, но есть проблема с финальным файлом TIFF. Значения пикселей в пикселях неправильные.
Shubham_geo
Возможно, вы видите странные несовместимости между ESRI WKT и OGC WKT: gis.stackexchange.com/questions/129764/…
Адам Эриксон
Одна вещь, с которой я столкнулся, это то, что упомянутый вами способ с легкостью изменит массив на растр. Но нам нужно установить географическую привязку этого растра с верхним левым и нижним правым координатами, используя gdal_translate. Один из способов сделать это - выполнить следующие два шага: 1) Сначала найти значения верхнего левого и нижнего правого лат-лона с помощью gdalinfo 2) Затем с помощью gdal_translate использовать геотиф (сгенерированный с помощью вышеупомянутого подхода преобразования массива в растр) привязать его к координатам левого верхнего и нижнего правого угла.
Shubham_geo
0

Альтернативой подходу, предложенному в других ответах, является использование rasterioпакета. У меня были проблемы с их созданием, gdalи я нашел этот сайт полезным.

Предполагая, что у вас есть другой файл tif ( other_file.tif) и массив numpy ( numpy_array), который имеет то же разрешение и размер, что и этот файл, для меня этот подход сработал:

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)
Тим Уильямс
источник