GLDAS/dist_average_temp.py at master · IFPRI/GLDAS · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124

import arcpy
from arcpy.sa import *
import glob, os, time
import multiprocessing

arcpy.CheckOutExtension("Spatial")
arcpy.env.workspace = "D:/Users/cmarciniak/tables"
arcpy.env.overwriteOutput = True


outputFolder = "D:/Users/cmarciniak/My Documents/GLDAS/rasters/"
tables = "D:/Users/cmarciniak/My Documents/GLDAS/tables/"


# takes in an array of rasters, outputs one file per day with the average surface temperature
def dailyAverage(rasters , filename):
	outCellStats = CellStatistics( rasters , "MEAN","DATA" )
	outCellStats.save("D:/Users/cmarciniak/My Documents/GLDAS/rasters/"+filename)

def makeRasterArray(files):
	#outputFolder = "D:/Users/cmarciniak/My Documents/GLDAS/rasters/"
	print("Making rasters...")
	rasters=[]
	for file in files:
		print("Making raster from "+file)
		tempFile = file[0:12]
		arcpy.MakeNetCDFRasterLayer_md( file , "AvgSurfT_GDS0_SFC", "g0_lon_1","g0_lat_0",tempFile)
		rasters.append(tempFile)
	return rasters


def cellStats(rasters):
	print(rasters)
	outputFolder = "D:/Users/cmarciniak/My Documents/GLDAS/rasters/"
	outCellStats = CellStatistics(rasters,"MEAN","DATA")
	try:
		outCellStats.save("in_memory/avg_temp")
	except RuntimeError:
		#os.rmdir(outputFolder+"avg_temp")
		outCellStats.save("in_memory/avg_temp")
	#arcpy.Delete_management("in_memory")

def deleteCellStats(rasters):
	for raster in rasters:
		arcpy.Delete_management(raster)

def resample(inRaster,outRaster):
	if os.path.isdir(outRaster):
		print("Removing "+outRaster)
		#os.rm("D:/Users/cmarciniak/My Documents/GLDAS/rasters/resample.aux.xml")
		#shutil.rmtree(outRaster)
		arcpy.Resample_management(inRaster,outRaster,.05)
	else:
		print("Resampling...")
		arcpy.Resample_management(inRaster,outRaster,.05)
	print("Done")

def ncFileIterator():
	files = glob.glob("*.nc")
	files.sort()
	rasters = []
	temp = []
	for i in range ( 0,len(files)-1 ):
		if files[i][0:7] == files[i+1][0:7]:
			print(files[i][0:7])
			temp.append(files[i])
		else:
			print("Done")
			temp.append(files[i])
			rasters.append(temp)
			temp = []
	return rasters


def mainLoop(ncFiles):
	districts = "D:/Users/cmarciniak/My Documents/ArcGIS/2011 Districts MAPS Shape Files/DISTRCT_BDY_geog.shp"
	#ncFiles = ncFileIterator()
	rasters = []
	for nc in ncFiles:
		day = nc[0][0:7]
		if(os.path.isfile("D:/Users/cmarciniak/tables/GLDAS_"+day+".dbf")):
			print("Skipping "+day+".dbf")
			pass
		else:
			rasters = makeRasterArray(nc)
			cellStats(rasters)
			print(outputFolder)
			resample("in_memory/avg_temp","in_memory/resample")
			outZonalStatistics = ZonalStatisticsAsTable(districts,"ID", "in_memory/resample","D:/Users/cmarciniak/tables/GLDAS_"+nc[0][0:7]+".dbf")
			deleteCellStats(rasters)
			arcpy.Delete_management("in_memory")
			#arcpy.Delete_management(outZonalStatistics)
			#del rasters

def f(ncFiles):
	districts = "D:/Users/cmarciniak/My Documents/ArcGIS/2011 Districts MAPS Shape Files/DISTRCT_BDY_geog.shp"
	rasters = makeRasterArray(ncFiles)
	cellStats(rasters)
	resample("in_memory/avg_temp","in_memory/resample")
	outZonalStatistics = ZonalStatisticsAsTable(districts,"ID", "in_memory/resample","D:/Users/cmarciniak/tables/GLDAS_"+ncFiles[0][0:7]+".dbf")


ncFiles = ncFileIterator()

'''if __name__== '__main__':
	pool = multiprocessing.Pool(processes=4)
	result = pool.map(f,ncFiles,4)
	pool.close()
	pool.join()
'''
mainLoop(ncFiles)

#rasters = makeRasterArray(files)
#cellStats(rasters)
#print(outputFolder+"avg_temp/w001001.adf")
#arcpy.env.workspace = ""
#outZonalStatistics = ZonalStatisticsAsTable(districts,"ID", "D:/Users/cmarciniak/resample4","D:/Users/cmarciniak/zones4.dbf")
#outZonalStatistics.save("D:/Users/cmarciniak/My Documents/GLDAS/zones")