Python提取常年存在的作物种植区格点并转换为shp文件导出

问题概述：现有多年的作物种植区格点Tif，想筛选多年一直存在的格点用于后期数据的掩膜处理；有种植区为1，无为0；

1.读取多年栅格数据后，通过numpy数组计算多年栅格点的和，可以将和大于某一数值的点认为是多年一直存在的格点；【其中存在一个疑惑，本身的Tif文件是行政区边界分明的，计算后是新行政区的外接最小矩形框，在矩形框内行政区外，原本没有数据的地方最终返回的值是1？原本希望对于最终满足求和条件返回的值的点直接导出shp，因为原本没有数据的地方最终返回的值也是1，考虑其他方法】

import numpy as np
import pandas as pd
from osgeo import gdal
import glob
import os
input_folder = 'E:/seasonyield_predict/datapre/mask_region/clip_db/'
search_term="*wheat*.tif"
tif_files = glob.glob(os.path.join(input_folder, search_term))

arrays = []
for tif_file in tif_files:
    ds = gdal.Open(tif_file)
    #num_bands = ds.RasterCount查看波段
    band = ds.GetRasterBand(1)
    arr = band.ReadAsArray()
    arrays.append(arr)   
stack = np.stack(arrays, axis=-1)
stack.shape
####5年中都有栅格点的位置
stack_sum=np.sum(stack,axis=2)
result = np.zeros_like(stack_sum)
result[stack_sum >= 5] = 1

2.将计算结果导出为tif格式；

# 将结果写入栅格，读取原始栅格文件，获取栅格信息
from osgeo import gdal, gdal_array, osr
###用ds读到的tif信息
geo_transform = ds.GetGeoTransform()
projection = ds.GetProjection()
driver = gdal.GetDriverByName("GTiff")
#driver = ds.GetDriver()
#colss = ds.RasterXSize#1065
#rowss = ds.RasterYSize##1805
# 创建输出栅格文件
output_file = 'E:/seasonyield_predict/datapre/mask_region/cliped_db_average/sprwheat5.tif'
rows,cols = result.shape
out_dataset = driver.Create(output_file, cols, rows, 1, gdal.GDT_Int16)
out_dataset.SetGeoTransform(geo_transform)
out_dataset.SetProjection(projection)

# 将结果数组写入栅格文件中
band = out_dataset.GetRasterBand(1) #获取波段i+1的地址，波段计数从1起
band.WriteArray(result,0,0) #写入数据
del(out_dataset)

3.因为外接矩形框的存在，tif格式再次按照行政区裁剪后筛选出1值栅格点，再转换为矢量：

python里没有操作出来，转arcpy操作（再次感叹，还是arcpy好用，代码如何写可以在工具帮助中找到，自己只需要加循环即可，唯一不方便的是没有脚本，只能一行一行输命令）：

import os
import arcpy
from arcpy.sa import *
arcpy.env.workspace="E:/seasonyield_predict/datapre/mask_region/cliped_db_average"
rasters=arcpy.ListRasters("*",'tif')#文件下的tif文件
mask=r'E:\seasonyield_predict\datapre\region\clip\clip_dbc.shp'#掩膜区域
outpathm=r'E:\seasonyield_predict\datapre\mask_region\cliped_crop_average'#提取属性值后的输出路径
outpaths=r'E:\seasonyield_predict\datapre\mask_region\cliped_crop_only'#转为矢量的输出路径


for raster in rasters:
   arcpy.CheckOutExtension("Spatial")
   tifm=ExtractByMask(raster, mask)
   outm=outpathm+'\\'+raster.split(".")[0]+".tif"
   crop_abs=arcpy.gp.ExtractByAttributes_sa(tifm, '"Value" =1', out)
   outs=outpaths+'\\'+raster.split(".")[0]+".shp"
   arcpy.RasterToPolygon_conversion(crop_abs, outs, "NO_SIMPLIFY","Value")

直接按属性值提取的栅格数据也可以用于裁剪其它因子的tif, 之前遇到的一个问题，类似于分类的栅格数据，比如有耕地面积的栅格是1，否则为0，我是想提取出有耕地面积的栅格以此来做掩膜提取其他数据；但是该栅格数据按属性值1提取后，它显示出来的确实是只有耕地面积的区域，但是用该栅格数据没法提取别的数据，提取的时候它还是包括了原本所有类型的范围，所以也是只能按属性提取后转为shp；

但是这次按栅格属性提取后可以掩膜提取其它数据，也没必要转为shp；

Python提取常年存在的作物种植区格点并转换为shp文件导出

猜你喜欢