ArcGIS栅格裁剪(多对多)

我有一百多个栅格和一个图幅框shapefle,shapefle有很多个矢量面,每一个矢量面对应一个栅格。我想要做的是,用这个图幅框shapefle对一百多个栅格按范围裁剪。前提是不能把栅格镶嵌在一起,因为镶嵌之后各图层之间有压盖。如图所示。

img

基于new BIng精确模式的回答:
您可以使用Python的gdal库来实现栅格裁剪。首先,您需要使用ogr模块读取矢量面的几何信息,然后使用gdal.Warp函数对每个栅格进行裁剪。这样,您就可以根据图幅框shapefile的范围对每个栅格进行裁剪,而不需要将它们镶嵌在一起。

下面是一个简单的示例代码:

from osgeo import gdal, ogr

# 打开矢量文件
driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open('your_shapefile.shp', 0)
layer = dataSource.GetLayer()

# 遍历每个矢量面
for feature in layer:
    geom = feature.GetGeometryRef()
    wkt = geom.ExportToWkt()

    # 设置裁剪选项
    warp_options = gdal.WarpOptions(cutlineDSName='your_shapefile.shp',
                                    cutlineWhere="id='{}'".format(feature.GetField('id')),
                                    cropToCutline=True,
                                    dstNodata=0)

    # 对每个栅格进行裁剪
    for raster in rasters:
        gdal.Warp('output_{}'.format(raster), raster, options=warp_options)

它首先使用ogr模块打开矢量文件,并遍历每个矢量面,获取其几何信息。然后,它设置裁剪选项,包括指定裁剪范围、裁剪后的空值等。最后,它遍历每个栅格,并使用gdal.Warp函数对其进行裁剪。

该回答通过自己思路及引用到GPTᴼᴾᴱᴺᴬᴵ搜索,得到内容具体如下:
可以使用ArcPy库中的ExtractByMask函数来实现栅格数据的裁剪操作。该函数可以将一个栅格数据集裁剪为一个指定区域内的子集,可以实现多对多的裁剪操作。

以下是使用Python编写的示例代码:

import arcpy
import os

# 设置工作空间
arcpy.env.workspace = r"C:\data\input"

# 获取所有栅格数据的文件名
raster_files = arcpy.ListRasters()

# 循环遍历每个矢量面,对应裁剪一个栅格数据
for shapefile in arcpy.ListFeatureClasses():
    # 获取矢量面的名称
    shapefile_name = os.path.splitext(shapefile)[0]
    # 循环遍历每个栅格数据
    for raster_file in raster_files:
        # 获取栅格数据的名称
        raster_file_name = os.path.splitext(raster_file)[0]
        # 如果矢量面名称和栅格数据名称相同,则进行裁剪操作
        if shapefile_name == raster_file_name:
            # 构建裁剪区域
            mask = shapefile
            # 构建裁剪后输出的栅格数据
            out_raster = os.path.join(r"C:\data\output", raster_file)
            # 执行裁剪操作
            arcpy.gp.ExtractByMask_sa(raster_file, mask, out_raster)

在上述代码中,首先设置工作空间为栅格数据所在的文件夹。然后使用ListRasters函数获取所有栅格数据的文件名,使用ListFeatureClasses函数获取所有矢量面的文件名。接着,循环遍历每个矢量面和每个栅格数据,如果矢量面名称和栅格数据名称相同,则使用ExtractByMask函数进行裁剪操作,并将裁剪后的栅格数据输出到指定文件夹中。

需要注意的是,裁剪操作会生成新的栅格数据,因此需要指定输出路径和文件名。另外,要确保裁剪后的栅格数据不会相互压盖,可以在输出时添加后缀或者使用不同的文件夹来存储裁剪后的数据。


如果以上回答对您有所帮助,点击一下采纳该答案~谢谢

以下答案由GPT-3.5大模型与博主波罗歌共同编写:
ArcGIS栅格裁剪(多对多)的步骤如下:

  1. 打开ArcMap软件,加载栅格数据和图幅框shapefile。

  2. 打开ArcToolbox工具箱,选择Spatial Analyst Tools工具集下的Extraction工具集,找到Extract by Mask工具,双击打开。

  3. 在Extract by Mask对话框中,选择栅格数据图层作为Input Raster,选择图幅框shapefile图层作为Input mask features。

  4. 在Output raster datasets中选择输出文件夹和文件名(可以用原文件名加后缀“_clip”)。

  5. 点击Environments按钮,在Processing Extent栏,选择“Same as Input mask features”选项,保证裁剪后的栅格数据范围与图幅框相同。

  6. 点击OK开始执行裁剪操作。此时,程序会自动遍历shapefile中的每一个矢量面,分别对应裁剪一个栅格文件。

