今天试了好多方法去除TIFF黑边,用Arcgis怎么都没法解决,也不知道是啥原因~
大概是因为这个TIFF太坚强了~捂脸~~
最后用Python GDAL终于将NoData设为了0,最后用Arcgis终于实现了Mosaic,在此做个笔记。
实现代码
以下为NoData设为0并存为tiff的代码:
from osgeo import gdal
import numpy as np
import os
'''
路径
'''
#path = "I:\\Folder\\"
##name_l = "landsat.tif"
#name_in = "lst5_30_129043_20010302Copy.tif"
#name_out = "test1.tif"
'''
处理
'''
def dealfunc(path, name_in, name_out):
'''
读取tif
'''
name_l = name_in
#np.set_printoptions(threshold=np.inf)#使print大量数据不用符号...代替而显示所有
dataset = gdal.Open(path + name_l)
print(dataset.GetDescription())#数据描述
print(dataset.RasterCount)#波段数
im_bands = dataset.RasterCount #波段数
print(dataset.GetGeoTransform())#获取仿射矩阵信息
im_geotrans = dataset.GetGeoTransform()#获取仿射矩阵信息
print(dataset.GetProjection())#获取投影信息
im_proj = dataset.GetProjection()#获取投影信息
cols=dataset.RasterXSize#图像宽度
im_width = dataset.RasterXSize #栅格矩阵的列数
rows=(dataset.RasterYSize)#图像长度
im_height = dataset.RasterYSize #栅格矩阵的行数
xoffset=cols/2
yoffset=rows/2
band = dataset.GetRasterBand(1)#取第1波段
#r = band.ReadAsArray(xoffset,yoffset,10,10)#从数据的中心位置位置开始,取10行10列数据
#im_data = dataset.ReadAsArray(0,0,cols,rows)#获取数据
g = band.ReadAsArray(0,0,cols,rows)#获取数据
#im_blueBand = im_data[0,0:cols,0:rows]#获取蓝波段
'''
显示
'''
import cv2
import matplotlib.pyplot as plt
#img2=cv2.merge([g,g,g])
plt.imshow(g)
plt.xticks([]),plt.yticks([]) # 不显示坐标轴
plt.show()
'''
替换nan
'''
import pandas as pd
df = pd.DataFrame(g)
df = df.fillna(0)
a = df.values
'''
输出替换后的tif
'''
# path = path + "bbb.tif"
path = path + name_out
im_data = a
if 'int8' in im_data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
elif len(im_data.shape) == 2:
im_data = np.array([im_data])
else:
im_bands, (im_height, im_width) = 1,im_data.shape
#创建文件
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(path, im_width, im_height, im_bands, datatype)
if(dataset!= None):
dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数
dataset.SetProjection(im_proj) #写入投影
for i in range(im_bands):
dataset.GetRasterBand(i+1).WriteArray(im_data[i])
del dataset
print("\n Finished!!! \n")
'''
主函数
'''
path = "I:\\Folder\\"
#单景图
#name_in = "lst5_30_129043_20010302Copy.tif"
#name_out = "test1.tif"
#dealfunc(path, name_in, name_out)
#多景图
#得到文件夹下的所有文件名称
names = []
files= os.listdir(path)
for i in files:
if os.path.splitext(i)[1] == '.tif':
names.append(i)
for j in names:
name_in = j
name_out = "z_" + j
dealfunc(path, name_in, name_out)
print("\n"+ j +"\n")
print("\n Finished!!! \n")
参考资料
https://blog.csdn.net/t46414704152abc/article/details/77482747
注意问题
经过以上转换,用Arcgis Mosaic的时候,记得修改Mosaic Operator,否则镶嵌还是没法达到预期。