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
我认为您可以按照以下方法生成表示光栅单元的地理坐标光栅:https://gis.stackexchange.com/questions/177061/ascii-file-with-latitude-longitude-and-data-to-geotiff-using-python
然后,如果创建了纬度和经度点的形状文件,则可以使用以下方法获取每个点的光栅单元ID:
相关问题 更多 >
编程相关推荐