我是GIS的新手。

我有一些代码将火星的红外图像转换成热惯性图，然后将其存储为2D numpy数组。我一直将这些地图另存为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中打开的栅格文件？

评论

您是指教程中的哪张幻灯片？
我喜欢Python，但在下面看到如此多的答案，这些答案将成为R中的衬里。
有没有人用AAIGrid选项在这里测试过Create方法？甚至AIG。它不会创建output_raster对象。也没有引发任何错误。元数据看起来可以创建Create对象，但是我总是得到NoneType。 GTiff可以正常工作。这是我尝试的代码段：```xmin，ymin，xmax，ymax = [lons.min（）。values，lats.min（）。values，\ lons.max（）。values ，lats.max（）。values] ncols = lons.shape [0] nrows = lats.shape [0] xres =（xmax-xmin）/ float（ncols）yres =（ymax-ymin）/ float（nrows）geotransform = [xmin，xres，0，ymin，0，yres] driver = gdal.GetDriverByName（'AAIGrid'）元数据= driver.GetMetadata（）outp

#1 楼

一种可能的解决方案：将其转换为ASCII栅格，其文档在这里。使用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将挂起。我将其报告为错误。如果您也遇到这种情况，那么那里还有许多其他免费的GIS。

评论

太好了，谢谢。而且我想像已将数组写为ascii文件，然后可以使用预先编写的转换函数将其转换为二进制格式。

– EddyTheB
2012年10月22日19:03
仅供参考，我没有QGIS的棘手问题，我也有版本1.8.0。

– EddyTheB
2012年10月22日19:10
嗨，有没有什么链接可以找到如何将numpy数组转换为ASCII的链接，以便最终文件可以得到geotiff进行绘制。

–sameer_nubia
20年11月2日，10：54


#2 楼

下面是我为一个使用numpy和gdal Python模块的研讨会编写的示例。它从一个.tif文件读取数据到一个numpy数组中，对该数组中的值进行重新分类，然后将其写回到.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


评论

+1表示冲洗-正在将我的头撞在墙上，试图弄清楚如何“保存”东西！

– Badgley
2015年8月7日在16:46
我必须添加outDs = None才能保存它

– JaakL
17年11月20日在11:44

#3 楼

我终于找到了这个解决方案，该解决方案是我从这次讨论中获得的（http://osgeo-org.1560.n6.nabble.com/gdal-dev-numpy-array-to-raster-td4354924.html）。我喜欢它，因为我可以直接从一个numpy数组转到一个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()


评论

“旋转两次”考虑了x上y的旋转位和y上x的旋转位的影响。请参阅lists.osgeo.org/pipermail/gdal-dev/2011-July/029449.html，它试图解释“旋转”参数之间的相互关系。

–戴夫X
2014年8月11日下午4:52
这篇文章真的很有用，谢谢。但是，就我而言，当我在ArcGIS外部将其作为图像打开时，我得到的是一个全黑色的tif文件。我的空间参考是英国国家网格（EPSG = 27700），单位是米。

–FaCoffee
17年3月16日在9:13
我在这里发布了一个问题：gis.stackexchange.com/questions/232301/…

–FaCoffee
17年3月16日在9:39
您是否了解如何设置IAU2000：49900火星编码？

– K.-Michael Aye
19年3月18日在7:38

#4 楼

适用于Python的GDAL / OGR官方食谱中也有一个不错的解决方案。


此配方根据数组创建栅格。


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文件有问题。像素的Latlon值不正确。

–Shubham_geo
19年8月16日在13:50
您可能会发现ESRI WKT与OGC WKT之间存在怪异的不兼容之处：gis.stackexchange.com/questions/129764/…

–亚当·埃里克森（Adam Erickson）
19年8月17日在17:13
我遇到的一件事是，您提到的方式肯定会轻松地将数组更改为栅格。但是我们需要使用gdal_translate协调左上角和右下角的栅格数据。一种方法是遵循以下两个步骤：1）首先通过gdalinfo查找左上和右下纬度值2）然后，通过gdal_translate利用geotiff（通过上述将数组转换为栅格的方法生成）使用左上和右下经纬度坐标对其进行地理配准。

–Shubham_geo
19年8月18日在7:59


#5 楼

其他答案中建议的方法的替代方法是使用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)