Skip to content

Commit 7e7cc6f

Browse files
lee1043lee1043
committed
python 2 style to python 3 style, pre-commit fix
1 parent 5c7b40c commit 7e7cc6f

File tree

2 files changed

+1281
-1055
lines changed

2 files changed

+1281
-1055
lines changed

‎pcmdi_metrics/sea_ice/generate_sector_masks.py

Lines changed: 193 additions & 170 deletions
Original file line numberDiff line numberDiff line change
@@ -1,195 +1,218 @@
1+
import os
2+
import string
3+
import sys
4+
15
import cdms2
26
import MV2 as MV
3-
import string
4-
import os,sys
5-
from pcmdi_metrics.pcmdi.pmp_parser import *
6-
import pcmdi_metrics
7-
import parser
8-
import collections
97
from sector_mask_defs import *
108

11-
import argparse
12-
from argparse import RawTextHelpFormatter
9+
from pcmdi_metrics.pcmdi.pmp_parser import *
1310

1411
P = PMPParser()
1512

16-
P.add_argument("--mp", "--modpath",
17-
type = str,
18-
dest = 'modpath',
19-
default = '',
20-
help = "Explicit path to model monthly PR climatology")
21-
P.add_argument("-o", "--obspath",
22-
type = str,
23-
dest = 'obspath',
24-
default = '',
25-
help = "Explicit path to obs monthly PR climatology")
26-
P.add_argument("--outpd", "--outpathdata",
27-
type = str,
28-
dest = 'outpathdata',
29-
default = '/export/gleckler1/processing/metrics_package/my_test/sea_ice/git_data/',
30-
help = "Output path for sector scale masks")
13+
P.add_argument(
14+
"--mp",
15+
"--modpath",
16+
type=str,
17+
dest="modpath",
18+
default="",
19+
help="Explicit path to model monthly PR climatology",
20+
)
21+
P.add_argument(
22+
"-o",
23+
"--obspath",
24+
type=str,
25+
dest="obspath",
26+
default="",
27+
help="Explicit path to obs monthly PR climatology",
28+
)
29+
P.add_argument(
30+
"--outpd",
31+
"--outpathdata",
32+
type=str,
33+
dest="outpathdata",
34+
default="/export/gleckler1/processing/metrics_package/my_test/sea_ice/git_data/",
35+
help="Output path for sector scale masks",
36+
)
3137

3238
args = P.parse_args(sys.argv[1:])
3339
sec_mask_dir = args.outpathdata
3440

35-
#Factors
36-
factor1=1.e-6 #model units are m^2, converting to km^2
37-
factor2=1.e-2 #model units are %, converting to non-dimen
38-
a=6371.009 #Earth radii in km
39-
pi=22./7.
40-
factor3=4.*pi*a*a # Earth's surface area in km2
41-
dc=0.15 #minimum ice concentration contour
41+
# Factors
42+
factor1 = 1.0e-6 # model units are m^2, converting to km^2
43+
factor2 = 1.0e-2 # model units are %, converting to non-dimen
44+
a = 6371.009 # Earth radii in km
45+
pi = 22.0 / 7.0
46+
factor3 = 4.0 * pi * a * a # Earth's surface area in km2
47+
dc = 0.15 # minimum ice concentration contour
4248

43-
pin = '/work/gleckler1/processed_data/cmip5clims-historical/sic/cmip5.MOD.historical.r1i1p1.mo.seaIce.OImon.sic.ver-1.1980-2005.SC.nc'
49+
pin = "/work/gleckler1/processed_data/cmip5clims-historical/sic/cmip5.MOD.historical.r1i1p1.mo.seaIce.OImon.sic.ver-1.1980-2005.SC.nc"
4450

45-
pins = string.replace(pin,'MOD','*')
51+
pins = string.replace(pin, "MOD", "*")
4652

47-
lst = os.popen('ls ' + pins).readlines()
53+
lst = os.popen("ls " + pins).readlines()
4854

4955
mods = []
5056
mods_failed = []
51-
for l in lst:
52-
mod = string.split(l,'.')[1]
53-
if mod not in mods: mods.append(mod)
57+
for li in lst:
58+
mod = string.split(li, ".")[1]
59+
if mod not in mods:
60+
mods.append(mod)
5461

