全国+各省ESA10 米土地利用镶嵌数据分享裁剪代码如下:
#裁剪函数
# -*- coding: utf-8 -*-
from geopandas import *
import rasterio as rio
import rasterio.mask
import rasterio
import os
from tqdm import tqdm
import _thread
#全局变量
# 读入矢量
shpdatafile = r'D:\资源\ESA\矢量\各省矢量.shp'
shpdata=GeoDataFrame.from_file(shpdatafile)
#创建各地文件夹函数
def mkdir(path):
# 引入模块
# 去除首位空格
path = path.strip()
# 去除尾部 \ 符号
path = path.rstrip("\\")
# 判断路径是否存在
# 存在 True
# 不存在 False
isExists = os.path.exists(path)
# 判断结果
if not isExists:
# 如果不存在则创建目录
# 创建目录操作函数
os.makedirs(path)
print(path + ' 创建成功')
return True
else:
# 如果目录存在则不创建,并提示目录已存在
print(path + ' 目录已存在')
return False
#ESA原始文件夹
files_path =r"D:\资源\ESA\CHINA"
pathDir= os.listdir(files_path)
,以后数据分享会设置限制条件✌
for i in tqdm(range(len(pathDir))):
# 读入栅格文件
rasterfile = files_path+"\\"+pathDir
rasterdata = rio.open(rasterfile)
#获取栅格信息
profile = rasterdata.profile
#标识符
biaoshi = pathDir
# 投影变换,使矢量数据与栅格数据投影参数一致
shpdata = shpdata.to_crs(rasterdata.crs)
# 按照所有矢量进行循环裁剪
for j in range(0, len(shpdata)):
try:
# 获取矢量数据的features
geo = shpdata.geometry[j]
#获取该要素的属性信息
data_shp_name_sheng=shpdata.NAME[j]
#文件保存位置的文件夹 各省
data_filepath=str(data_shp_name_sheng)
feature = [geo.__geo_interface__]
# 通过feature裁剪栅格影像
out_image, out_transform = rio.mask.mask(rasterdata, feature, all_touched=True, crop=True, nodata=0)
profile.update(
height=out_image.shape[1],
width=out_image.shape[2],
shape=(out_image.shape[1],out_image.shape[2]),
nodata=0,
bounds=[],
transform=out_transform,
)
# 定义要创建的目录
mkpath = "D:\资源\ESA\PRO_CLIP\\" + data_shp_name_sheng+ "\\"
mkdir(mkpath)
#文件名字
name=mkpath+"\\"+data_shp_name_sheng+"_"+biaoshi
with rasterio.open(name, mode='w', **profile) as dst:
dst.write(out_image)
except:
pass
镶嵌代码如下:
# 镶嵌一个文件夹中的所有文件
import os
from geopandas import *
from osgeo import gdal, gdalconst
from tqdm import tqdm
tifPath = r'D:\资源\ESA\PRO_CLIP\\' # 待融合的图像所在的文件夹
tifPaths_folder = os.listdir(tifPath)
# 输出文件位置
path_save = r"D:\\资源\\ESA\\PRO_CLIP_mosaic\\"
# 循环目录
for path in tqdm(tifPaths_folder):
try:
DEM_SMALL_PATH = os.path.join(tifPath, path)
for (pathname, dirs, files) in os.walk(DEM_SMALL_PATH):
son_Paths_file = files
# 如果影像数量大于一
if len(son_Paths_file) >= 2:
DEM_SMALL_PATH2 = DEM_SMALL_PATH + "\\"
# 循环子目录,进行镶嵌
# 循环同一个文件下的tif文件
inputFiles = []
for path_small in son_Paths_file:
# 每一个栅格的路径
son_Paths_PATH = os.path.join(DEM_SMALL_PATH2, path_small)
# 读取影像
inputrasfile = gdal.Open(son_Paths_PATH, gdal.GA_ReadOnly) # 读取影像
inputProj = inputrasfile.GetProjection() # 获取坐标系
inputFiles.append(inputrasfile) # 推入列表
options = gdal.WarpOptions(srcSRS=inputProj, # 输入坐标系
dstSRS=inputProj, # 输出坐标系
format='GTiff', # 图像格式
resampleAlg=gdalconst.GRIORA_NearestNeighbour, # 重采样算法,这里是双线性内插
dstNodata=0, # 缺省值
creationOptions=['COMPRESS=LZW'],
outputType=gdalconst.GDT_Byte)
# 创建文件夹
Landcover_to_save = path_save + path
mkdir(Landcover_to_save)
# 输出文件名
outputfilePath = Landcover_to_save + "\\" + path + "_土地利用数据" + "_10m分辨率_ESA数据_2020年" + ".tif"
# 写栅格
gdal.Warp(outputfilePath, inputFiles, options=options) # 图像镶嵌
except:
pass
数据如下图