问题描述
我一直在研究从栅格数据集到栅格数据集（即tiff文件）进行栅格重投影和重采样的不同方法。
尽管做出了努力，但我还是无法正确理解Rasterio webpage中的某些可用选项。
在下面的列表中，我已订购了感兴趣的主题以供以后讨论：
- 虚拟包装器重新投影
- Rasterio标准驱动程序重新投影
我还订购了一些小片段，为清单的每个主题提供示例，以供进一步讨论：
# Importing main libraries
import os,sys
import zipfile
import numpy as np
import glob
import concurrent.futures
import Rasterio
import affine
from Rasterio.crs import CRS
from Rasterio.enums import resampling
from Rasterio import shutil as rio_shutil
from Rasterio.vrt import WarpedVRT
from Rasterio.warp import reproject,calculate_default_transform
# Defining some helpful functions
def get_transform_res(transform):
return (transform.a,transform.e)
def get_bounds(filename):
with Rasterio.open(filename) as dst:
return dst.bounds
def get_crs(filepath):
if isinstance(filepath,Rasterio.DatasetReader):
dataset = filepath
else:
dataset = Rasterio.open(filepath)
kwargs = dataset.Meta.copy()
crs = kwargs['crs']
dataset.close()
return crs
def get_dim_sizes_from_ds(filepath):
if isinstance(filepath,Rasterio.DatasetReader):
ds = filepath
else:
ds = Rasterio.open(filepath)
height,width = ds.height,ds.width
ds.close()
return height,width
def check_ds_resolution(filepath):
if isinstance(filepath,Rasterio.DatasetReader):
dataset = filepath
else:
dataset = Rasterio.open(filepath)
kwargs = dataset.Meta.copy()
transform = kwargs['transform']
dataset.close()
return get_transform_res(transform)
def resample_and_save_via_vrt(filepaths,xres=None,yres=None,crs_from_epsg=None,directory = None,stack=True,stacked_filename = 'stacked.tif',windowing=False):
'''
Description:
Function that does
1) resamples
2) reprojects (if crs is provided),3) stacks multiple files into a single one: (if stack ==True)
4) exports files to user-defined CRS resolution.
'''
if not os.path.exists(directory):
os.makedirs(directory)
input_files = filepaths
# Destination CRS being taken from one of the datasets
if crs_from_epsg == None:
dst_crs = get_crs(filepaths[0])
else:
dst_crs = CRS.from_epsg(crs_from_epsg)
# standard bounds that are in CRS coordinate
origin_bounds = get_bounds( filepaths[0] )
# Output standard image dimensions
origin_height,origin_width = get_dim_sizes_from_ds(filepaths[0])
# Output image transform
left,bottom,right,top = origin_bounds
origin_xres = (right - left) / origin_width
origin_yres = (top - bottom) / origin_height
if xres == None or yres == None:
# same heights and widths of the original dataset
dst_height,dst_width = origin_height,origin_width
# standard destination transform
dst_transform = affine.Affine(origin_xres,0.0,left,-origin_yres,top)
else:
# ensuring that all resolutions are being passed correctly by the user
if xres == None:
xres = yres
elif yres == None:
yres == xres
else:
pass
dst_width = abs(int( (right - left) / xres ))
dst_height = abs(int( (top - bottom) / yres ))
# transform with the new width resolutions
dst_transform = affine.Affine(xres,-yres,top)
vrt_options = {
'resampling': resampling.cubic,'crs': dst_crs,'transform': dst_transform,'height': dst_height,'width': dst_width
}
print('Files being saved in: ',directory,'\n')
if stack is False:
for path in input_files:
with Rasterio.open(path) as src:
# https://Rasterio.readthedocs.io/en/latest/topics/virtual-warping.html
with WarpedVRT(src,**vrt_options) as vrt:
# At this point 'vrt' is a full dataset with dimensions,# CRS,and spatial extent matching 'vrt_options'.
name = os.path.basename(path).split('.')[0] + '.tif'
outfile = os.path.join(directory,name + '_resampled_{0}x_{1}_y_reprojected_to_epsg{2}'.format(xres,yres,dst_crs.to_epsg()))
# Read all data into memory.
if not windowing:
data = vrt.read()
rio_shutil.copy(vrt,outfile,driver='GTiff')
else:
# Process the dataset in chunks. Likely not very efficient.
# Dump the aligned data into a new file. A VRT representing
# this transformation can also be produced by switching
# to the VRT driver.
for _,window in vrt.block_windows():
data = vrt.read(window=window)
rio_shutil.copy(vrt,driver='GTiff',window=window)
print('{0}'.format(name),'\t\t is complete')
if stack == True:
# Adapted from Ref: https://gis.stackexchange.com/questions/223910/using-Rasterio-or-gdal-to-stack-multiple-bands-without-using-subprocess-commands
outfile = os.path.join(directory,stacked_filename)
vrt_options.update( driver = 'GTiff',transform = dst_transform,height = dst_height,width = dst_width,count = len(input_files)
)
with Rasterio.open(outfile,'w',dtype = np.float32,**vrt_options) as dst:
for idd,filename in enumerate(input_files,start=1):
with Rasterio.open(filename,'r') as src:
with WarpedVRT(src,**vrt_options) as vrt:
# At this point 'vrt' is a full dataset with dimensions,and spatial extent matching 'vrt_options'.
# making sure that the data only contains one band,therefore an 2Darray (xsize,ysize)
if src.count == 1:
if windowing:
# Alternatie for processing the dataset in chunks. Likely not very efficient.
for _,window in vrt.block_windows():
data = vrt.read(window=window).astype(np.float32).squeeze(0)
dst.write_band(idd,data,window=window)
# Dump the aligned data into a new file. A VRT representing
# this transformation can also be produced by switching
# to the VRT driver.
else:
# Read all data into memory.
data = vrt.read().astype(np.float32).squeeze(0)
dst.write_band(idd,data)
print('{0}'.format(stacked_filename),'\t is complete')
通过应用虚拟驱动程序进行重采样和重投影（请参见上面的函数“ resample_and_save_via_vrt”），我注意到了两个主要问题：
-
首先，生成的Tiff文件后没有元文件，该文件应包含Tiff数据集的投影和CRS信息。因此，我了解到所有地理转换和CRS信息都直接存储在生成的Tiff文件中。
-
第二，该脚本需要Rasterio的rio_shutil函数。这是为什么？为什么不直接使用标准（即“ tiff”）驱动程序来编写结果数据集？它对结果数据集有什么区别？看起来，它不会创建tiff的元数据。
以下是初始列表的第二项：
用于创建数据集的虚拟驱动程序和标准（即：Tiff）驱动程序之间有什么区别。
如果使用以下代码，则结果数据集将遵循Tiff文件的标准约定，其中将有一个结果tiff文件，后跟一个元数据文件。该元数据将包含该生成的Tiff文件的地理转换和CRS信息。
以下是代码段：
filepath = r'C:\original_file.tif'
dirpath = r'C:\stacked'
dst_crs = 'epsg:5880'
destinations,dst_transforms = {},{}
with Rasterio.open(filepath) as src:
old_transform = src.transform
transform,width,height = calculate_default_transform(
src.crs,dst_crs,src.width,src.height,*src.bounds,resolution=(40,40))
kwargs = src.Meta.copy()
kwargs.update({
'crs': dst_crs,'transform': transform,'width': width,'height': height
})
name,ending = os.path.basename(filepath).split('.')
new_name = name + '_reprojected_to_epsg{0}_using_conventional_driver'.format(dst_crs.split(':')[1]) + '.tif'
print(new_name)
with Rasterio.open(os.path.join(dirpath,new_name),**kwargs) as dst:
for i in range(1,src.count + 1):
reproject(
source=Rasterio.band(src,i),destination=Rasterio.band(dst,src_transform=src.transform,src_crs=src.crs,dst_transform=transform,dst_crs=dst_crs,resampling=resampling.nearest)
对于创建的结果数据集，使用虚拟和物理光栅驱动程序进行重新投影和重采样有什么区别？
此致
解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）