12章代码

12.1 栅格基本信息获得

#coding=utf8
import arcpy

import os
import sys
import math


def getbandinfo(band1):
    desc = arcpy.Describe(band1)
    if hasattr(desc,"height"):
        arcpy.AddMessage("Height: %d" % desc.height)
    if hasattr(desc,"width"):
        arcpy.AddMessage("Width:  %d" % desc.width)
    if hasattr(desc,"isInteger"):
        arcpy.AddMessage("Integer Raster: %s" % desc.isInteger)
    if hasattr(desc,"noDataValue"):
        arcpy.AddMessage("noDataValue: %s" % desc.noDataValue)

    if hasattr(desc,"meanCellHeight"):
        arcpy.AddMessage(u"Y方向分辨率: %s" % desc.meanCellHeight)
    if hasattr(desc,"meanCellWidth"):
        arcpy.AddMessage(u"X方向分辨率: %s" % desc.meanCellWidth)
    if hasattr(desc,"pixelType"):
        arcpy.AddMessage("pixelType: %s" % desc.pixelType)
    if hasattr(desc,"primaryField"):
        arcpy.AddMessage("primaryField: %d" % desc.primaryField)
    if hasattr(desc,"tableType"):
        arcpy.AddMessage("tableType: %s" % desc.tableType)

##    U1 —1 bit
##    U2 —2 bits
##    U4 —4 bits
##    U8 —Unsigned 8 bit integers
##    S8 —8 bit integers
##    U16 —Unsigned 16 bit integers
##    S16 —16 bit integers
##    U32 —Unsigned 32 bit integers
##    S32 —32 bit integers
##    F32 —Single precision floating point
##    F64 —Double precision floating point


#获得栅格数据集的信息
def getRasterDataset():
    desc = arcpy.Describe(sRaster)

    arcpy.AddMessage("Band Count:       %d" % desc.bandCount)
    arcpy.AddMessage("Compression Type: %s" % desc.compressionType)
    arcpy.AddMessage("Raster Format:    %s" % desc.format)
    arcpy.AddMessage("Raster sensorType:    %s" % desc.sensorType)
    arcpy.AddMessage("Raster permanent: " + str(desc.permanent))

    my_raster = arcpy.Raster(sRaster)
    arcpy.AddMessage("Raster maximum: " + str(my_raster.maximum))
    arcpy.AddMessage("Raster minimum: " + str(my_raster.minimum))
    arcpy.AddMessage("Raster mean: " + str(my_raster.mean))
    arcpy.AddMessage("Raster name: " + my_raster.name)
    arcpy.AddMessage("Raster path: " + my_raster.path)

    band1=sRaster+os.sep+"Band_1"
    if not arcpy.Exists(band1):
        band1 = sRaster + os.sep + "Layer_1"  #在mdb中是Layer_1
        if not arcpy.Exists(band1):
            return



    getbandinfo(band1)
    CELLSIZEX=arcpy.GetRasterProperties_management(sRaster, "CELLSIZEX") #分辨率
    arcpy.AddMessage("CELLSIZEX: " + str(CELLSIZEX))

def Main():
    getRasterDataset()


sRaster  = arcpy.GetParameterAsText(0)

Main()

12.2   栅格计算

def Main2():#第一种方法
    outRas = (arcpy.Raster(sRaster) > myValue-0.01) & ( arcpy.Raster(sRaster) < myValue+0.01)
    myCon=arcpy.sa.Con(outRas,sRaster)
    try:
        arcpy.RasterToPoint_conversion(myCon,outPoint,raster_field="Value")
    except Exception as e:
        arcpy.AddError(e)
        arcpy.AddMessage(u"没有对应数据".encode('gbk')) #汉字乱码的解决
def Main():#第二种方法
    yl999="in_memory/yl999"
    arcpy.gp.RasterCalculator_sa('Con(Abs("'+sRaster+'" - '+str(myValue)+')<0.01,"'+sRaster+'")',yl999)
    arcpy.RasterToPoint_conversion(yl999,outPoint,raster_field="Value")

sRaster  = arcpy.GetParameterAsText(0)
myValue  = arcpy.GetParameter(1)
outPoint=arcpy.GetParameterAsText(2)
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput=True

arcpy.env.extent = "MAXOF"
Main2()

12.3 影像合并

#coding=utf8
import arcpy
import arcpy.sa
import os
import sys
import math

sdir  = arcpy.GetParameterAsText(0)
outRaster  = arcpy.GetParameterAsText(1)
bandnum=arcpy.GetParameter(2)
arcpy.CheckOutExtension("Spatial")

arcpy.env.extent = "MAXOF"
arcpy.env.workspace = sdir

rasters = arcpy.ListRasters()
num=len(rasters)
i=1
dras=""
for raster in rasters:
    dras=dras+raster
    if i<num:
        dras=dras+";"

    i=i+1
a,b=os.path.split(outRaster)

arcpy.MosaicToNewRaster_management(dras, a, 
                                   b, "",
                                   "8_BIT_UNSIGNED", "",str(bandnum) , "LAST","FIRST")

12.4  多个栅格统计

sRaster  = arcpy.GetParameterAsText(0)
tjlx  = arcpy.GetParameterAsText(1) #统计类型
outRaster  = arcpy.GetParameterAsText(2)
arcpy.CheckOutExtension("Spatial")
inRaster=sRaster.split(";")

