我正在研究整个景观中大量野生传粉媒介的计算模型。该模型本身已经完成,现在我正在努力进行后处理步骤。



我有一个OGR形状文件,表示景观上的样本位置:我正在尝试对这些点下的像素进行一些分析,但要这样做,我需要能够提取出一个点下的像素的值。通过Python仅使用OGR和GDAL的点下的像素数?我宁愿避免通过ReadAsArray()将整个栅格读入内存,因为我的输出栅格非常非常大(太大而无法容纳到内存中)。需要命令行调用。

评论

那么ReadAsArray()只读一点呢?因此,仅读取您感兴趣的单个单元格?您需要将点坐标转换为像素空间,并提取必要的像元。

查看gdalsrsinfo的代码,它向您展示如何使用GDALInvertGeoTransform()以及如何在地理空间和像素空间之间切换:trac.osgeo.org/gdal/browser/trunk/gdal/apps/gdalsrsinfo.cpp

如果您不介意使用PostGIS,请参阅此内容。它非常快,只有1条SQL行。

如果遇到此问题并可以访问PostGIS数据库,请记住这一点!我没有解决这个特定的问题,因此下面的GDAL解决方案可以解决问题。不过谢谢!

@kyle我不知道情况是否已更改,但它看起来像是GDALInvGeoTransform没有反转,这是一个示例。

#1 楼

您可以使用gdal.Dataset或gdal.Band ReadRaster方法。请参阅GDAL和OGR API教程以及以下示例。 ReadRaster不使用/不需要numpy,返回值是原始二进制数据,需要使用标准python struct模块解压缩。

示例:

from osgeo import gdal,ogr
import struct

src_filename = '/tmp/test.tif'
shp_filename = '/tmp/test.shp'

src_ds=gdal.Open(src_filename) 
gt=src_ds.GetGeoTransform()
rb=src_ds.GetRasterBand(1)

ds=ogr.Open(shp_filename)
lyr=ds.GetLayer()
for feat in lyr:
    geom = feat.GetGeometryRef()
    mx,my=geom.GetX(), geom.GetY()  #coord in map units

    #Convert from map to pixel coordinates.
    #Only works for geotransforms with no rotation.
    px = int((mx - gt[0]) / gt[1]) #x pixel
    py = int((my - gt[3]) / gt[5]) #y pixel

    structval=rb.ReadRaster(px,py,1,1,buf_type=gdal.GDT_UInt16) #Assumes 16 bit int aka 'short'
    intval = struct.unpack('h' , structval) #use the 'short' format code (2 bytes) not int (4 bytes)

    print intval[0] #intval is a tuple, length=1 as we only asked for 1 pixel value


或者,由于您不使用numpy的原因是避免使用ReadAsArray()读取整个数组,因此下面是使用numpy而不读取整个栅格的示例。

from osgeo import gdal,ogr
import struct

src_filename = '/tmp/test.tif'
shp_filename = '/tmp/test.shp'

src_ds=gdal.Open(src_filename) 
gt=src_ds.GetGeoTransform()
rb=src_ds.GetRasterBand(1)

ds=ogr.Open(shp_filename)
lyr=ds.GetLayer()
for feat in lyr:
    geom = feat.GetGeometryRef()
    mx,my=geom.GetX(), geom.GetY()  #coord in map units

    #Convert from map to pixel coordinates.
    #Only works for geotransforms with no rotation.
    px = int((mx - gt[0]) / gt[1]) #x pixel
    py = int((my - gt[3]) / gt[5]) #y pixel

    intval=rb.ReadAsArray(px,py,1,1)
    print intval[0] #intval is a numpy array, length=1 as we only asked for 1 pixel value


评论


以及如何保存为csv,表或其他对象?具有相同坐标长度的对象

–gonzalez.ivan90
18年1月26日在18:12

这是问题,卢克。谢谢! gis.stackexchange.com/questions/269603/…

–gonzalez.ivan90
18年1月29日在18:24

如果mx / my位于rb的边界之外,则分配px / py的行是错误的,因为int(-0.5)==0。您需要floor(...),然后需要检查两者px / py的值小于零(或在调用int()之前执行),因为负索引有效(它们获取数组的另一侧)。我很想知道是否有一种更整洁的方法来解决这个问题。另外,如何重新编写这些行,以便它们正确处理旋转?

–naught101
19-3-7在3:55