以下是对上述步骤的代码实现:

import arcpy
from arcpy.sa import *

#设置工作空间
arcpy.env.workspace = r"C:\data"
arcpy.CheckOutExtension("Spatial")

#定义栅格数据文件夹和图幅框shapefile文件
raster_folder = "rasters"
shapefile = "masks.shp"

#遍历shapefile中每一个矢量面,对应裁剪一个栅格文件
with arcpy.da.SearchCursor(shapefile, ['OID@', 'SHAPE@']) as cursor:
    for row in cursor:
        #设置输出栅格文件名
        output_raster = raster_folder + "/raster_" + str(row[0]) + "_clip.tif"
        #裁剪栅格
        out_extract = ExtractByMask(row[1], "input_raster.tif")
        #保存裁剪后的栅格
        out_extract.save(output_raster)

arcpy.CheckInExtension("Spatial")
print("裁剪完成!")

其中,SearchCursor函数用于遍历图幅框shapefile中的每一个矢量面,ExtractByMask函数用于裁剪栅格文件,save函数用于保存裁剪后的栅格文件。
如果我的回答解决了您的问题,请采纳!

以下内容部分参考ChatGPT模型:


Python解决思路:

  1. 读取图幅框shapefile和栅格数据,使用geopandas和rasterio库实现。
  2. 遍历图幅框shapefile中的所有矢量面,对每个矢量面进行裁剪操作。
  3. 对每个矢量面,使用geopandas的geometry属性获取其几何形状。
  4. 对每个栅格数据,使用rasterio库的mask函数进行裁剪操作,裁剪范围为当前矢量面的几何形状。
  5. 将裁剪后的栅格数据保存到指定路径。

Python示例代码:

import geopandas as gpd
import rasterio
from rasterio.mask import mask

# 读取图幅框shapefile和栅格数据
grid_files = ["grid1.tif", "grid2.tif", ...]
grid_shapes = gpd.read_file("grid_shapes.shp")

# 遍历图幅框shapefile中的所有矢量面,对每个矢量面进行裁剪操作
for i, shape in grid_shapes.iterrows():
    # 获取当前矢量面的几何形状
    geom = shape.geometry
    
    # 遍历所有栅格数据,对当前矢量面进行裁剪操作
    for grid_file in grid_files:
        # 读取栅格数据
        with rasterio.open(grid_file) as src:
            # 裁剪栅格数据
            out_image, out_transform = mask(src, [geom], crop=True)
            out_meta = src.meta.copy()

        # 保存裁剪后的栅格数据到指定路径
        with rasterio.open("clipped_" + grid_file, "w", **out_meta) as dest:
            dest.write(out_image)

C++解决思路:

  1. 读取图幅框shapefile和栅格数据,使用GDAL库实现。
  2. 遍历图幅框shapefile中的所有矢量面,对每个矢量面进行裁剪操作。
  3. 对每个矢量面,使用GDAL库的OGRGeometry类获取其几何形状。
  4. 对每个栅格数据,使用GDAL库的GDALRasterBand类进行裁剪操作,裁剪范围为当前矢量面的几何形状。
  5. 将裁剪后的栅格数据保存到指定路径。

C++示例代码:

#include "gdal_priv.h"
#include "ogr_geometry.h"