outCellStats = arcpy.sa.CellStatistics(inRaster, tjlx, "NODATA")
outCellStats.save(outRaster)

12.5  批量格式转换

def FormatConvert():
    arcpy.env.workspace = sdir
    rasters = arcpy.ListRasters()
    num=len(rasters)
    i=1
    dras=""
    for raster in rasters:
        dras=dras+raster
        if i<num:
            dras=dras+";"
        i=i+1
    Raster_Format="TIFF"
    if picFormat==".tif":
        Raster_Format="TIFF"
    elif picFormat==".png":
        Raster_Format="PNG"
    elif picFormat==".img":
        Raster_Format="IMAGINE Image"
    elif picFormat==".jpg":
        Raster_Format="JPEG"

    try:
        arcpy.RasterToOtherFormat_conversion(dras, outdir,Raster_Format)
    except Exception as e:
        arcpy.AddError(e)

sdir  = arcpy.GetParameterAsText(0)
outdir  = arcpy.GetParameterAsText(1)
picFormat=arcpy.GetParameterAsText(2).lower()
arcpy.CheckOutExtension("Spatial")

arcpy.env.extent = "MAXOF"
FormatConvert()
#使用复制栅格 第二种方法
def FormatConvert2():
    arcpy.env.workspace = sdir
    rasters = arcpy.ListRasters()
    num=len(rasters)
    initProgress(u"正在更新",num)

    for raster in rasters:
        step()
        try:
            file_name = os.path.basename(raster)
            #print(file_name)
            # 输出为 test.py
            file_name = file_name.split('.')[0]
            #print(file_name),不加扩展名
            #arcpy.AddMessage("file_name: " + file_name)
            out_Raster=outdir+os.sep+file_name+picFormat
            arcpy.AddMessage("outRaster: " + out_Raster)

            arcpy.CopyRaster_management(raster,out_Raster)
        except Exception as e:
            arcpy.AddError(e)

    freeProgress()

12.6 栅格彩色转黑白

1.numpy使用
def copy8(in_ras,outras):
    arcpy.CopyRaster_management(in_raster=in_ras,out_rasterdataset=outras,pixel_type="8_BIT_UNSIGNED")
    desc=arcpy.Describe(outras)
    sr=desc.spatialReference
    if sr.name == "Unknown":
        desc=arcpy.Describe(inFeature)
        sr=desc.spatialReference
        arcpy.DefineProjection_management(outras, sr)

def colortoBlack():
    CELLSIZEX=arcpy.GetRasterProperties_management(inFeature, "CELLSIZEX") #分辨率
    cellsize=float(CELLSIZEX.getOutput(0))
    X=arcpy.GetRasterProperties_management(inFeature, "LEFT") #LEFT
    Y=arcpy.GetRasterProperties_management(inFeature, "BOTTOM") #TOP,BOTTOM
    x=float( X.getOutput(0))
    y=float( Y.getOutput(0))

    n=arcpy.RasterToNumPyArray(inFeature)

    m = n[0] * 0.299 + n[1] * 0.587 + n[2] * 0.114
    #mint=numpy.trunc(m) #取整数
    point = arcpy.Point(x, y)

    my_raster = arcpy.NumPyArrayToRaster(m,point,cellsize,cellsize)
    copy8(my_raster,outFeature)
2.使用计算器
def Main():

    desc = arcpy.Describe(inFeature)
    filepath=desc.catalogPath
    #arcpy.AddMessage(u"filepath:" + filepath)
    if desc.bandCount!=3:
        arcpy.AddMessage(u"数据:" + inFeature+"波段不为3")
        return
    myRaster1 = arcpy.Raster(filepath+"/Band_1")
    myRaster2 = arcpy.Raster(filepath+"/Band_2")
    myRaster3 = arcpy.Raster(filepath+"/Band_3")

    myRaster=myRaster1*0.299 + myRaster2*0.587 + myRaster1*0.114
    #myRaster.save(outFeature)
    copy8(myRaster,outFeature)
    arcpy.BuildPyramids_management(outFeature)
    mxd = arcpy.mapping.MapDocument("CURRENT")
    df = arcpy.mapping.ListDataFrames(mxd)[0]

    addLayer = arcpy.mapping.Layer(outFeature)
    arcpy.mapping.AddLayer(df, addLayer, "AUTO_ARRANGE")

inFeature  = arcpy.GetParameterAsText(0)
outFeature  = arcpy.GetParameterAsText(1)
arcpy.CheckOutExtension("Spatial")
#colortoBlack()
Main()

12.7 栅格重分类

def Main2():
    arcpy.gp.Reclassify_sa(sRaster,"Value","0 2 1;2 6 2;6 15 3;15 25 4;25 90 5",outRaster,"DATA")

def Main():
    outReclass = arcpy.sa.Reclassify(sRaster, "Value", arcpy.sa.RemapRange([[0,2,1],[2,6,2],[6,15,3], [15,25,4] ,[25,90,5]]))
    outReclass.save(outRaster)

sRaster  = arcpy.GetParameterAsText(0)
outRaster  = arcpy.GetParameterAsText(1)
arcpy.CheckOutExtension("Spatial")
arcpy.CheckOutExtension("3D")
arcpy.env.overwriteOutput=True

arcpy.env.extent = "MAXOF"
#Main()
Main2()