55-
#w =sys.stdin.readline()
62+
# w =sys.stdin.readline()
5663

57-
var = 'sic'
64+
var = "sic"
5865
factor2 = 1
5966

60-
mods = ['ACCESS1-3']
67+
mods = ["ACCESS1-3"]
6168

6269
for mod in mods:
63-
try:
64-
fc = string.replace(pin,'MOD',mod)
65-
f=cdms2.open(fc)
66-
sic=f(var, squeeze=1)
67-
sic_grid=sic.getGrid()
68-
lat=sic.getLatitude()
69-
lon=sic.getLongitude()
70-
sic=MV.multiply(sic,factor2)
71-
72-
print 'CMIP5-native= ',MV.max(sic)
73-
if MV.rank(lat)==1:
74-
tmp2d=f(var,time=slice(0,1),squeeze=1)
75-
lats=MV.zeros(tmp2d.shape)
76-
for ii in range (0,len(lon)):
77-
lats[:,ii]=lat[:]
78-
else:
79-
lats=lat
80-
81-
if MV.rank(lon)==1:
82-
tmp2d=f(var,time=slice(0,1),squeeze=1)
83-
lons=MV.zeros(tmp2d.shape)
84-
for ii in range (0,len(lat)):
85-
lons[ii,:]=lon[:]
86-
else:
87-
lons=lon
88-
89-
f.close()
90-
91-
92-
#######################################################
93-
### areacello
94-
area_dir = '/work/cmip5/fx/fx/areacello/'
95-
alist = os.listdir(area_dir) # LIST OF ALL AREACELLO FILES
96-
97-
for a in alist:
98-
if string.find(a,mod) != -1:
99-
areaf = a
100-
print mod,' ', a
101-
# w = sys.stdin.readline()
102-
103-
g = cdms2.open(area_dir + areaf)
104-
105-
try:
106-
area = g('areacello')
107-
except:
108-
area = g('areacella')
109-
area = MV.multiply(area,factor1)
110-
area = MV.multiply(area,factor1)
111-
112-
g.close()
113-
114-
land_mask=MV.zeros(area.shape)
115-
116-
# Reading the ocean/land grid cell masks (sftof/sftlf)
117-
mask_dir = '/work/cmip5/fx/fx/sftof/'
118-
mlist = os.listdir(mask_dir) # LIST OF ALL AREACELLO FILES
119-
120-
for m in mlist:
121-
if string.find(m,mod) != -1:
122-
maskf = m
123-
print mod,' ', m
124-
125-
s = cdms2.open(mask_dir + maskf)
126-
try:
127-
frac = s('sftof')
128-
if (mod != 'MIROC5' and mod!= 'GFDL-CM2p1' and mod!='GFDL-CM3' and mod!='GFDL-ESM2M'):
129-
frac = MV.multiply(frac,factor2)
130-
print mod,MV.max(frac)
131-
area = MV.multiply(area,frac)
132-
land_mask = MV.multiply(1,(1-frac))
133-
except:
134-
frac = s('sftlf')
135-
if (mod != 'MIROC5' or mod!= 'GFDL-CM2p1' or mod!='GFDL-CM3' or mod!='GFDL-ESM2M'):
136-
frac = MV.multiply(frac,factor2)
137-
area = MV.multiply(area,(1-frac))
138-
land_mask = MV.multiply(1,frac)
139-
s.close()
140-
141-
# Creating regional masks
142-
# GFDL and bcc model grids have shift of 80
143-
if (mod[0:4] == 'GFDL' or mod[0:3] == 'bcc') :
144-
lons_a=MV.where(MV.less(lons,-180.),lons+360.,lons)
145-
lons_p=MV.where(MV.less(lons,0.),lons+360.,lons)
146-
else:
147-
lons_a=MV.where(MV.greater(lons,180.),lons-360.,lons)
148-
lons_p=lons
149-
150-
print 'CMIP5'
151-
print 'lons_na= ',MV.min(lons_a),MV.max(lons_a)
152-
print 'lons_np= ',MV.min(lons_p),MV.max(lons_p)
153-
# mask_arctic=MV.zeros(area.shape)
154-
# mask_antarctic=MV.zeros(area.shape)
155-
mask_ca=MV.zeros(area.shape)
156-
mask_na=MV.zeros(area.shape)
157-
mask_np=MV.zeros(area.shape)
158-
mask_sa=MV.zeros(area.shape)
159-
mask_sp=MV.zeros(area.shape)
160-
mask_io=MV.zeros(area.shape)
161-
162-
###############################################################
163-
164-
sectors = ['ca','na','np','sp','sa','io']
165-
166-
for sector in sectors:
167-
mask = getmask(sector,lats,lons,lons_a,lons_p,land_mask)
168-
169-
g = cdms2.open(sec_mask_dir + 'mask_' + mod + '_' + sector + '.nc','w+')
170-
mask.id = 'mask'
171-
g.write(mask)
172-
# land_mask.id = 'sftof'
173-
# g.write(land_mask)
174-
g.close()
175-
176-
print 'got it!',' ', mask.shape
177-
w = sys.stdin.readline()
178-
179-
# Calculate the Total Sea Ice Concentration/Extent/Area
180-
# ice_area = MV.multiply(sic,1) #SIC
181-
area = MV.multiply(1,area) #SIE
182-
ice_area = MV.multiply(sic,area) #SIA
183-
ice_area = MV.where(MV.greater_equal(sic,0.15),ice_area,0.) #Masking out the sic<0.15
184-
185-
# arctic=MV.logical_and(MV.greater_equal(lats,35.),MV.less(lats,87.2)) #SSM/I limited to 87.2N
186-
mask_arctic=MV.logical_and(MV.greater_equal(lats,45.),MV.less(lats,90.)) #Adding currently in SSM/I 100% in the area >87.2N
187-
mask_antarctic=MV.logical_and(MV.greater_equal(lats,-90.),MV.less(lats,-55.))
188-
189-
except:
190-
'Failed for ', mod
191-
mods_failed.append(mod)
192-
193-
print 'failed for ', mods_failed
70+
try:
71+
fc = string.replace(pin, "MOD", mod)
72+
f = cdms2.open(fc)
73+
sic = f(var, squeeze=1)
74+
sic_grid = sic.getGrid()
75+
lat = sic.getLatitude()
76+
lon = sic.getLongitude()
77+
sic = MV.multiply(sic, factor2)
78+
79+
print("CMIP5-native= ", MV.max(sic))
80+
if MV.rank(lat) == 1:
81+
tmp2d = f(var, time=slice(0, 1), squeeze=1)
82+
lats = MV.zeros(tmp2d.shape)
83+
for ii in range(0, len(lon)):
84+
lats[:, ii] = lat[:]
85+
else:
86+
lats = lat
87+
88+
if MV.rank(lon) == 1:
89+
tmp2d = f(var, time=slice(0, 1), squeeze=1)
90+
lons = MV.zeros(tmp2d.shape)
91+
for ii in range(0, len(lat)):
92+
lons[ii, :] = lon[:]
93+
else:
94+
lons = lon
95+
96+
f.close()
97+
98+
#######################################################
99+
### areacello
100+
area_dir = "/work/cmip5/fx/fx/areacello/"
101+
alist = os.listdir(area_dir) # LIST OF ALL AREACELLO FILES
102+
103+
for a in alist:
104+
if string.find(a, mod) != -1:
105+
areaf = a
106+
print(mod, " ", a)
107+
108+
# w = sys.stdin.readline()
109+
110+
g = cdms2.open(area_dir + areaf)
111+
112+
try:
113+
area = g("areacello")
114+
except:
115+
area = g("areacella")
116+
area = MV.multiply(area, factor1)
117+
area = MV.multiply(area, factor1)
118+
119+
g.close()
120+
121+
land_mask = MV.zeros(area.shape)
122+
123+
# Reading the ocean/land grid cell masks (sftof/sftlf)
124+
mask_dir = "/work/cmip5/fx/fx/sftof/"
125+
mlist = os.listdir(mask_dir) # LIST OF ALL AREACELLO FILES
126+
127+
for m in mlist:
128+
if string.find(m, mod) != -1:
129+
maskf = m
130+
print(mod, " ", m)
131+
132+
s = cdms2.open(mask_dir + maskf)
133+
try:
134+
frac = s("sftof")
135+
if (
136+
mod != "MIROC5"
137+
and mod != "GFDL-CM2p1"
138+
and mod != "GFDL-CM3"
139+
and mod != "GFDL-ESM2M"
140+
):
141+
frac = MV.multiply(frac, factor2)
142+
print(mod, MV.max(frac))
143+
area = MV.multiply(area, frac)
144+
land_mask = MV.multiply(1, (1 - frac))
145+
except:
146+
frac = s("sftlf")
147+
if (
148+
mod != "MIROC5"
149+
or mod != "GFDL-CM2p1"
150+
or mod != "GFDL-CM3"
151+
or mod != "GFDL-ESM2M"
152+
):
153+
frac = MV.multiply(frac, factor2)
154+
area = MV.multiply(area, (1 - frac))
155+
land_mask = MV.multiply(1, frac)
156+
s.close()
157+
158+
# Creating regional masks
159+
# GFDL and bcc model grids have shift of 80
160+
if mod[0:4] == "GFDL" or mod[0:3] == "bcc":
161+
lons_a = MV.where(MV.less(lons, -180.0), lons + 360.0, lons)
162+
lons_p = MV.where(MV.less(lons, 0.0), lons + 360.0, lons)
163+
else:
164+
lons_a = MV.where(MV.greater(lons, 180.0), lons - 360.0, lons)
165+
lons_p = lons
166+
167+
print("CMIP5")
168+
print("lons_na= ", MV.min(lons_a), MV.max(lons_a))
169+
print("lons_np= ", MV.min(lons_p), MV.max(lons_p))
170+
# mask_arctic=MV.zeros(area.shape)
171+
# mask_antarctic=MV.zeros(area.shape)
172+
mask_ca = MV.zeros(area.shape)
173+
mask_na = MV.zeros(area.shape)
174+
mask_np = MV.zeros(area.shape)
175+
mask_sa = MV.zeros(area.shape)
176+
mask_sp = MV.zeros(area.shape)
177+
mask_io = MV.zeros(area.shape)
178+
179+
###############################################################
180+
181+
sectors = ["ca", "na", "np", "sp", "sa", "io"]
182+
183+
for sector in sectors:
184+
mask = getmask(sector, lats, lons, lons_a, lons_p, land_mask)
185+
186+
g = cdms2.open(sec_mask_dir + "mask_" + mod + "_" + sector + ".nc", "w+")
187+
mask.id = "mask"
188+
g.write(mask)
189+
# land_mask.id = 'sftof'
190+
# g.write(land_mask)
191+
g.close()
192+
193+
print("got it!", " ", mask.shape)
194+
w = sys.stdin.readline()
195+
196+
# Calculate the Total Sea Ice Concentration/Extent/Area
197+
# ice_area = MV.multiply(sic,1) #SIC
198+
area = MV.multiply(1, area) # SIE
199+
ice_area = MV.multiply(sic, area) # SIA
200+
ice_area = MV.where(
201+
MV.greater_equal(sic, 0.15), ice_area, 0.0
202+
) # Masking out the sic<0.15
203+
204+
# arctic=MV.logical_and(MV.greater_equal(lats,35.),MV.less(lats,87.2)) #SSM/I limited to 87.2N
205+
mask_arctic = MV.logical_and(
206+
MV.greater_equal(lats, 45.0), MV.less(lats, 90.0)
207+
) # Adding currently in SSM/I 100% in the area >87.2N
208+
mask_antarctic = MV.logical_and(
209+
MV.greater_equal(lats, -90.0), MV.less(lats, -55.0)
210+
)
211+
212+
except:
213+
"Failed for ", mod
214+
mods_failed.append(mod)
215+
216+
print("failed for ", mods_failed)
194217

195218
# Calculate the Sea Ice Covered Area

0 commit comments

Comments
 (0)