查找最近的地理柵格元素

Question

我有一個經度緯度點的大列表，並且想要在給定的地理座標柵格中找到最近的矩形（以便哪個矩形包含點）。

但是，對於柵格，我只有柵格中每個矩形（多邊形）的質心。我知道矩形的尺寸是250米x 250米。

只檢查到中心的絕對距離或地理距離不起作用，因爲矩形不一定對齊。我很高興得到想法。

Answer 1

我想你可以生成你的代表採取這一方式的柵格單元的地理座標的柵格數據：https://gis.stackexchange.com/questions/177061/ascii-file-with-latitude-longitude-and-data-to-geotiff-using-python

，然後如果你創建你的latitute和經度點shape文件，你可以使用獲得每個點光柵小區ID這種方法：

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