查找最近的地理栅格元素

Question

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

但是，对于光栅，我只有光栅中每个矩形（多边形）的质心。我知道矩形的尺寸是250米x 250米

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

Answer 1

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

然后，如果创建了纬度和经度点的形状文件，则可以使用以下方法获取每个点的光栅单元ID：

def GetRasterValueAtPoints(rasterfile, shapefile, fieldname):
'''
__author__ =   "Marc Weber <weber.marc@epa.gov>"
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