【发布时间】：2021-07-21 14:35:25
【问题描述】：
我正在堆叠大量栅格来计算 2 个月的卫星数据的中位数，当数据的分辨率为 10m 时，这可以正常工作。
由于我在 20m 分辨率数据上运行相同的函数（-> 栅格应该是列和行的 1/2），所以我遇到了内存错误。
除了初始卫星数据的波段之外，我没有改变任何东西。
我知道有大量数据，因为这既是很长的时间，也是很大的空间范围，但它仍然适用于较小的分辨率。
我正在 Anaconda 中使用 python3.6 开发一个虚拟机，该机器有 128 GB RAM 和 16 个 VCPus。
错误信息：
<class 'numpy.core._exceptions._ArrayMemoryError'>, ((563, 256, 55296), dtype('int64')) -> always
<class 'MemoryError'>, ((506, 256, 54528), dtype('bool')) -> sometimes
以下是合并代码，其中 file_list 是链接：
import os
from typing import List
from osgeo import gdal
import numpy as np
import glob
def build_vrt(vrt: str, files: List[str], resample_name: str) -> None:
"""builds .vrt file which will hold information needed for overlay
Args:
vrt (:obj:`string`): name of vrt file, which will be created
files (:obj:`list`): list of file names for merging
resample_name (:obj:`string`): name of resampling method
"""
options = gdal.BuildVRTOptions(srcNodata=-9999)
gdal.BuildVRT(destName=vrt, srcDSOrSrcDSTab=files, options=options)
add_pixel_fn(vrt, resample_name)
def add_pixel_fn(filename: str, resample_name: str) -> None:
"""inserts pixel-function into vrt file named 'filename'
Args:
filename (:obj:`string`): name of file, into which the function will be inserted
resample_name (:obj:`string`): name of resampling method
"""
header = """ <VRTRasterBand dataType="uInt16" band="1" subClass="VRTDerivedRasterBand">"""
contents = """
<PixelFunctionType>{0}</PixelFunctionType>
<PixelFunctionLanguage>Python</PixelFunctionLanguage>
<PixelFunctionCode><![CDATA[{1}]]>
</PixelFunctionCode>"""
lines = open(filename, 'r').readlines()
lines[3] = header # FIX ME: 3 is a hand constant
lines.insert(4, contents.format(resample_name,
get_resample(resample_name)))
open(filename, 'w').write("".join(lines))
def get_resample(name: str) -> str:
"""retrieves code for resampling method
Args:
name (:obj:`string`): name of resampling method
Returns:
method :obj:`string`: code of resample method
"""
methods = {
"median":
"""
import numpy as np
def median(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,raster_ysize, buf_radius, gt, **kwargs):
div = np.zeros((len(in_ar),in_ar[0].shape[0],in_ar[0].shape[1]), dtype=np.float16)
for i in range(len(in_ar)):
div[i,:,:] = np.where(in_ar[i] != 0,in_ar[i],np.nan)
y = np.nanmedian(div, axis=0)
np.clip(y,y.min(),y.max(), out = out_ar)
"""}
if name not in methods:
raise ValueError(
"ERROR: Unrecognized resampling method (see documentation): '{}'.".
format(name))
return methods[name]
def merge(files: List[str], output_file: str, resample: str = "average") -> None:
"""merges list of files using specific resample method for overlapping parts
Args:
files (:obj:`list[string]`): list of files to merge
output_file (:obj:`string`): name of output file
resample (:obj:`string`): name of resampling method
"""
#des=r"E:\naser\code_de\output\test\_vrt.vrt"
des=os.getcwd() + "/_vrt.vrt"
print("1")
build_vrt(des, files, resample)
print("2")
gdal.SetConfigOption('GDAL_VRT_ENABLE_PYTHON', 'YES')
print("3")
translateoptions = gdal.TranslateOptions(gdal.ParseCommandLine("-of Gtiff -ot UINT16 -co TILED=YES -co COMPRESS=LZW BIGTIFF=YES NUM_THREADS=ALL_CPUS -a_nodata 0"))
gdal.SetConfigOption("GDAL_CACHEMAX","512")
gdal.Translate(destName=output_file, srcDS=des, options=translateoptions)
print("4")
gdal.SetConfigOption('GDAL_VRT_ENABLE_PYTHON', None)
print("5")
if os.path.isfile(des):
os.remove(des)
def mergeAll(file_list,Outname,resample):
merge(file_list,Outname,resample)
对于为什么会以这种方式发生，是否有合理的解释？或者我能做什么？
【问题讨论】：
-
一个 int64 形状数组 (563, 256, 55296) 已经需要 60 GB 内存。如果你有另一个这样的数组，那么你的内存几乎已经耗尽，所以仔细检查你有多少个大数组。
-
这是一个好点，我将不得不再次检查。谢谢。