tsg_create_3d.py - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tsg_create_3d.py (7323B)
---
1 #!/usr/bin/env python3
2 #
3 # Copyright (C) 2011, 2014, 2018, 2019 Andy Aschwanden
4 #
5 # This file is part of PISM.
6 #
7 # PISM is free software; you can redistribute it and/or modify it under the
8 # terms of the GNU General Public License as published by the Free Software
9 # Foundation; either version 3 of the License, or (at your option) any later
10 # version.
11 #
12 # PISM is distributed in the hope that it will be useful, but WITHOUT ANY
13 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 # FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
15 # details.
16 #
17 # You should have received a copy of the GNU General Public License
18 # along with PISM; if not, write to the Free Software
19 # Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20
21
22 import sys
23 import time
24 import numpy as np
25 from pyproj import Proj
26 from sys import stderr
27
28 write = stderr.write
29
30 from netCDF4 import Dataset as CDF
31
32 from optparse import OptionParser
33
34 ice_density = 910.0
35
36 __author__ = "Andy Aschwanden"
37
38 # Create PISM-readable input file from Storglaciaren DEM
39
40 parser = OptionParser()
41 parser.usage = "usage: %prog [options]"
42 parser.description = "Preprocess Storglaciaren files."
43
44
45 (options, args) = parser.parse_args()
46
47
48 # Create PISM-readable input file from Storglaciaren DEM
49
50
51 write('------------------------------\n')
52 write('PISM-Storglaciaren example\n')
53 write('------------------------------\n')
54
55 # data dir
56 data_dir = './'
57 # X is Northing (http://en.wikipedia.org/wiki/Swedish_grid)
58 XFile = data_dir + 'X.txt.gz'
59 # Y is Easting (http://en.wikipedia.org/wiki/Swedish_grid)
60 YFile = data_dir + 'Y.txt.gz'
61 # Bed and Surface DEMs for Storglaciaren
62 zBaseFile = data_dir + 'zBase.txt.gz'
63 zSurfFile = data_dir + 'zSurf.txt.gz'
64
65 # load coordinate information. Note: Swedish grid (RT90) uses inverse notation
66 # X -> northing, Y -> easting
67 try:
68 write('Reading northing coordinate infos from %s: ' % XFile)
69 X = np.loadtxt(XFile)
70 write('Done.\n')
71 write('Reading easting coordinate infos from %s: ' % YFile)
72 Y = np.loadtxt(YFile)
73 write('Done.\n')
74 except IOError:
75 write('ERROR: File %s or %s could not be found.\n' % XFile % YFile)
76 exit(2)
77
78 # load Bed DEM
79 try:
80 write('Reading DEM from %s: ' % zBaseFile)
81 zBase = np.loadtxt(zBaseFile)
82 write('Done.\n')
83 except IOError:
84 write('ERROR: File %s could not be found.\n' % zBaseFile)
85 exit(2)
86
87 # load Surface DEM
88 try:
89 write('Reading DEM from %s: ' % zSurfFile)
90 zSurf = np.loadtxt(zSurfFile)
91 write('Done.\n')
92 except IOError:
93 write('ERROR: File %s could not be found.\n' % zSurfFile)
94 exit(2)
95
96 # Grid size. DEM has 10m spacing.
97 N = zBase.shape[1] - 1
98 M = zBase.shape[0] - 1
99 Ne = N + 30
100 Me = M
101 n0 = X[-1, 0]
102 e0 = Y[0, 0]
103 de = 10 # m
104 dn = 10 # m
105
106 e1 = e0 + (Ne - 1) * de
107 n1 = n0 + (Me - 1) * dn
108
109 easting = np.linspace(e0, e1, Ne)
110 northing = np.linspace(n0, n1, Me)
111
112 # convert to lat/lon
113
114 ee, nn = np.meshgrid(easting, northing)
115
116 projection = Proj(init="epsg:3021")
117 longitude, latitude = projection(ee, nn, inverse=True)
118
119 write("Coordinates of the lower-left grid corner:\n"
120 " easting = %.0f\n"
121 " northing = %.0f\n"
122 "Grid size:\n"
123 " rows = %d\n"
124 " columns = %d\n" % (e0, n0, N, M))
125
126 # Fill value
127 fill_value = -9999
128 bed_valid_min = -5000.0
129 thk_valid_min = 0.0
130
131 bed = np.flipud(zBase)
132 dem = np.flipud(zSurf) # ignored by bootstrapping
133 thk = np.flipud(zSurf-zBase) # used for bootstrapping
134
135 m_bed = np.zeros((Me, Ne)) + 1100
136 m_bed[0:M, 0:N] = bed[0:M, 0:N]
137
138 m_dem = np.zeros((Me, Ne))
139 m_dem[0:M, 0:N] = dem[0:M, 0:N]
140
141 m_thk = np.zeros((Me, Ne))
142 m_thk[0:M, 0:N] = thk[0:M, 0:N]
143
144 # Replace NaNs with zeros
145 thk = np.nan_to_num(thk)
146 m_thk = np.nan_to_num(m_thk)
147 dem = np.nan_to_num(dem)
148 m_dem = np.nan_to_num(m_dem)
149
150 # There are some negative thickness values
151 # Quick and dirty: set to zero
152 # some inconsistencies in the original data still needs to be sorted out
153 # (filtering)
154 thk[thk < 0] = 0
155 m_thk[m_thk < 0] = 0
156
157
158 # Output filename
159 ncfile = 'pism_storglaciaren_3d.nc'
160
161 # Write the data:
162 nc = CDF(ncfile, "w", format='NETCDF3_CLASSIC') # for netCDF4 module
163
164 # Create dimensions x and y
165 nc.createDimension("x", size=easting.shape[0])
166 nc.createDimension("y", size=northing.shape[0])
167
168 x = nc.createVariable("x", 'f4', dimensions=("x",))
169 x.units = "m"
170 x.long_name = "easting"
171 x.standard_name = "projection_x_coordinate"
172
173 y = nc.createVariable("y", 'f4', dimensions=("y",))
174 y.units = "m"
175 y.long_name = "northing"
176 y.standard_name = "projection_y_coordinate"
177
178 x[:] = easting
179 y[:] = northing
180
181 from scipy.interpolate import griddata
182
183
184 def def_var(nc, name, units, fillvalue):
185 var = nc.createVariable(name, 'f', dimensions=("y", "x"), fill_value=fillvalue)
186 var.units = units
187 return var
188
189
190 lon_var = def_var(nc, "lon", "degrees_east", None)
191 lon_var.standard_name = "longitude"
192 lon_var[:] = longitude
193
194 lat_var = def_var(nc, "lat", "degrees_north", None)
195 lat_var.standard_name = "latitude"
196 lat_var[:] = latitude
197
198 bed_var = def_var(nc, "topg", "m", fill_value)
199 bed_var.valid_min = bed_valid_min
200 bed_var.standard_name = "bedrock_altitude"
201 bed_var.coordinates = "lat lon"
202 bed_var[:] = m_bed
203
204 thk_var = def_var(nc, "thk", "m", fill_value)
205 thk_var.valid_min = thk_valid_min
206 thk_var.standard_name = "land_ice_thickness"
207 thk_var.coordinates = "lat lon"
208 thk_var[:] = m_thk
209
210 dem_var = def_var(nc, "usurf_from_dem", "m", fill_value)
211 dem_var.standard_name = "surface_altitude"
212 dem_var.coordinates = "lat lon"
213 dem_var[:] = m_dem
214
215 ftt_mask = def_var(nc, "ftt_mask", "", fill_value)
216 ftt_mask[:] = 1
217
218
219 # generate (somewhat) reasonable acab
220 acab_max = 2.5 # m/a
221 acab_min = -3.0 # m/a
222 acab_up = easting.min() + 200 # m; location of upstream end of linear acab
223 acab_down = easting.max() - 900 # m;location of downstream end of linear acab
224 acab = np.ones_like(longitude)
225
226 acab[:] = acab_min
227 acab[:] = acab_max - (acab_max-acab_min) * (easting - acab_up) / (acab_down - acab_up)
228 acab[m_thk < 1] = acab_min
229
230 acab_var = def_var(nc, "climatic_mass_balance", "kg m-2 year-1", fill_value)
231 acab_var.standard_name = "land_ice_surface_specific_mass_balance"
232 acab_var[:] = acab * ice_density
233
234 # Set boundary conditions for Scandinavian-type polythermal glacier
235 # ------------------------------------------------------------------------------
236 #
237 # (A) Surface temperature for temperature equation bc
238 T0 = 273.15 # K
239 Tma = -6.0 # degC, mean annual air temperature at Tarfala
240 zcts = 1300 # m a.s.l.; altitude where CTS is at the surface, projected to topg
241 slope = 100 # m; range around which surface temp transition happens
242
243
244 # smoothed version; FIXME: can't we at least have it depend on initial DEM?
245 # additional lapse rate?
246 artm = T0 + Tma * (zcts + slope - m_bed) / (2.0 * slope)
247 artm[m_bed < zcts-slope] = T0 + Tma
248 artm[m_bed > zcts+slope] = T0
249
250 artm_var = def_var(nc, "ice_surface_temp", "K", fill_value)
251 artm_var[:] = artm
252
253 # set global attributes
254 nc.Conventions = "CF-1.4"
255 historysep = ' '
256 historystr = time.asctime() + ': ' + historysep.join(sys.argv) + '\n'
257 setattr(nc, 'history', historystr)
258 # PROJ string equivalent to EPSG 3021
259 nc.proj = "+proj=tmerc +lat_0=0 +lon_0=15.80827777777778 +k=1 +x_0=1500000 +y_0=0 +ellps=bessel +towgs84=419.384,99.3335,591.345,0.850389,1.81728,-7.86224,-0.99496 +units=m +no_defs"
260 nc.close()
261 write('Done writing NetCDF file %s!\n' % ncfile)