0

蓝光：：

我在使用rasterstats包计算带有旋转数组的区域统计数据时遇到问题。我猜问题出在我的仿射矩阵上，但我不完全确定。下面是仿射变换矩阵和输出：

| 951.79, 0.45, 2999993.57|
| 0.00,-996.15,-1985797.84|
| 0.00, 0.00, 1.00|

在此处输入图像描述

背景：

我正在为地下水流模型创建文件，并且需要使用来自波多黎各农业水资源管理门户网站的一些 .csv 数据来计算每个模型网格单元的区域统计数据。这些数据可用于许多参数（例如 ET、tmax、tmin、precip 等）的每日时间步长。这些文件没有地理参考，但是可以使用辅助文件来指定每个单元格的经度/纬度，然后可以使用以下方式进行投影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()

在此处输入图像描述

4

0 回答 0