我有一百多个栅格和一个图幅框shapefle,shapefle有很多个矢量面,每一个矢量面对应一个栅格。我想要做的是,用这个图幅框shapefle对一百多个栅格按范围裁剪。前提是不能把栅格镶嵌在一起,因为镶嵌之后各图层之间有压盖。如图所示。
基于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栅格裁剪(多对多)的步骤如下:
打开ArcMap软件,加载栅格数据和图幅框shapefile。
打开ArcToolbox工具箱,选择Spatial Analyst Tools工具集下的Extraction工具集,找到Extract by Mask工具,双击打开。
在Extract by Mask对话框中,选择栅格数据图层作为Input Raster,选择图幅框shapefile图层作为Input mask features。
在Output raster datasets中选择输出文件夹和文件名(可以用原文件名加后缀“_clip”)。
点击Environments按钮,在Processing Extent栏,选择“Same as Input mask features”选项,保证裁剪后的栅格数据范围与图幅框相同。
点击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函数用于保存裁剪后的栅格文件。
如果我的回答解决了您的问题,请采纳!
Python解决思路:
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++解决思路:
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解决思路:
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