python遥感数据有偿处理_Python处理遥感影像总结

由于以前从没有用python处理过遥感影像数据，总觉得很难，这次有个算法用arcgis处理起来实在是过于麻烦，于是进行了一次尝试，没想到居然一点一点的啃下来了。

因为平时python代码写的少，所以就用sublime来编写，调试全靠print大法，虽然也能实现，但是效率要低很多。我寻思用python写小任务确实很方便，索性就着手用VSCode搭建了python开发环境。(在这里安利一波VSCode，真是太好用了，开源、轻量化、跨平台、兼容多种语言、完善的插件、完美集成其他编辑器的快捷键，简直就是开发者的福音)

具体搭建过程我参考了VSCode官方文档及下边这个博客

在搭建过程中涉及到了一个我没听过的新名词：python虚拟环境解释器，说白了就是在当前项目下构建了一个目录，其中放置了简易的python解释环境，这个环境和你本机安装python的环境不同，它非常干净，可能用到的大部分库均需重新下载。

环境搭建好以后就要着手开始编写了，我看了很多博客，大家对于arcgis数据库操作及GIS分析方法大部分都调用了ArcPy库，我就寻思也用这个库试试，然而发现在虚拟环境下没办法安装，原因是国内外各大站点都找不到这个库，通过检索分析，我猜可能是由于esri考虑到吃饭问题没有将arcgis相关的算法库开源出来，毕竟以国人的模仿能力分分钟开发出一个ChinaGIS。

所以用虚拟环境解释器的想法只好作罢，于是我用了arcgis的python解释环境，果然可以用arcpy库中的函数。但是又出现了其他问题，由于pandas库包含了非常广泛的科学计算函数，很好用，我就想调用一下，然而将pandas和arcpy库集成的过程中我遇到了大麻烦。arcgis的python环境是没有pip工具的，因此在线快捷安装就别想了。然后我将anaconda中的pandas库放到python环境下，提示仍然没有找到该模块，但是我用pip去安装的时候命令行提示已经安装好了。偶然之间我发现了下面这一篇博客

究其原因是由于pandas库依赖numpy，而且版本必须相对应，否则就会报错。于是我通过arcgis版本找到了numpy版本，再根据numpy版本找到对应pandas版本，到pandas官网下载对应EXE文件(此处注意x86和x64区别)，安装，顺利解决，over.

