Oto funkcja, którą napisałem używając Pythona i Gdal. Funkcja przyjmuje listę rastrów i ramkę danych pandy zawierającą współrzędne punktu i zwraca ramkę danych pandy ze współrzędnymi punktu, centroidami dla odpowiednich komórek rastrowych i odpowiednich wartości komórek. Ta funkcja jest częścią pakietu rozwijanego pakietu chorospy (znajdującego się tutaj ).
import pandas
import numpy
from osgeo import gdal
def getValuesAtPoint(indir, rasterfileList, pos, lon, lat):
#gt(2) and gt(4) coefficients are zero, and the gt(1) is pixel width, and gt(5) is pixel height.
#The (gt(0),gt(3)) position is the top left corner of the top left pixel of the raster.
for i, rs in enumerate(rasterfileList):
presValues = []
gdata = gdal.Open('{}/{}.tif'.format(indir,rs))
gt = gdata.GetGeoTransform()
band = gdata.GetRasterBand(1)
nodata = band.GetNoDataValue()
x0, y0 , w , h = gt[0], gt[3], gt[1], gt[5]
data = band.ReadAsArray().astype(numpy.float)
#free memory
del gdata
if i == 0:
#iterate through the points
for p in pos.iterrows():
x = int((p[1][lon] - x0)/w)
Xc = x0 + x*w + w/2 #the cell center x
y = int((p[1][lat] - y0)/h)
Yc = y0 + y*h + h/2 #the cell center y
try:
if data[y,x] != nodata:
presVAL = [p[1][lon],p[1][lat], '{:.6f}'.format(Xc), '{:.6f}'.format(Yc), data[y,x]]
presValues.append(presVAL)
except:
pass
df = pandas.DataFrame(presValues, columns=['x', 'y', 'Xc', 'Yc', rs])
else:
#iterate through the points
for p in pos.iterrows():
x = int((p[1][lon] - x0)/w)
y = int((p[1][lat] - y0)/h)
try:
if data[y,x] != nodata:
presValues.append(data[y,x])
except:
pass
df[rs] = pandas.Series(presValues)
del data, band
return df
Przykład tego, jak go uruchomić, biorąc pod uwagę, że rastry znajdują się w bieżącym katalogu roboczym:
rasDf = getValuesAtPoint('.', ['raster1', 'raster2'], inPoints, 'x', 'y')