-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathRead NetCDF in batch (ArcGIS)
56 lines (43 loc) · 1.93 KB
/
Read NetCDF in batch (ArcGIS)
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
###################################################################################
######## Read NOAA's NARR meteorolgy data (NetCDF format) in ArcGIS ###############
###################################################################################
import arcpy
# Variable name
variable = "prmsl"
# Longitude
x_dimension = "x"
# Latitude
y_dimension = "y"
band_dimension = ""
dimension = "time"
valueSelectionMethod = "BY_VALUE"
outLoc = "C:/Users/######/AOD/SLP/SLP_out/2012/"
inNetCDF = "C:/Users/######/AOD/SLP/prmsl.2012.nc"
nc_FP = arcpy.NetCDFFileProperties(inNetCDF)
nc_Dim = nc_FP.getDimensions()
for dimension in nc_Dim:
top = nc_FP.getDimensionSize(dimension)
for i in range(365,366):
if dimension == "time":
dimension_values = nc_FP.getDimensionValue(dimension, i)
nowFile = "prmsl"+str(i)
#THIS IS THE NEW CODE HERE
dv1 = ["time", dimension_values]
dimension_values = [dv1]
#END NEW CODE
arcpy.MakeNetCDFRasterLayer_md(inNetCDF, variable, x_dimension, y_dimension, nowFile,band_dimension, dimension_values, valueSelectionMethod)
arcpy.CopyRaster_management(nowFile, outLoc + nowFile + ".img", "", "", "", "NONE", "NONE", "")
print dimension_values, i
# The number of output layers will equal to the number of "dimension"
import arcpy
# Conver raster to vector points
in_path = "C:/Users/######//AOD/SLP/SLP_out/2012/"
out_path = "C:/Users/######//GIS/AOD/SLP/2012_TX/"
# Clip raster
for i in range(365,366):
in_raster = in_path +"prmsl"+str(i)+".img"
out_feature = "prmsl_2012_tx"+str(i)
Clip_feature="C:/Users/xzhang10/Desktop/Infill/GIS/AOD/Boundary/Texas_extend.shp"
arcpy.Clip_management(in_raster,"#",out_feature,Clip_feature,"#","ClippingGeometry","MAINTAIN_EXTENT")
# Conver raster to vector points
arcpy.RasterToPoint_conversion("prmsl_2012_tx"+str(i), out_path+"prmsl_2012_tx"+str(i)+".shp", "VALUE")