# Author： zhb
# Email: zhanghongbo@itpcas.ac.cn
# Task: to use Multiday Auqa and Terra modis data to extract snow cover
import arcpy
import os
from arcpy import env
from arcpy.sa import *
import arcgisscripting
import re
# mostPath: path of Terra
mostPath = “G:\MODIS\\TerraOut”
# mosaPath: path of Aqua
mosaPath = “G:\MODIS\\AquaOut”
Tfns = os.listdir(mostPath)
Afns = os.listdir(mosaPath)
woPath =”G:\MODIS\wo7″
totalSize = len(Tfns)
outFile = open(woPath+”\\”+”result.txt”, ‘w’)
outFile.write(“日期,”+”雪盖,”+”云盖,”+”云量,”+”最低高程”+’\n’)
env.workspace = woPath
env.extent = “G:\MODIS\WO2\mdomin.shp”
gp=arcgisscripting.create(10.0)
ma=gp.MultiOutputMapAlgebra_sa
mask = “G:\卡曲数据准备\kalu_bound.shp”
dem = “dem”
# sini: 是否第一次做雪盖统计
# cini: 是否第一次做云盖统计
# sfailed: 雪盖栅格数据是否不存在
# cfailed: 云盖栅格数据是否不存在
# maxDays: 最大合并日期
# multiCount: 多日整合计数器
# isCombine: 是否整合多日数据
istart = 0
sini = True
cini = True
sfailed = False
cfailed = False
count = 1
sot = “smerge.dbf”
cot = “cmerge.dbf”
maxDays = 1
multiCount = 1
isCombine = True
totalArea = 296487900
startDay = “”
cloudThreshold = 0.1
cloudRate = 0.0
special = False
Tfns.sort()
for tfn in Tfns :
special = False
if multiCount == 1 :
isCombine = False
else :
isCombine = True
i = tfn.find(“MOD10A1.A”)
if (i == -1 or (not tfn.endswith(“.tif”))):
count = count+1
continue
print tfn + ” is processing”
i = i + len(“MOD10A1.A”)
fflag = tfn[i:i+7]
# 在mosa目录中找匹配文件
isFound = False
for afn in Afns :
ia = afn.find(“MYD10A1.A”)
if ia == -1 :
continue
ia = ia + 9
at = afn[i:i+7]
if at == fflag :
isFound = True
break
if not isFound :
print “AquaData for “+tfn+” is not found!”
count = count+1
continue
name = tfn[i+2:i+7]
iname = int(name)
# 预处理Terra数据
TRasterIn = mostPath+”\\”+tfn
TRasterOut = woPath+”\\”+”T”+fflag
tmp = woPath + “\\” + “linshi”
arcpy.CheckOutExtension(“Spatial”)
out = ExtractByMask(TRasterIn, mask)
out.save(tmp)
tmp1 = woPath + “\\” + “linshi1″
arcpy.Resample_management(tmp, tmp1, “30”, “NEAREST”)
out2 = ExtractByMask(tmp1, mask)
out2.save(TRasterOut)
arcpy.Delete_management(tmp)
arcpy.Delete_management(tmp1)
Tsname = “Tsnow”+fflag
Tcname = “Tcloud”+fflag
env.extent = TRasterOut
ma(Tsname+”=con(T”+fflag+” == 200 , “+str(iname)+”)”)
gp=arcgisscripting.create(10.0)
ma=gp.MultiOutputMapAlgebra_sa
ma(Tcname+”=con(T”+fflag+” == 50 , “+str(iname)+”)”)
gp=arcgisscripting.create(10.0)
ma=gp.MultiOutputMapAlgebra_sa
# 预处理Aqua数据
ARasterIn = mosaPath+”\\”+afn
ARasterOut = woPath+”\\”+”A”+fflag
tmp = woPath + “\\” + “linshi”
arcpy.CheckOutExtension(“Spatial”)
out = ExtractByMask(ARasterIn, mask)
out.save(tmp)
tmp1 = woPath + “\\” + “linshi1″
arcpy.Resample_management(tmp, tmp1, “30”, “NEAREST”)
out2 = ExtractByMask(tmp1, mask)
out2.save(ARasterOut)
arcpy.Delete_management(tmp)
arcpy.Delete_management(tmp1)
Asname = “Asnow”+fflag
Acname = “Acloud”+fflag
env.extent = ARasterOut
ma(Asname+”=con(A”+fflag+” == 200 , “+str(iname)+”)”)
gp=arcgisscripting.create(10.0)
ma=gp.MultiOutputMapAlgebra_sa
ma(Acname+”=con(A”+fflag+” == 50 , “+str(iname)+”)”)
gp=arcgisscripting.create(10.0)
ma=gp.MultiOutputMapAlgebra_sa
# 提取雪盖和云盖
sname = “snow”+fflag
cname = “cloud”+fflag
env.extent = TRasterOut
ma(sname+”=con((not(IsNull(“+Tsname+”))) | (not(IsNull(“+Asname+”))),”+str(iname)+”)”)
gp=arcgisscripting.create(10.0)
ma=gp.MultiOutputMapAlgebra_sa
ma(cname+”=con((not(IsNull(“+Tcname+”))) & (not(IsNull(“+Acname+”))),”+str(iname)+”)”)
gp=arcgisscripting.create(10.0)
ma=gp.MultiOutputMapAlgebra_sa
if multiCount == 1:
snowCombine1 = sname
cloudCombine1 = cname
startDay = fflag
elif multiCount == 2:
snowCombine = snowCombine1
cloudCombine = cloudCombine1
snowCombine1 = “tmpsnow1″
cloudCombine1 = “tmpcloud1″
else :
snowCombine = “tmpsnow”
cloudCombine = “tmpcloud”
snowCombine1 = “tmpsnow1″
cloudCombine1 = “tmpcloud1″
# 合并提取
if isCombine :
ma(snowCombine1+”=con((not(IsNull(“+sname+”))) | (not(IsNull(“+snowCombine+”))),”+str(int(startDay[2:7]))+”)”)
gp=arcgisscripting.create(10.0)
ma=gp.MultiOutputMapAlgebra_sa
ma(cloudCombine1+”=con((not(IsNull(“+cname+”))) & (not(IsNull(“+cloudCombine+”))),”+str(int(startDay[2:7]))+”)”)
gp=arcgisscripting.create(10.0)
ma=gp.MultiOutputMapAlgebra_sa
# 计算面积
sizeX = 30.0
sizeY = 30.0
ncount = 0
try :
rows = arcpy.SearchCursor(cloudCombine1,””,””,”COUNT”,””)
row = rows.next();
ncount = row.getValue(“COUNT”)
rows = None
row = None
except :
rows = None
row = None
ncount = 0
cloudArea = ncount * sizeX * sizeY
cloudRate = cloudArea / totalArea
print ‘合成天数： ‘+ str(multiCount) + “, 云量: “+ str(cloudRate)
if( multiCount != 1) :
arcpy.Delete_management(snowCombine)
arcpy.Delete_management(cloudCombine)
arcpy.CopyRaster_management(snowCombine1, “tmpsnow”)
arcpy.Delete_management(snowCombine1)
arcpy.CopyRaster_management(cloudCombine1, “tmpcloud”)
arcpy.Delete_management(cloudCombine1)
# 判断输出
if (multiCount == maxDays or cloudRate <= cloudThreshold or count == totalSize) :
if multiCount != 1 :
arcpy.Copy_management(“tmpsnow”, “xsnow” + startDay)
arcpy.Copy_management(“tmpcloud”, “xcloud” + startDay)
arcpy.Delete_management(“tmpsnow”)
arcpy.Delete_management(“tmpcloud”)
else :
arcpy.Copy_management(sname, “xsnow” + startDay)
arcpy.Copy_management(cname, “xcloud” + startDay)
arcpy.Delete_management(sname)
arcpy.Delete_management(cname)
special = True
try :
rows = arcpy.SearchCursor(“xsnow” + startDay,””,””,”COUNT”,””)
row = rows.next()
ncount = row.getValue(“COUNT”)
rows = None
row = None
except :
rows = None
row = None
ncount = 0
snowArea = ncount * sizeX * sizeY
gp=arcgisscripting.create(10.0)
ma=gp.MultiOutputMapAlgebra_sa
ma(“sdem”+startDay+”=con((not(IsNull(“+”xsnow” + startDay+”))),”+dem+”)”)
gp=arcgisscripting.create(10.0)
ma=gp.MultiOutputMapAlgebra_sa
try :
rows = arcpy.SearchCursor(“sdem” + startDay,””,””,”VALUE”,””)
i0 = (rows.next()).getValue(“VALUE”)
for row in rows :
i1 = row.getValue(“VALUE”)
if( i1 < i0) :
i0 = i1
minElevation = i0
rows = None
row = None
except :
minElevation = 0
outFile.write(startDay + ‘,’ + str(snowArea) + “,” + str(cloudArea) + “,” + str(cloudRate) + “,” + str(minElevation) + ‘\n’)
multiCount = 0
arcpy.Delete_management( Tsname )
arcpy.Delete_management( Tcname )
arcpy.Delete_management( Asname )
arcpy.Delete_management( Acname )
arcpy.Delete_management( “T”+fflag)
arcpy.Delete_management( “A”+fflag )
multiCount = multiCount + 1
print fflag + ” is finished. “+”achieved : “+ str(count) + “/”+ str(totalSize)
count = count+1
outFile.close()
print ‘Finished !’