BLUF:
У меня проблемы с вычислением зональной статистики с помощью повернутого массива с использованием пакета rasterstats . Я предполагаю, что проблема в моей аффинной матрице, но я не совсем уверен. Ниже приведена матрица аффинного преобразования и вывод:
| 951.79, 0.45, 2999993.57|
| 0.00,-996.15,-1985797.84|
| 0.00, 0.00, 1.00|
Справка:
Я создаю файлы для модели потока подземных вод, и мне нужно вычислить зональную статистику для каждой ячейки сетки модели, используя некоторые данные .csv с веб-портала управления водными ресурсами в Пуэрто-Рико. Эти данные доступны на ежедневном временном шаге для многочисленных параметров (например, ET, tmax, tmin, осадок и т. Д.). Эти файлы не имеют географической привязки, но доступны вспомогательные файлы, в которых указывается lon / lat для каждой ячейки, которую затем можно спроецировать с помощью pyproj:
import pandas as pd
import numpy as np
import pyproj
url_base = 'http://academic.uprm.edu/hdc/GOES-PRWEB_RESULTS'
# Load some data
f = '/'.join([url_base, 'actual_ET', 'actual_ET20090101.csv'])
array = pd.read_csv(f, header=None).values
# Read longitude values
f = '/'.join([url_base, 'NON_TRANSIENT_PARAMETERS', 'LONGITUDE.csv'])
lon = pd.read_csv(f, header=None).values
# Read latitude values
f = '/'.join([url_base, 'NON_TRANSIENT_PARAMETERS', 'LATITUDE.csv'])
lat = np.flipud(pd.read_csv(f, header=None).values)
# Project to x/y coordinates (North America Albers Equal Area Conic)
aea = pyproj.Proj('+init=ESRI:102008')
x, y = aea(lon, lat)
Прежде чем я смогу вычислить зональную статистику, мне нужно создать аффинное преобразование, которое связывает координаты строки / столбца с проецируемыми координатами x / y. Я задаю 6 параметров для создания объекта Affine, используя аффинную библиотеку:
import math
from affine import Affine
length_of_degree_longitude_at_lat_mean = 105754.71 # 18.25 degrees via http://www.csgnetwork.com/degreelenllavcalc.html
length_of_degree_latitude_at_lat_mean = 110683.25 # 18.25 degrees via http://www.csgnetwork.com/degreelenllavcalc.html
# Find the upper left x, y
xul, yul = aea(lon[0][0], lat[0][0])
xll, yll = aea(lon[-1][0], lat[-1][0])
xur, yur = aea(lon[0][-1], lat[0][-1])
xlr, ylr = aea(lon[-1][-1], lat[-1][-1])
# Compute pixel width
a = abs(lon[0][1] - lon[0][0]) * length_of_degree_longitude_at_lat_mean
# Row rotation
adj = abs(xlr - xll)
opp = abs(ylr - yll)
b = math.atan(opp/adj)
# x-coordinate of the upper left corner
c = xul - a / 2
# Compute pixel height
e = -abs(lat[1][0] - lat[0][0]) * length_of_degree_latitude_at_lat_mean
# Column rotation
d = 0
# y-coordinate of the upper left corner
f = yul - e / 2
affine = Affine(a, b, c, d, e, f)
где:
- a = ширина пикселя
- b = вращение строки (обычно ноль)
- c = x-координата левого верхнего угла левого верхнего пикселя
- d = вращение столбца (обычно ноль)
- e = высота пикселя (обычно отрицательная)
- f = y-координата верхнего левого угла верхнего левого пикселя
(из https://www.perrygeo.com/python-affine-transforms.html)
Полученная аффинная матрица выглядит разумно, и я попытался передать вращение строк и столбцов как радианы, так и градусы с небольшим изменением результата. Ссылка на сетку: grid_2km.geojson
import rasterstats
import matplotlib.pyplot as plt
grid_f = 'grid_2km.geojson'
gdf = gpd.read_file(grid_f)
zs = rasterstats.zonal_stats(gdf,
array,
affine=affine,
stats=['mean'])
df = pd.DataFrame(zs).fillna(value=np.nan)
fig = plt.figure(figsize=(14, 6))
ax = fig.add_subplot(131, aspect='equal')
ax.pcolormesh(x, y, np.zeros_like(array))
ax.pcolormesh(x, y, array)
ax.set_title('Projected Data')
ax = fig.add_subplot(132, aspect='equal')
gdf.plot(ax=ax)
ax.set_title('Projected Shapefile')
ax = fig.add_subplot(133, aspect='equal')
ax.imshow(df['mean'].values.reshape((gdf.row.max(), gdf.col.max())))
ax.set_title('Zonal Statistics Output')
plt.tight_layout()
plt.show()
Кроме того, существует несоответствие между парами значений x, y, преобразованными с использованием аффинного объекта, и теми, которые получены из собственных значений lon, lat с использованием pyproj:
rr = np.array([np.ones(x.shape[1], dtype=np.int) * i for i in range(x.shape[0])])
cc = np.array([np.arange(x.shape[1]) for i in range(x.shape[0])])
x1, y1 = affine * (cc, rr)
fig = plt.figure(figsize=(14, 6))
ax = fig.add_subplot(111, aspect='equal')
ax.scatter(x, y, s=.2)
ax.scatter(x1, y1, s=.2)
plt.show()