int main()
{
    // 读取图幅框shapefile和栅格数据
    std::vector<std::string> grid_files = {"grid1.tif", "grid2.tif", ...};
    GDALDataset *grid_shapes = (GDALDataset*)GDALOpenEx("grid_shapes.shp", GDAL_OF_VECTOR, NULL, NULL, NULL);

    // 遍历图幅框shapefile中的所有矢量面,对每个矢量面进行裁剪操作
    for (int i = 0; i < grid_shapes->GetLayerCount(); i++)
    {
        OGRLayer *layer = grid_shapes->GetLayer(i);

        // 获取当前矢量面的几何形状
        OGRFeature *feature = layer->GetNextFeature();
        OGRGeometry *geom = feature->GetGeometryRef();

        // 遍历所有栅格数据,对当前矢量面进行裁剪操作
        for (std::string grid_file : grid_files)
        {
            // 读取栅格数据
            GDALDataset *grid_data = (GDALDataset*)GDALOpenEx(grid_file.c_str(), GDAL_OF_RASTER, NULL, NULL, NULL);
            GDALRasterBand *grid_band = grid_data->GetRasterBand(1);

            // 裁剪栅格数据
            GDALWarpOptions *options = GDALWarpOptionsNew(NULL, NULL);
            options->papszWarpOptions = CSLSetNameValue(options->papszWarpOptions, "CUTLINE", OGR_G_ExportToWkt(geom));
            options->papszWarpOptions = CSLSetNameValue(options->papszWarpOptions, "CROP_TO_CUTLINE", "YES");
            GDALDataset *clipped_data = GDALWarp("clipped_" + grid_file, NULL, 1, &grid_data, options, NULL);
            GDALClose(clipped_data);
            GDALWarpOptionsFree(options);

            // 释放资源
            GDALClose(grid_data);
        }

        // 释放资源
        OGRFeature::DestroyFeature(feature);
    }

    // 释放资源
    GDALClose(grid_shapes);

    return 0;
}

Java解决思路:

  1. 读取图幅框shapefile和栅格数据,使用GeoTools库实现。
  2. 遍历图幅框shapefile中的所有矢量面,对每个矢量面进行裁剪操作。
  3. 对每个矢量面,使用GeoTools的JTS类获取其几何形状。
  4. 对每个栅格数据,使用GeoTools的GridCoverage2D类进行裁剪操作,裁剪范围为当前矢量面的几何形状。
  5. 将裁剪后的栅格数据保存到指定路径。

Java示例代码:

import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.processing.CoverageProcessor;
import org.geotools.geometry.jts.JTS;
import org.geotools.geometry.jts.JTSFactoryFinder;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.gce.arcgrid.ArcGridReader;
import org.geotools.gce.geotiff.GeoTiffReader;
import org.geotools.geometry.Envelope2D;
import org.geotools.referencing.CRS;
import org.opengis.coverage.grid.GridCoverage;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.simple.SimpleFeatureType;
import org.opengis.geometry.Geometry;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.util.ProgressListener;

import java.io.File;
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;

public class Main {
    public static void main(String[] args) throws IOException {
        // 读取图幅框shapefile和栅格数据
        List<String> grid_files = List.of("grid1.tif", "grid2.tif", ...);
        File grid_shapes = new File("grid_shapes.shp");
        CoordinateReferenceSystem crs = CRS.decode("EPSG:4326");

        // 遍历图幅框shapefile中的所有矢量面,对每个矢量面进行裁剪操作
        SimpleFeatureType schema = DataUtilities.createType("grid_shapes", "geometry:Polygon:srid=4326");
        FeatureReader<SimpleFeatureType, SimpleFeature> reader = new ShapefileReader(grid_shapes, schema, true);
        while (reader.hasNext()) {
            SimpleFeature feature = reader.next();

            // 获取当前矢量面的几何形状
            Geometry geom = (Geometry) feature.getDefaultGeometry();
            ReferencedEnvelope envelope = new ReferencedEnvelope(geom.getEnvelopeInternal(), crs);

            // 遍历所有栅格数据,对当前矢量面进行裁剪操作
            for (String grid_file : grid_files) {
                // 读取栅格数据
                GridCoverage2D grid = null;
                if (grid_file.endsWith(".tif")) {
                    grid = new GeoTiffReader(new File(grid_file)).read(null);
                } else if (grid_file.endsWith(".asc")) {
                    grid = new ArcGridReader(new File(grid_file)).read(null);
                }

                // 裁剪栅格数据
                CoverageProcessor processor = new CoverageProcessor();
                GridCoverage coverage = processor.crop(grid, envelope, null);

                // 保存裁剪后的栅格数据到指定路径
                File clipped_file = new File("clipped_" + grid_file);
                GridCoverageWriter writer = new GeoTiffWriter(clipped_file);
                writer.write(coverage, null);

                // 释放资源
                writer.dispose();
            }
        }

        // 释放资源
        reader.close();
    }
}

如果我的建议对您有帮助、请点击采纳、祝您生活愉快

ArcGIS:矢量、栅格文件裁剪(批量处理)
可以借鉴下
https://blog.csdn.net/iamWinnie_/article/details/127095524