今天试了下arcgis 矢量裁剪栅格，首先在arcmap上面操作，发现有好几种方式都能实现：
方法1，矢量转栅格，把像元值设置1(方法是添加一个字段，将字段值全部赋值1）然后调用toolbox的时候，选择该字段；这样就生成了一个所有cell的值为1的栅格RasterA。
接下来，将RasterA与我们要裁剪的栅格对象RasterB进行栅格计算(栅格计算也用toolbox里面的，不要用菜单里面的，因为我们最后代码里面调用的是toolbox)，将RasterA*RasterB就是最终结果啦。
方法2，和方法1的区别就是，不用添加字段，而是重分类一下，Reclassify重分类，将某个字段的值重分类为1。
方法3，也是最屌的方法，直接用Extract by Polygon工具。
3种方法区别：
方法1和方法2算是同一种，他们比方法3稳定，因为他们只是单纯栅格计算，在大数据量的时候比较靠谱。方法3虽然简单，但面对大数据量处理的时候，容易出错，不太稳定。当然小数据的话，还是木有问题的。
下边直接贴代码啦，简单的封装了一个类
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.Geometry;
using ESRI.ArcGIS.Geoprocessor;
using ESRI.ArcGIS.Geoprocessing;
using ESRI.ArcGIS.esriSystem;
using ESRI.ArcGIS.SpatialAnalystTools;
using ESRI.ArcGIS.SpatialAnalyst;
using ESRI.ArcGIS.DataSourcesRaster;
using ESRI.ArcGIS.Carto;
using ESRI.ArcGIS.DataSourcesGDB;
namespace GPOperator
{
public class ClipHelper : GPTool
{
/// <summary>
/// 栅格裁剪
/// </summary>
/// <param name="inputRaster">要裁剪的栅格(支持对象和栅格文件路径)</param>
/// <param name="maskData">掩膜数据(支持对象和shp文件路径)</param>
/// <returns></returns>
public IGeoDataset ClipRaster2(object inputRaster, object maskData, out string resultInfo)
{
string outRaster = System.IO.Path.Combine(System.IO.Path.GetTempPath(), "ClipRaster_" +DateTime.Now.ToString("yyyy-MM-dd HH-mm-ss")+ Guid.NewGuid().ToString
("D").Substring(0,4) + ".tif");
//参数：对象或路径都支持
ESRI.ArcGIS.SpatialAnalystTools.ExtractByMask mask = new ESRI.ArcGIS.SpatialAnalystTools.ExtractByMask(inputRaster, maskData, outRaster);//@"C:\tmp\新建文件夹
\ret.tif"
bool isok = this.Excute(mask, out resultInfo);
//IWorkspaceFactory wsf = new RasterWorkspaceFactory();
//IRasterWorkspace rasterWS = wsf.OpenFromFile(System.IO.Path.GetDirectoryName(outRaster), 0) as IRasterWorkspace;
//IRasterDataset rds = rasterWS.OpenRasterDataset(System.IO.Path.GetFileName(outRaster));
IRasterDataset rasterDataset = new RasterDataset();
rasterDataset.OpenFromFile(outRaster);
return rasterDataset as IGeoDataset;
}
/// <summary>
/// 栅格裁剪
/// </summary>
/// <param name="inputRaster">要裁剪的栅格(支持对象和栅格文件路径)</param>
/// <param name="maskData">裁剪多边形</param>
/// <returns>裁剪得到的栅格</returns>
public IGeoDataset ClipRaster(object inputRaster, IPolygon polygon, out string resultInfo)
{
var maskData = CreateFeatureClass(polygon);
string outRaster = System.IO.Path.Combine(System.IO.Path.GetTempPath(), "ClipRaster_" + DateTime.Now.ToString("yyyy-MM-dd HH-mm-ss") + Guid.NewGuid().ToString
("D").Substring(0, 4) + ".tif");
ESRI.ArcGIS.SpatialAnalystTools.ExtractByMask mask = new ESRI.ArcGIS.SpatialAnalystTools.ExtractByMask(inputRaster, maskData, outRaster);
if (this.Excute(mask, out resultInfo))
{
return null;
}
IRasterDataset rasterDataset = new RasterDataset();
try
{
rasterDataset.OpenFromFile(outRaster);
}
catch (Exception ex)
{
resultInfo = ex.ToString();
return null;
}
finally
{
try
{
ESRI.ArcGIS.ADF.ComReleaser.ReleaseCOMObject(maskData);
}
catch { }
}
return rasterDataset as IGeoDataset;
}
/// <summary>
/// 栅格裁剪
/// Op方法，直接返回对象
/// </summary>
/// <param name="inputRaster">要裁剪的栅格</param>
/// <param name="maskData">掩膜数据</param>
/// <returns></returns>
public IGeoDataset ClipRaster3(IGeoDataset inputRaster, IGeoDataset maskData,out string resultInfo)
{
resultInfo = "执行成功";
IExtractionOp2 op = new RasterExtractionOpClass();
try
{
return op.Raster(inputRaster, maskData);
}
catch (Exception ex)
{
resultInfo = ex.ToString();
return null;
}
}
public IFeatureClass CreateFeatureClass(IPolygon polygon)
{
IFeatureClassDescription fcDescription = new FeatureClassDescriptionClass();
IObjectClassDescription ocDescription = (IObjectClassDescription)fcDescription;
IFeatureClass fc= createFeatureClassInmemeory("temp", "temp", polygon.SpatialReference, esriGeometryType.esriGeometryPolygon, ocDescription.RequiredFields);
if (null == fc)
{
return null;
}
IFeature f= fc.CreateFeature();
f.Shape = polygon;
f.Store();
return fc;
}
//IFeatureCache
/// <summary>
/// 在内存中创建临时要素类
/// </summary>
/// <param name="DataSetName">数据集名称</param>
/// <param name="AliaseName">别名</param>
/// <param name="SpatialRef">空间参考</param>
/// <param name="GeometryType">几何类型</param>
/// <param name="PropertyFields">属性字段集合</param>
/// <returns>IfeatureLayer</returns>
private IFeatureClass createFeatureClassInmemeory(string DataSetName, string AliaseName, ISpatialReference SpatialRef, esriGeometryType GeometryType, IFields
PropertyFields)
{
IWorkspaceFactory workspaceFactory = new InMemoryWorkspaceFactoryClass();
ESRI.ArcGIS.Geodatabase.IWorkspaceName workspaceName = workspaceFactory.Create("", "MyWorkspace", null, 0);
ESRI.ArcGIS.esriSystem.IName name = (IName)workspaceName;
ESRI.ArcGIS.Geodatabase.IWorkspace inmemWor = (IWorkspace)name.Open();
IField oField = new FieldClass();
IFields oFields = new FieldsClass();
IFieldsEdit oFieldsEdit = null;
IFieldEdit oFieldEdit = null;
IFeatureClass oFeatureClass = null;
try
{
oFieldsEdit = oFields as IFieldsEdit;
oFieldEdit = oField as IFieldEdit;
for (int i = 0; i < PropertyFields.FieldCount; i++)
{
IField field = PropertyFields.get_Field(i);
if (field.Type != esriFieldType.esriFieldTypeGeometry)
{
oFieldsEdit.AddField(field);
}
}
IGeometryDef geometryDef = new GeometryDefClass();
IGeometryDefEdit geometryDefEdit = (IGeometryDefEdit)geometryDef;
geometryDefEdit.AvgNumPoints_2 = 5;
geometryDefEdit.GeometryType_2 = GeometryType;
geometryDefEdit.GridCount_2 = 1;
geometryDefEdit.HasM_2 = false;
geometryDefEdit.HasZ_2 = false;
geometryDefEdit.SpatialReference_2 = SpatialRef;
oFieldEdit.Name_2 = "SHAPE";
oFieldEdit.Type_2 = esriFieldType.esriFieldTypeGeometry;
oFieldEdit.GeometryDef_2 = geometryDef;
oFieldEdit.IsNullable_2 = true;
oFieldEdit.Required_2 = true;
oFieldsEdit.AddField(oField);
oFeatureClass = (inmemWor as IFeatureWorkspace).CreateFeatureClass(DataSetName, oFields, null, null, esriFeatureType.esriFTSimple, "SHAPE", "");
(oFeatureClass as IDataset).BrowseName = DataSetName;
return oFeatureClass;
}
catch(Exception ex)
{
return null;
}
}
}
/// <summary>
/// GP工具基类，实现执行工具公共方法和释放资源方法的定义。
/// </summary>
public class GPTool
{
private Geoprocessor gp = null;
private IGeoProcessorResult gpResult = null;
public GPTool()
{
gp = new Geoprocessor();
gp.OverwriteOutput = true;
}
/// <summary>
/// 释放资源
/// </summary>
public void Release()
{
if (gp != null)
{
ESRI.ArcGIS.ADF.ComReleaser.ReleaseCOMObject(gp);
gp = null;
}
if (gpResult != null)
{
ESRI.ArcGIS.ADF.ComReleaser.ReleaseCOMObject(gpResult);
gpResult = null;
}
}
/// <summary>
/// 执行已定义参数的GP工具。
/// </summary>
protected bool Excute(IGPProcess gpProcess, out string resultInfo)
{
resultInfo = "执行成功";
gpResult = (IGeoProcessorResult)gp.Execute(gpProcess, null);
if (gpResult == null || gpResult.Status != esriJobStatus.esriJobSucceeded)
{
// 执行错误，输入错误报告。
StringBuilder message = new StringBuilder(gpProcess.ToolName);
for (int i = 0; i < gp.MessageCount; i++)
{
message.AppendLine(gp.GetMessage(i));
}
resultInfo = message.ToString();
return false;
}
ESRI.ArcGIS.ADF.ComReleaser.ReleaseCOMObject(gpProcess);
return true;
}
}
}
下边是调用代码：
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using ESRI.ArcGIS.esriSystem;
using GPOperator;
using ESRI.ArcGIS.DataSourcesRaster;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.DataSourcesFile;
using ESRI.ArcGIS.Geometry;
namespace ConsoleApplication1
{
class Program
{
private static LicenseInitializer m_AOLicenseInitializer = new ConsoleApplication1.LicenseInitializer();
static void Main(string[] args)
{
//ESRI License Initializer generated code.
m_AOLicenseInitializer.InitializeApplication(new esriLicenseProductCode[] { esriLicenseProductCode.esriLicenseProductCodeBasic,
esriLicenseProductCode.esriLicenseProductCodeAdvanced },
new esriLicenseExtensionCode[] { esriLicenseExtensionCode.esriLicenseExtensionCodeNetwork, esriLicenseExtensionCode.esriLicenseExtensionCodeSpatialAnalyst,
esriLicenseExtensionCode.esriLicenseExtensionCodeSchematics, esriLicenseExtensionCode.esriLicenseExtensionCodeMLE,
esriLicenseExtensionCode.esriLicenseExtensionCodeDataInteroperability, esriLicenseExtensionCode.esriLicenseExtensionCodeTracking,
esriLicenseExtensionCode.esriLicenseExtensionCodeArcScan, esriLicenseExtensionCode.esriLicenseExtensionCodeBusiness, esriLicenseExtensionCode.esriLicenseExtensionCodeCOGO,
esriLicenseExtensionCode.esriLicenseExtensionCodeGeoStats, esriLicenseExtensionCode.esriLicenseExtensionCodePublisher });
//ESRI License Initializer generated code.
test();
//Do not make any call to ArcObjects after ShutDownApplication()
m_AOLicenseInitializer.ShutdownApplication();
}
public static void test()
{
//准备数据
string polygon = @"D:\quickly\其它\矢量裁剪栅格\shp\LD.shp";
string inputRastger = @"D:\quickly\其它\矢量裁剪栅格\temp\risk11.tif";
ClipHelper c = new ClipHelper();//必须要在打开栅格之前调用，不然打开栅格会出错，还不知道原因。
string errInfo = "";
IRasterDataset rasterDataset = new RasterDataset();
rasterDataset.OpenFromFile(inputRastger);
IWorkspaceFactory wsf = new ShapefileWorkspaceFactory();
IFeatureWorkspace fws = wsf.OpenFromFile(System.IO.Path.GetDirectoryName(polygon), 0) as IFeatureWorkspace;
IFeatureClass f = fws.OpenFeatureClass("LD");
var feat = f.Search(null, false).NextFeature();
var ret = c.ClipRaster(rasterDataset, feat.Shape as IPolygon, out errInfo);
}
}
}