我正在尝试获取一个人口光栅，并对其进行重新采样和重新投影，以匹配降水光栅的形状和分辨率。
数据链接：人口数据：https://figshare.com/ndownloader/files/10257111Precipitation数据：https://www.ncei.noaa.gov/data/nclimgrid-monthly/access/nclimgrid_prcp.nc
人口数据是覆盖美国大陆的5个不同人口模型的每十年一系列栅格。如果你只需选择其中一个栅格，我就可以计算出其余的栅格（无论如何，我已经组合成一个多频带栅格）。例如，如果使用pop_m4_2010光栅将有所帮助。分辨率为1x1km，投影为阿尔伯斯等面积圆锥NAD 83 ESRI:102003。
降水数据是涵盖美国大陆月降水数据的netcdf文件。分辨率为5x5km，投影为WGS84 EPSG:4326。
我使用以下代码将netcdf转换为tiff：
import xarray as xr
import rioxarray as rio
prcp_file = xr.open_dataset('nclimgrid_prcp.nc')
prp = prcp_file['prcp']
prp = prp.rio.set_spatial_dims(x_dim='lon', y_dim='lat')
prp.rio.write_crs("epsg:4326", inplace=True)
prp.rio.to_raster('prp_raster.tiff')
我还使用QGIS打开填充文件（添加光栅层，导航到pop_m4_2010的下载文件夹并选择“w001001.adf”文件）。当我在WGS84项目中这样做时，QGIS会自动显示为强制重新投影，但我是新手，所以我不确定它是否正确。
从这一点开始，我尝试了几种方法来重新采样人口光栅，以匹配降水光栅的5x5分辨率。
- 在QGIS处理工具箱GRASS r.resample
- 中QGIS处理工具包Raster Layer Zonal Statistics
- 在Python中，老实说，我已经忘记了我在GDAL.Warp、Rasterio.Warp、仿射变换、rio.reproject_match等上关注的所有不同论坛帖子和教程。下面是一些代码尝试的示例
其中许多似乎都有效（尤其是rio.reproject_match似乎简单有效。然而，这些似乎都没有达到预期效果。当我通过传递县矢量形状文件的分区统计数据来测试生成的人口栅格的准确性时，该区域的人口总数要么为0，要么非常不准确。
我做错了什么？
重新投影匹配（_M）：
import rioxarray # for the extension to load
import xarray
import matplotlib.pyplot as plt
%matplotlib inline
def print_raster(raster):
print(
f"shape: {raster.rio.shape}\n"
f"resolution: {raster.rio.resolution()}\n"
f"bounds: {raster.rio.bounds()}\n"
f"sum: {raster.sum().item()}\n"
f"CRS: {raster.rio.crs}\n"
)
xds = rioxarray.open_rasterio('pop_m4_2010.tif')
xds_match = rioxarray.open_rasterio('prp_raster.tiff')
fig, axes = plt.subplots(ncols=2, figsize=(12,4))
xds.plot(ax=axes[0])
xds_match.plot(ax=axes[1])
plt.draw()
print("Original Raster:\n----------------\n")
print_raster(xds)
print("Raster to Match:\n----------------\n")
print_raster(xds_match)
xds_repr_match = xds.rio.reproject_match(xds_match)
print("Reprojected Raster:\n-------------------\n")
print_raster(xds_repr_match)
print("Raster to Match:\n----------------\n")
print_raster(xds_match)
xds_repr_match.rio.to_raster("reproj_pop.tif")
Rasterio的另一种方式。扭曲：
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
#open source raster
srcRst =rasterio.open('pop_m4_2010.tif')
print("source raster crs:")
print(srcRst.crs)
dstCrs = {'init': 'EPSG:4326'}
print("destination raster crs:")
print(dstCrs)
#calculate transform array and shape of reprojected raster
transform, width, height = calculate_default_transform(
srcRst.crs, dstCrs, srcRst.width, srcRst.height, *srcRst.bounds)
print("transform array of source raster")
print(srcRst.transform)
print("transform array of destination raster")
print(transform)
#working of the meta for the destination raster
kwargs = srcRst.meta.copy()
kwargs.update({
'crs': dstCrs,
'transform': transform,
'width': width,
'height': height
})
#open destination raster
dstRst = rasterio.open('pop_m4_2010_reproj4326.tif', 'w', **kwargs)
#reproject and save raster band data
for i in range(1, srcRst.count + 1):
reproject(
source=rasterio.band(srcRst, i),
destination=rasterio.band(dstRst, i),
#src_transform=srcRst.transform,
src_crs=srcRst.crs,
#dst_transform=transform,
dst_crs=dstCrs,
resampling=Resampling.bilinear)
print(i)
#close destination raster
dstRst.close()
这里是Rasterio的第二次尝试。Warp：
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
prcp = rasterio.open('prp_raster.tiff', mode = 'r')
with rasterio.open('pop_m4_2010.tif') as dataset:
# resample data to target shape
data = dataset.read(out_shape=(dataset.count,prcp.height,prcp.width), resampling=Resampling.bilinear)
# scale image transform
transform = dataset.transform * dataset.transform.scale((dataset.width / data.shape[-1]),
(dataset.height / data.shape[-2]))
# Register GDAL format drivers and configuration options with a
# context manager.
with rasterio.Env():
profile = src.profile
profile.update(
dtype=rasterio.float32,
count=1,
compress='lzw')
with rasterio.open('pop_m4_2010_resampledtoprcp.tif', 'w', **profile) as dst:
dst.write(data.astype(rasterio.float32))