接下来就是编码问题，刚开始的时候由于没用过这个库，每写一行代码都需要去检索应该用哪个函数(一点都不夸张)，而且由于不知道如何描述问题，只能乱七八糟的找一些博客，看了半天都和自己的情况有所出入，偶然之下发现了arcgis官方文档(留下链接

在代码里边我做的注释非常详细，所以在这里就不分段落解释了。就将主要的函数列一下吧，方便检索。ListDatasets(wild_card, feature_type) 返回数据集列表

GetRasterProperties_management(in_raster, {property_type}, {band_index}) 获取栅格各种属性，包括波段、像元值等

ExtractByAttributes(InRaster, SQL) 按属性提取

RasterToPolygon_conversion(in_raster, out_polygon_features, {simplify}, {raster_field}) 栅格转面

Dissolve_management(in_features, out_feature_class, {dissolve_field}, {statistics_fields}, {multi_part}, {unsplit_lines}) 融合(Data Management)

AddField_management(in_table, field_name, field_type, {field_precision}, {field_scale}, {field_length}, {field_alias}, {field_is_nullable}, {field_is_required}, {field_domain}) 添加字段

CalculateGeometryAttributes_management (in_features, field, geometry_property, {length_unit}, {area_unit}, {coordinate_system}) 计算几何属性

SearchCursor(in_table, field_names, {where_clause}, {spatial_reference}, {explode_to_points}, {sql_clause}) 游标访问属性数据

# coding:utf-8

import arcpy

from arcpy import env

from arcpy.sa import * # ArcGIS Spatial Analyst

import os

import pandas as pd

gdbPath = "D:/Codemanagement/PythonCode/*****.gdb"

arcpy.env.workspace = gdbPath

featureClasses = []

datasetNum = 0

# ListDatasets(wild_card, feature_type) 返回数据集列表

# wild_card：通配符，限制返回结果

# feature_type: 限制返回结果的数据集类型 Feature(地理数据集数据集) Raster(栅格数据集)

for lightsImg in arcpy.ListDatasets("", "") + [""]:

if not lightsImg == "":

datasetNum = datasetNum + 1

# GetRasterProperties_management(in_raster, {property_type}, {band_index}) 获取栅格各种属性

# in_raster：输入栅格

# {property_type}：属性类型

# {band_index}：波段索引，默认使用第一个波段

# imgBandCount:获取栅格数据波段数

# imgBandCount = arcpy.GetRasterProperties_management(lightsImg, "BANDCOUNT") 数据为单波段 暂不执行

# pixelMax:获取栅格所有像元最大值

pixelMax = int(float(arcpy.GetRasterProperties_management(lightsImg, "MAXIMUM").getOutput(0)))

# pixelMin:获取栅格所有像元最小值

# pixelMin = int(float(arcpy.GetRasterProperties_management(lightsImg, "MINIMUM").getOutput(0))) 自定义最小值，该步骤暂不执行

# print(str(lightsImg) + "Info bandCount:" + str(imgBandCount) + "; pixelMax:" + str(pixelMax) + "; pixelMin:" + str(pixelMin))

print(str(lightsImg) + "Info bandCount: pixelMax:" + str(pixelMax))

# 前期变化较小，自定义亮度值从30开始

def_pixelMin = 30

# 用于存储某年份灯光结果数据

lightTable = []

for lightVal in range(def_pixelMin, pixelMax):

SQLClause = "Value > {0}".format(lightVal)

# ExtractByAttributes(InRaster, SQL) 按属性提取

# InRaster: 输入栅格数据

# SQL: 属性条件

# resPixel存储按属性提取的结果

resPixel = ExtractByAttributes(lightsImg, SQLClause)

outResPixel = "D:/Codemanagement/PythonCode/*****.gdb/{0}_Extra_{1}".format(str(lightsImg),lightVal)

outPixelUrlNum = len(outResPixel.split("/")) - 1

outResPixelName = outResPixel.split("/")[outPixelUrlNum]

if outResPixelName not in arcpy.ListDatasets("", ""):

# arcpy.DeleteFeatures_management(gdbPath + "/" + outResPixelName)

resPixel.save(outResPixel)

# RasterToPolygon_conversion(in_raster, out_polygon_features, {simplify}, {raster_field}) 栅格转面

# in_raster: 输入栅格

# out_polygon_features：输出要素集

# {simplify}： 平滑属性 SIMPLEFY(输出结果面为平滑简单形状); NO_SIMPLIFY(输出结果与栅格像元边缘保持完全一致)

# {raster_field}：将输入栅格字段制定输出数据集的面

# outVecSurf = "D:/Codemanagement/PythonCode/LightsMutDet/outVecSurf/resultTmp.shp"

# 保持工作空间不变，仍在当前gdb下

outVecSurf = "D:/Codemanagement/PythonCode/*****.gdb/{0}_Vec_{1}".format(str(lightsImg),lightVal)

outVecUrlNum = len(outVecSurf.split("/")) - 1

outVecSurfName = outVecSurf.split("/")[outVecUrlNum]

if outVecSurfName not in arcpy.ListDatasets("", ""):

arcpy.RasterToPolygon_conversion(resPixel, outVecSurf, "SIMPLIFY", "Value")

# Dissolve_management(in_features, out_feature_class, {dissolve_field}, {statistics_fields}, {multi_part}, {unsplit_lines}) 融合(Data Management)

# in_features: 输入要素

# out_feature_class：输出要素

vecSurfDissolve = "D:/Codemanagement/PythonCode/*****.gdb/{0}_Dissolve_{1}".format(str(lightsImg),lightVal)

outDisUrlNum = len(vecSurfDissolve.split("/")) - 1

vecSurfDisName = vecSurfDissolve.split("/")[outDisUrlNum]

if vecSurfDisName not in arcpy.ListDatasets('',''):

arcpy.Dissolve_management(outVecSurf, vecSurfDissolve)

# AddField_management(in_table, field_name, field_type, {field_precision}, {field_scale}, {field_length}, {field_alias}, {field_is_nullable}, {field_is_required}, {field_domain}) 添加字段

# in_table: 添加字段的表，包括地理数据库要素类、coverage、shapefile、栅格目录、独立表等

# field_name：添加的字段名

# field_type：新增的字段类型

# 其余参数可选

arcpy.AddField_management(vecSurfDissolve, "Perimeter", "DOUBLE")

# CalculateField_management 字段计算器，，计算周长

arcpy.CalculateField_management(vecSurfDissolve, "Perimeter", "!shape.length!", "PYTHON_9.3")

# CalculateGeometryAttributes_management (in_features, field, geometry_property, {length_unit}, {area_unit}, {coordinate_system}) 计算几何属性

# in_features：要素图层

# field: 将要更新的字段

# geometry_property：要计算的几何属性

# arcpy.CalculateGeometryAttributes_management(vecSurfDissolve, "Perimeter", "DOUBLE", "KILOMETERS") # 模块内不含该函数，暂且屏蔽

# SearchCursor(in_table, field_names, {where_clause}, {spatial_reference}, {explode_to_points}, {sql_clause}) 游标访问属性数据

cursorRows = arcpy.da.SearchCursor(vecSurfDissolve, "Perimeter")

Perimeter = 0

for cursorRow in cursorRows:

# print(cursorRow)

Perimeter = cursorRow[0]

# 每一行的数据，将阈值和周长结果记入

lightRow = []

lightRow.append(lightVal)

lightRow.append(Perimeter)

# 将每一行结果记入表中

lightTable.append(lightRow)

print(lightTable)

field = ["lightVal", "Perimeter"]

resLightTable = pd.DataFrame(columns=field, data = lightTable) #数据为两列

print(resLightTable)

resLightTable.to_csv("D:/Codemanagement/PythonCode/*****/{0}".format(lightsImg))

print("dataSet_Num:" + str(datasetNum))

从搭开发环境到最终出结果整个过程花了三四天的时间吧，其实作出成果的时候成就感还是非常强烈的，这个事情从零开始学习还是具有一定的挑战性，但是人生不就是这样从无到有积累经验吗？

在这里想表达一个想法：工作中遇到最大的困难其实就是自己的恐惧心理！共勉。