Python3.8: 最后一个输出文件没有正确存储在磁盘上
Python3.8: The last output file is not stored properly on disk
我有一个大约 300m 分辨率的 tif 全球数据集。我想将它升级到 9km 分辨率(在下面你可以看到我的代码)。由于高分辨率数据和大量计算时间,我决定分段进行放大。所以我把整个全球数据分成10块,做upscaling,把每块单独存到一个tif文件里。现在我的问题出现了:最后一段全局数据没有完全保存在磁盘上。每张地图应该是 2M,但 piece#10 是 1.7M。奇怪的是,我的脚本运行两次后,第10块就完成了,从1.7M变成了2M。但是目前的piece10又是不完整的
import numpy as np
from osgeo import gdal
from osgeo import osr
from osgeo.gdalconst import *
import pandas as pd
#
#%%
#-----converting--------#
df_new = pd.read_excel("input_attribute_table.xlsx",sheet_name='Global_data')
listvar = ['var1']
number = df_new['data_number'][:]
##The size of global array is 129599 x 51704. The pieces should be square
xoff = np.array([0, 25852.00, 51704.00, 77556.00, 103408.00])
yoff = np.array([0, 25852.00])
xcount = 25852
ycount = 25852
o = 1
for q in range(len(yoff)):
for p in range(len(xoff)):
src = gdal.Open('Global_database.tif')
ds_xform = src.GetGeoTransform()
ds_driver = gdal.GetDriverByName('Gtiff')
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
data =src.GetRasterBand(1).ReadAsArray(xoff[p],yoff[q],xcount,ycount).astype(np.float32)
Var = np.zeros(data.shape, dtype=np.float32)
Variable_load = df_new[listvar[0]][:]
for m in range(len(number)):
Var[data==number[m]] = Variable_load[m]
#-------rescaling-----------#
Var[np.where(np.isnan(Var))]=0
ds_driver = gdal.GetDriverByName('Gtiff')
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
sz = Var.itemsize
h,w = Var.shape
bh, bw = 36, 36
shape = (h/bh, w/bw, bh, bw)
shape2 = (int(shape[0]),int(shape[1]),shape[2],shape[3])
strides = sz*np.array([w*bh,bw,w,1])
blocks = np.lib.stride_tricks.as_strided(Var,shape=shape2,strides=strides)
resized_array=ds_driver.Create(str(listvar[0])+'_resized_to_9km_glob_piece'+str(o)+'.tif',shape2[1],shape2[0],1,gdal.GDT_Float32) resized_array.SetGeoTransform((ds_xform[0],ds_xform[1]*bw,ds_xform[2],ds_xform[3],ds_xform[4],ds_xform[5]*bh))
resized_array.SetProjection(srs.ExportToWkt())
band = resized_array.GetRasterBand(1)
zero_array = np.zeros([shape2[0],shape2[1]], dtype=np.float32)
for z in range(len(blocks)):
for k in range(len(blocks)):
zero_array[z][k] = np.mean(blocks[z][k])
band.WriteArray(zero_array)
band.FlushCache()
band = None
del zero_array
del Var
o=o+1
通常,您应该确保在文件上调用 close
,或者使用 with
语句。但是,gdal
.
似乎都不支持这些
相反,您应该删除对该文件的所有引用。您已经设置了 band = None
,但您还需要设置 src = None
。
这是一个糟糕的非 Pythonic 界面,但这显然是 Python gdal 库的作用。除了本身就是一个奇怪的陷阱之外,它与异常的交互也很差;任何未处理的异常也可能导致文件未保存(或部分保存或损坏)。
不过,对于眼前的问题,添加 src = None
或 del src
应该可以解决问题。
PS(来自评论):另一种选择是将 for
循环的主体移动到一个函数中;这将自动删除所有变量,而您不必将它们全部列出并可能遗漏一个。有异常还是会出问题,但至少正常情况下应该开始工作了...
我有一个大约 300m 分辨率的 tif 全球数据集。我想将它升级到 9km 分辨率(在下面你可以看到我的代码)。由于高分辨率数据和大量计算时间,我决定分段进行放大。所以我把整个全球数据分成10块,做upscaling,把每块单独存到一个tif文件里。现在我的问题出现了:最后一段全局数据没有完全保存在磁盘上。每张地图应该是 2M,但 piece#10 是 1.7M。奇怪的是,我的脚本运行两次后,第10块就完成了,从1.7M变成了2M。但是目前的piece10又是不完整的
import numpy as np
from osgeo import gdal
from osgeo import osr
from osgeo.gdalconst import *
import pandas as pd
#
#%%
#-----converting--------#
df_new = pd.read_excel("input_attribute_table.xlsx",sheet_name='Global_data')
listvar = ['var1']
number = df_new['data_number'][:]
##The size of global array is 129599 x 51704. The pieces should be square
xoff = np.array([0, 25852.00, 51704.00, 77556.00, 103408.00])
yoff = np.array([0, 25852.00])
xcount = 25852
ycount = 25852
o = 1
for q in range(len(yoff)):
for p in range(len(xoff)):
src = gdal.Open('Global_database.tif')
ds_xform = src.GetGeoTransform()
ds_driver = gdal.GetDriverByName('Gtiff')
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
data =src.GetRasterBand(1).ReadAsArray(xoff[p],yoff[q],xcount,ycount).astype(np.float32)
Var = np.zeros(data.shape, dtype=np.float32)
Variable_load = df_new[listvar[0]][:]
for m in range(len(number)):
Var[data==number[m]] = Variable_load[m]
#-------rescaling-----------#
Var[np.where(np.isnan(Var))]=0
ds_driver = gdal.GetDriverByName('Gtiff')
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
sz = Var.itemsize
h,w = Var.shape
bh, bw = 36, 36
shape = (h/bh, w/bw, bh, bw)
shape2 = (int(shape[0]),int(shape[1]),shape[2],shape[3])
strides = sz*np.array([w*bh,bw,w,1])
blocks = np.lib.stride_tricks.as_strided(Var,shape=shape2,strides=strides)
resized_array=ds_driver.Create(str(listvar[0])+'_resized_to_9km_glob_piece'+str(o)+'.tif',shape2[1],shape2[0],1,gdal.GDT_Float32) resized_array.SetGeoTransform((ds_xform[0],ds_xform[1]*bw,ds_xform[2],ds_xform[3],ds_xform[4],ds_xform[5]*bh))
resized_array.SetProjection(srs.ExportToWkt())
band = resized_array.GetRasterBand(1)
zero_array = np.zeros([shape2[0],shape2[1]], dtype=np.float32)
for z in range(len(blocks)):
for k in range(len(blocks)):
zero_array[z][k] = np.mean(blocks[z][k])
band.WriteArray(zero_array)
band.FlushCache()
band = None
del zero_array
del Var
o=o+1
通常,您应该确保在文件上调用 close
,或者使用 with
语句。但是,gdal
.
相反,您应该删除对该文件的所有引用。您已经设置了 band = None
,但您还需要设置 src = None
。
这是一个糟糕的非 Pythonic 界面,但这显然是 Python gdal 库的作用。除了本身就是一个奇怪的陷阱之外,它与异常的交互也很差;任何未处理的异常也可能导致文件未保存(或部分保存或损坏)。
不过,对于眼前的问题,添加 src = None
或 del src
应该可以解决问题。
PS(来自评论):另一种选择是将 for
循环的主体移动到一个函数中;这将自动删除所有变量,而您不必将它们全部列出并可能遗漏一个。有异常还是会出问题,但至少正常情况下应该开始工作了...