查找最近的地理栅格元素

Posted 2023-02-23

【中文标题】查找最近的地理栅格元素

【英文标题】：Finding nearest geographic raster element

【发布时间】：2017-05-04 12:34:29

【问题描述】：

我有一个大的经度纬度点列表，想在给定的地理坐标栅格中找到最近的矩形（因此哪个矩形包含该点）。

但是，对于栅格，我只有栅格中每个矩形（多边形）的质心。我知道矩形的大小为 250m x 250m。

仅检查到中心的绝对距离或地理距离是行不通的，因为矩形不一定是对齐的。我很高兴得到想法。

【参考方案1】：

我认为您可以按照这种方法生成代表栅格单元的地理坐标栅格：https://gis.stackexchange.com/questions/177061/ascii-file-with-latitude-longitude-and-data-to-geotiff-using-python

然后，如果您创建了纬度和经度点的 shapefile，您可以使用这种方法获取每个点的栅格单元 ID：

def GetRasterValueAtPoints(rasterfile, shapefile, fieldname):
'''
__author__ =   "Marc Weber <[e-mail]>"
Original code attribution: https://gis.stackexchange.com/a/46898/2856
returns raster values at points in a point shapefile
assumes same projection in shapefile and raster file
Arguments
---------
rasterfile        : a raster file with full pathname and extension
shapefile         : a shapefile with full pathname and extension
fieldname         : field name in the shapefile to identify values
'''
src_ds=gdal.Open(rasterfile)
no_data = src_ds.GetRasterBand(1).GetNoDataValue()
gt=src_ds.GetGeoTransform()
rb=src_ds.GetRasterBand(1)
df = pd.DataFrame(columns=(fieldname, "RasterVal"))
i = 0
ds=ogr.Open(shapefile)
lyr=ds.GetLayer()

for feat in lyr:
    geom = feat.GetGeometryRef()
    name = feat.GetField(fieldname)
    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)
    if intval == no_data:
        intval = -9999
    df.set_value(i,fieldname,name) 
    df.set_value(i,"RasterVal",float(intval)) 
    i+=1
return df