使用矢量裁剪栅格,并统计栅格的和(排除noda与小于0的值)
# -*- coding: utf-8 -*- import arcpy import os import sys import argparse from osgeo import gdal import xlsxwriter # 导入模块 def clip(shp,raster): volume={} cursor = arcpy.da.SearchCursor(shp,["OID@","SHAPE@","XZQHDM"]) for row in cursor: mask = row[1] # 让每个斑块成为掩膜 outraster = arcpy.sa.ExtractByMask(raster, mask) im_volume = cal_volume(outraster) volume[row[2]]=im_volume # try: # outraster = arcpy.sa.ExtractByMask(raster, mask) # 根据掩膜提取栅格,就是裁剪 # # # im_volume=cal_volume(outraster) # volume[cursor[row]]=im_volume # # # MEAN = arcpy.GetRasterProperties_management(outraster, "MEAN") # 读取平均值 # # except: # x.append(np.nan) # 如果报错则在列表中添加一个nan值,不然FID值会乱掉,不利于后面的相关分析 # print 'blank part' # 这个部分报错主要是因为有的斑块与栅格数据没有相交的部分,因而就会没法裁剪,从而出错 del row del cursor def cal_volume(outraster): no_data_value = outraster.noDataValue arr = arcpy.RasterToNumPyArray(outraster, nodata_to_value=no_data_value) row_num, col_num = arr.shape sum=0 count=0 pixelsize=2 for i in range(0, row_num): for j in range(0, col_num): if arr[i][j] != no_data_value and arr[i][j] > 0: # 找到相应的值之后处理,我这里是设为NoData值,也可以使用算法来修补 count=count+1 sum=sum+arr[i][j] volume=sum*pixelsize*pixelsize return volume |