-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathsubsetGlobalGrids.py
119 lines (90 loc) · 6.02 KB
/
subsetGlobalGrids.py
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
# import libraries
from mbgw import master_grid
import tables as tb
from numpy import *
# set some parameters
#resRatio=5 # ratio between resolution of grids defines in master_grid, and the ones being subsetted here
datafolder = "/home/pwg/mbg-world/datafiles/auxiliary_data/GridsForCS/"
################################################################################################################################################
# using row/col numbers of a subset specified for 5km grids (and stored in master_grid), extract the same subset for corresponding 1km grids
def subsetGlobalGrids (region,lims,gridname,resRatio = 1):
# build input path for the global grid
input_path = datafolder+gridname+".hdf5"
# build ouput paths for this region and grid
output_path = datafolder+gridname+"_"+region+".hdf5"
# open link to full sized 1km file
fullHDF5 = tb.openFile(input_path, mode = "r")
# get 1km row/cols for this subset from the pre-defined object in mbgw 'master_grid' (NB - the values therein were entered manually after running the R script DefineSubsetCoordinates.R)
br=(lims['bottomRow'])*resRatio
tr=(lims['topRow']-1)*resRatio
lc=(lims['leftCol']-1)*resRatio
rc=(lims['rightCol'])*resRatio
# define subsetted lat and long vector
long = fullHDF5.root.long[lc:rc:1]
lat = fullHDF5.root.lat[tr:br:1]
nrows = len(lat)
ncols = len(long)
#print(ncols)
#print(nrows)
# define subsetted grid
gridsubset = fullHDF5.root.data[tr:br:1,lc:rc:1]
#print(mean(gridsubset))
#print(shape(gridsubset))
#print (type(gridsubset))
##
## now build new hdf5 file for subsetted grid
##
# Initialize hdf5 archive
outHDF5 = tb.openFile(output_path, mode='w', title=gridname+'_'+region)
# build grid metadata
outHDF5.root._v_attrs.asc_file = 'subsetted from hdf5 file: '+input_path
outHDF5.root._v_attrs.ncols = ncols
outHDF5.root._v_attrs.nrows = nrows
outHDF5.root._v_attrs.missing = fullHDF5.root._v_attrs.missing
outHDF5.root._v_attrs.order = fullHDF5.root._v_attrs.order
outHDF5.root._v_attrs.cellsize = fullHDF5.root._v_attrs.cellsize
outHDF5.root._v_attrs.minx = long.min()
outHDF5.root._v_attrs.maxx = long.max()
outHDF5.root._v_attrs.miny = lat.min()
outHDF5.root._v_attrs.maxy = lat.max()
# Add longitude and latitude to archive, uncompressed.
outHDF5.createArray('/','long',long)
outHDF5.createArray('/','lat',lat)
# Add data to archive, heavily compressed, in a chunk array row-by-row (if a row won't fit in memory, the whole array won't fit on disk).
outHDF5.createCArray('/', 'data', tb.Float64Atom(), (nrows, ncols), filters = tb.Filters(complevel=9, complib='zlib'),chunkshape = (1,ncols))
for i in xrange(nrows):
outHDF5.root.data[i,:] = gridsubset[i,:]
# close the files
fullHDF5.close()
outHDF5.close()
################################################################################################################################################
#subsetGlobalGrids(region="AM",lims = master_grid.AM_lims,gridname = "gr071km_y-x+",resRatio=5)
#subsetGlobalGrids(region="AM",lims = master_grid.AM_lims,gridname = "un_mask1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="AM",lims = master_grid.AM_lims,gridname = "salblim1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="AM",lims = master_grid.AM_lims,gridname = "lims1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="AM",lims = master_grid.AM_lims,gridname = "ur1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="AM",lims = master_grid.AM_lims,gridname = "st_mask1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="AM",lims = master_grid.AM_lims,gridname = "salb1km-e2_y-x+",resRatio=5)
#subsetGlobalGrids(region="AF",lims = master_grid.AF_lims,gridname = "gr071km_y-x+",resRatio=5)
#subsetGlobalGrids(region="AF",lims = master_grid.AF_lims,gridname = "un_mask1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="AF",lims = master_grid.AF_lims,gridname = "salblim1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="AF",lims = master_grid.AF_lims,gridname = "lims1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="AF",lims = master_grid.AF_lims,gridname = "ur1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="AF",lims = master_grid.AF_lims,gridname = "st_mask1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="AF",lims = master_grid.AF_lims,gridname = "salb1km-e2_y-x+",resRatio=5)
#subsetGlobalGrids(region="AS",lims = master_grid.AS_lims,gridname = "gr071km_y-x+",resRatio=5)
#subsetGlobalGrids(region="AS",lims = master_grid.AS_lims,gridname = "un_mask1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="AS",lims = master_grid.AS_lims,gridname = "salblim1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="AS",lims = master_grid.AS_lims,gridname = "lims1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="AS",lims = master_grid.AS_lims,gridname = "ur1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="AS",lims = master_grid.AS_lims,gridname = "st_mask1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="AS",lims = master_grid.AS_lims,gridname = "salb1km-e2_y-x+",resRatio=5)
#subsetGlobalGrids(region="KE",lims = master_grid.KE_lims,gridname = "gr071km_y-x+",resRatio=5)
#subsetGlobalGrids(region="KE",lims = master_grid.KE_lims,gridname = "un_mask1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="KE",lims = master_grid.KE_lims,gridname = "salblim1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="KE",lims = master_grid.KE_lims,gridname = "lims1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="KE",lims = master_grid.KE_lims,gridname = "ur1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="KE",lims = master_grid.KE_lims,gridname = "st_mask1km-e_y-x+",resRatio=5)
#subsetGlobalGrids(region="KE",lims = master_grid.KE_lims,gridname = "salb1km-e2_y-x+",resRatio=5)
subsetGlobalGrids(region="KE",lims = master_grid.KE_lims,gridname = "st_mask5km-e_y-x+",resRatio=1)
subsetGlobalGrids(region="KE",lims = master_grid.KE_lims,gridname = "gr075km_y-x+",resRatio=1)