TIFF裁切
步骤如下：
- 读取xml：读取tiff对应的xml文件，得到tiff文件的经纬度坐标和wkt范围；
- 相交分析：将用户输入的wkt范围与tiff源文件范围进行相交，得到新的wkt范围；
- wkt转shp：将相交分析后的wkt转为shp文件；
- rpc校正：某些tiff中缺少坐标系信息，故需要进行rpc校正，将rpb文件中的信息嵌入到tiff中，得到新的tiff；
- tiff裁切：用shp去裁切rpc校正后的tiff；
- 生成缩略图：tiff转图片，图片格式为jpg或png。
读取xml
def parse_xml(eles):
"""
适用于除GF3以外的tiff以及部分H1D类型的h5
:param eles: 文档元素树
:return: wkt_range: wkt范围
"""
TopLeftLatitude = eles.getElementsByTagName("TopLeftLatitude")[0].childNodes[0].data # 左上 纬度
TopLeftLongitude = eles.getElementsByTagName("TopLeftLongitude")[0].childNodes[0].data # 左上 经度
TopRightLatitude = eles.getElementsByTagName("TopRightLatitude")[0].childNodes[0].data # 右上 纬度
TopRightLongitude = eles.getElementsByTagName("TopRightLongitude")[0].childNodes[0].data # 右上 经度
BottomRightLatitude = eles.getElementsByTagName("BottomRightLatitude")[0].childNodes[0].data # 右下 纬度
BottomRightLongitude = eles.getElementsByTagName("BottomRightLongitude")[0].childNodes[0].data # 右下经度
BottomLeftLatitude = eles.getElementsByTagName("BottomLeftLatitude")[0].childNodes[0].data # 左下 纬度
BottomLeftLongitude = eles.getElementsByTagName("BottomLeftLongitude")[0].childNodes[0].data # 左下 经度
wkt_range = f"POLYGON (({TopLeftLongitude} {TopLeftLatitude},{BottomLeftLongitude} {BottomLeftLatitude}," \
f"{BottomRightLongitude} {BottomRightLatitude},{TopRightLongitude} {TopRightLatitude}," \
f"{TopLeftLongitude} {TopLeftLatitude}))"
min_lon = min(TopLeftLongitude, BottomLeftLongitude, BottomRightLongitude, TopRightLongitude)
max_lon = max(TopLeftLongitude, BottomLeftLongitude, BottomRightLongitude, TopRightLongitude)
min_lat = min(TopLeftLatitude, BottomLeftLatitude, BottomRightLatitude, TopRightLatitude)
max_lat = max(TopLeftLatitude, BottomLeftLatitude, BottomRightLatitude, TopRightLatitude)
lon_range = float(max_lon) - float(min_lon)
lat_range = float(max_lat) - float(min_lat)
return wkt_range, lon_range, lat_range
相交分析
def xmlIntersect(wkt_tiff, wkt_range):
"""
相交分析: 将图像范围与输入的待裁剪范围进行相交, 得到交集部分
:param wkt_tiff: 图像的wkt范围
:param wkt_range: 待裁切的wkt范围
:return: intersection.ExportToWkt(): 经纬度范围, 为wkt格式
"""
poly1 = ogr.CreateGeometryFromWkt(wkt_tiff) # 图像范围wkt
poly2 = ogr.CreateGeometryFromWkt(wkt_range) # 裁剪范围wkt
intersection = poly1.Intersection(poly2)
logger.info(f"图像范围: {poly1}")
logger.info(f"裁剪范围: {poly2}")
logger.info(f"相交分析成功: {intersection.ExportToWkt()}")
return intersection.ExportToWkt()
wkt转shp
def wktToShp(wkt_range_intersect, output_shp_file):
"""
根据wkt字符串生成shp文件, 以便用于裁切
:param wkt_range_intersect: wkt字符串, 包含经纬度, 根据此范围进行裁切
:param output_shp_file: 输出的shp文件
"""
wkt_range_string = wkt.loads(wkt_range_intersect)
logger.info("根据wkt生成矢量")
points = geopandas.GeoSeries(wkt_range_string, crs='EPSG:4326') # 坐标系: WGS 1984
gdal.SetConfigOption("GTIFF_SRS_SOURCE", "EPSG")
points.to_file(output_shp_file, driver="ESRI Shapefile", encoding='utf-8')
logger.info(f"根据wkt生成矢量成功: {output_shp_file}")
rpc校正
def rpcCorrect(input_tiff_file, lon_range, lat_range, output_rpc_file):
"""
原tiff文件没有坐标系及相关信息, 故需要进行rpc校正
读取tiff文件路径下的.rpb文件, 输出新的tiff文件, 带有坐标系信息
若无rpb文件则不进行rpc校正, 继续后续裁切
param input_tiff: 输入的原始影像
param output_tiff: rpc校正后输出的正射影像
"""
logger.info(f"开始进行rpc校正: {input_tiff_file}")
dataset = gdal.Open(input_tiff_file, gdal.GA_Update) # 读入影像
rpc = dataset.GetMetadata("RPC") # 读入影像，rpc
# 设置分辨率
x_size, y_size = dataset.RasterXSize, dataset.RasterYSize
x_res = lon_range / x_size
y_res = lat_range / y_size
if rpc != {}:
gdal.Warp(output_rpc_file, dataset, dstSRS='EPSG:4326',
xRes=x_res, # x方向正射影像分辨率
yRes=y_res, # y方向正射影像分辨率
resampleAlg=gdal.GRIORA_Bilinear, # 插值方式，默认为最邻近，我们一般使用双线性
rpc=True, # 使用RPC模型进行校正
)
dataset = None # 关闭文件
logger.info(f"rpc校正成功: {output_rpc_file}")
else:
logger.info("读取rpc为空，跳过rpc校正")
tiff裁切
def cutTiff(tiff_file, shp_file, cut_tiff_file):
"""
核心函数: 实现tiff裁切功能
:param tiff_file: 经过rpc校正后的tiff文件路径
:param shp_file: 根据wkt生成的shp文件路径
:param cut_tiff_file: 裁切后的tiff文件路径
"""
logger.info("开始裁切")
input_raster = gdal.Open(tiff_file)
if input_raster is None:
logger.warning(f"tiff文件无法加载: {tiff_file}")
gdal.Warp(cut_tiff_file,
input_raster,
format='GTiff',
# dstSRS='EPSG:4326',
cutlineDSName=shp_file,
cropToCutline=False,
# 保证裁剪后影像大小跟矢量文件的图框大小一致（设置为False时，结果图像大小会跟待裁剪影像大小一样，则会出现大量的空值区域）
copyMetadata=True,
dstNodata=0,
creationOptions=["COMPRESS=LZW", "TILED=True", "BIGTIFF=YES"],
)
input_raster = None # 关闭文件
logger.info(f"裁切成功: {cut_tiff_file}")
生成缩略图
def tiffToPng(tiff_file, png_file, scale=0.3):
"""
将tiff文件转为png图片
:param tiff_file: tiff文件路径
:param png_file: 生成的图片路径
:param scale: 图像缩放比例
"""
logger.info("开始生成缩略图")
tiff_name = os.path.split(tiff_file)[1]
img = io.imread(tiff_file) # 读取文件名
logger.debug(f"图片尺寸: {img.shape}; 数据类型: {img.dtype}") # 显示图片大小和类型
# GF3类型部分文件变量类型 float16 溢出，需要调整类型为 float32
if tiff_name.startswith("GF3"):
my_type = np.float32
else:
my_type = np.float16
logger.debug(f"ndarray的数据类型: {my_type}")
if len(img.shape) == 3:
img = np.divide(img[:, :, 0:3], img.max(), dtype=my_type) # 使其所有值不大于一
logger.debug(f"处理后的图片尺寸: {img.shape}")
# 几乎全黑，解决色差
if tiff_name.startswith("GF3_KSC_FSI") or tiff_name.startswith("GF3_KRN_SS"):
img = img * 2550
elif tiff_name.startswith("GF3_SYC_QPSI"):
img = img * 30000
else:
img = img * 255
img = img.astype(np.uint8) # 强制转换成8位整型
# 得到图片长宽高
height, width, _ = img.shape
# 图片像素点太多，网页无法正常显示，需要减小像素点数量
img_resize = cv2.resize(img, dsize=(round(width * scale), round(height * scale)))
b = img_resize[:, :, 0] # 读取蓝通道
g = img_resize[:, :, 1] # 读取绿通道
r = img_resize[:, :, 2] # 读取红通道
rgb_img = cv2.merge([b, g, r]) # 通道拼接
# 三通道转四通道，设置黑色为透明
gray_image = cv2.cvtColor(rgb_img, cv2.COLOR_BGR2GRAY)
_, alpha = cv2.threshold(gray_image, 0, 255, cv2.THRESH_BINARY)
rgba_image = cv2.merge([b, g, r, alpha], 4)
# 识别轮廓，裁剪图像
contours, _ = cv2.findContours(alpha, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
x, y, w, h = cv2.boundingRect(contours[0])
# 绘制轮廓
cv2.drawContours(rgba_image, contours, 0, (0, 0, 255), 5)
# 提取图像范围
cropped_image = rgba_image[y:y + h, x:x + w]
# 保存图像
cv2.imwrite(png_file, cropped_image, [int(cv2.IMWRITE_PNG_COMPRESSION), 9])
elif len(img.shape) == 2:
img = np.divide(img, img.max(), dtype=my_type) # 使其所有值不大于一
logger.debug(f"处理后的图片尺寸: {img.shape}")
# 部分GF3类型存在问题: 几乎全黑，解决色差
if tiff_name.startswith("GF3_KSC_FSI") or tiff_name.startswith("GF3_KRN_SS"):
img = img * 2550
elif tiff_name.startswith("GF3_SYC_QPSI"):
img = img * 30000
else:
img = img * 255
img = img.astype(np.uint8) # 强制转换成8位整型
# 得到图片的高度、宽度
height, width = img.shape
img_resize = cv2.resize(img, dsize=(round(width * scale), round(height * scale)))
rgb_image = cv2.cvtColor(img_resize, cv2.COLOR_GRAY2RGB)
b = rgb_image[:, :, 0] # 读取蓝通道
g = rgb_image[:, :, 1] # 读取绿通道
r = rgb_image[:, :, 2] # 读取红通道
_, alpha = cv2.threshold(img_resize, 0, 255, cv2.THRESH_BINARY)
rgba_image = cv2.merge([b, g, r, alpha], 4)
# 识别轮廓，裁剪图像
contours, _ = cv2.findContours(alpha, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
x, y, w, h = cv2.boundingRect(contours[0])
# 绘制轮廓
cv2.drawContours(rgba_image, contours, 0, (0, 0, 255), 5)
# 提取图像范围
cropped_image = rgba_image[y:y + h, x:x + w]
# 保存图像
cv2.imwrite(png_file, cropped_image, [int(cv2.IMWRITE_PNG_COMPRESSION), 9])
try:
del img, b, g, r, rgba_image # 释放内存
gc.collect()
logger.debug(f"释放内存成功，删除变量: img, b, g, r, rgba_image")
except Exception as e:
logger.error(e)
logger.info(f"生成缩略图成功: {png_file}")
参考文章
- 一般使用 - python GDAL/OGR 中文手册 (luolingchun.github.io)
- 2. GDAL python教程（1）——用OGR读写矢量数据 — 开放地理空间实验室: 犹他州立大学——开源GIS类库GDAL资料：Python GDAL课程笔记 (osgeo.cn)
- 30分钟学会shapely空间几何分析_shapely空间计算近邻分析-CSDN博客
- 用Python剪裁栅格数据，以及如何用Python读取shapefile - 简书 (jianshu.com)
- python批量遥感影像裁剪的三种方法（gdal.Warp,arcpy）-CSDN博客
- python-gdal-test: 基于Python-GDAL的实践总结 (gitee.com)
- Polygon—ArcGIS Pro | 文档
- GTiff—GeoTIFF文件格式 — GDAL 文档 (osgeo.cn)
- Python与开源GIS (osgeo.cn)
浙公网安备 33010602011771号