查找最近的地理栅格元素

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