Siqi Dong's Blog

# -*- coding: utf-8 -*-
import os

def addAsciiHeadFile(asciiFolder):
    rootdir = os.path.join(asciiFolder)
    for (dirpath, dirnames, filenames) in os.walk(rootdir):
        for filename in filenames:
            if os.path.splitext(filename)[1] == '.txt':
                asciiPath=asciiFolder+"\\"+filename
                with open(asciiPath, 'r+') as f:
                    content = f.read()
                    f.seek(0, 0)
                    # head file content
                    f.write('XXX\n' + 'XXX\n' + 'XXX\n' + 'XXX\n' + 'XXX\n' + 'XXX\n' + content)
                    print asciiPath + "done!"

# -*- coding: utf-8 -*-
import arcpy
import os

arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True
def mosaicToNewRaster(inputFolder,outputFolder,tempfolder,outputname):
    infolder = inputFolder + "\\"
    outfolder = outputFolder + "\\"
    tempfolder = tempfolder + "\\"
    rasters = []
    datanames = os.listdir(infolder)
    for dataname in datanames:
        if os.path.splitext(dataname)[1] == '.tif':
            print  dataname
            inputfullname = infolder + dataname
            print "inputfullname: " + inputfullname
            raster = arcpy.Raster(inputfullname)
            # remove black border
            deldarkfullname = tempfolder + "temp" + dataname
            arcpy.CopyRaster_management(raster, deldarkfullname, "DEFAULTS", 0, 0, "", "", "8_BIT_UNSIGNED")
            raster2 = arcpy.Raster(deldarkfullname)
            rasters.append(raster2)
    # mosaic raster
    arcpy.MosaicToNewRaster_management(rasters, outfolder, outputname, "", "8_BIT_UNSIGNED", "30", "1", "MINIMUM", "FIRST")
    print("Mosaic success")

import arcpy
import os

def  GDBConvertToShp(GDBfolder ,SHPfolder):
    for root, dirs, files in os.walk(GDBfolder):
        for d in dirs:
            path= root+"\\"+d
            if path[-4:]==".gdb":
                name=path.split(".")[0].split("\\")[-1]
                arcpy.env.workspace =path
                gdbshp=arcpy.ListFeatureClasses()
                for fc in gdbshp:
                    # the name of the feature in GDB you want to convert to shpfile
                    if str(fc)=="XXXX":
                        arcpy.FeatureClassToFeatureClass_conversion(fc, SHPfolder,name+str(fc)+".shp")

# -*- coding: utf-8 -*-
import arcpy
import  os
from arcpy.sa import *
import envipy
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True

name=["h23v04","h24v04","h25v04","h23v05","h24v05","h25v05","h26v05","h25v06","h26v06","h27v06"]
# parameter
path="E:\\MODIS"
shp="E:\\MODIS\\tibetALB_Dissolve.shp"
Sprj="E:\\MODIS\\Sinusoidal.prj"
Albprj="E:\\MODIS\\Albers_Conic_Equal_Area.prj"

for i in range(len(name)):
    img=name[i]
    imgarray=[]
    f_list = os.listdir(path)
    for file in f_list:
        filename = path + "\\" + file
        if os.path.splitext(file)[1] == '.hdf' and filename.split(".")[2][0:6] == img:
            hdf = path + "\\" + file
            tif = path+"\\"+os.path.splitext(file)[0] + ".tif"
            print hdf
            print tif
            # Extract 0-band subdataset from HDF
            arcpy.ExtractSubDataset_management(hdf, tif, "0")
            # DefineProjection
            arcpy.DefineProjection_management(tif, Sprj)
            imgarray.append(tif)
    outCellStatistics = CellStatistics(imgarray, "MEAN", "DATA")
    # Save the output
    outCellStatistics.save("E:\\MODIS\\result\\"+img+"_2015NDVI_MEAN.tif")
f_list = os.listdir("E:\\MODIS\\result")
rasterarray = []
for file in f_list:
    filename = path + "\\result\\" + file
    if os.path.splitext(file)[1] == ".tif":
        rasterarray.append(filename)
print rasterarray
arcpy.MosaicToNewRaster_management(rasterarray,"E:\\MODIS\\result","NDVI2015_m.tif","", "16_BIT_SIGNED", "", "1", "LAST", "FIRST")
# Project transfer
arcpy.ProjectRaster_management("E:\\MODIS\\result\\NDVI2015_m.tif", "E:\\MODIS\\result\\NDVI2015_alb.tif", Albprj, "","250", "","","","")

# Extract
# outExtractByMask = ExtractByMask("E:\\MODIS\\result\\NDVI2015_alb.tif", shp)
# Save the output
# outExtractByMask.save("E:\\MODIS\\result\\Tibet_NDVI2015_int_alb.tif")
# intRaster=arcpy.Raster("E:\\MODIS\\result\\Tibet_NDVI2015_int_alb.tif")

import arcpy
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")

def getValueSum(rasterPath):
    raster=arcpy.Raster(rasterPath)
    raster2 = SetNull(raster == 0,raster)
    sum=0
    rstArray = arcpy.RasterToNumPyArray(raster2)
    rows, cols = rstArray.shape
    for rowNum in xrange(rows):
        for colNum in xrange(cols):
            s=rstArray.item(rowNum, colNum)
            sum=sum+s
    return sum