全部版块 我的主页
论坛 提问 悬赏 求职 新闻 读书 功能一区 经管文库(原现金交易版)
1325 0
2021-12-21
全国+各省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

数据如下图
1.jpg
5.jpg 3.jpg 2.jpg
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

相关推荐
栏目导航
热门文章
推荐文章

分享

扫码加好友,拉您进群
各岗位、行业、专业交流群