pism_python.py (2829B)
1 #!/usr/bin/env python3 2 3 # Copyright (C) 2009-2015, 2018 the PISM Authors 4 5 # @package pism_python 6 # \author the PISM authors 7 # \brief Creates "from scratch" a boring dataset with the right format 8 # to use as a PISM bootstrapping file. 9 # \details Example use of Python for this purpose. 10 # 11 # Usage, including a minimal PISM call to bootstrap from this file: 12 # 13 # \verbatim $ pism_python.py # creates foo.nc \endverbatim 14 # \verbatim $ pismr -i foo.nc -bootstrap -Mx 41 -My 41 -Mz 21 -Lz 4000 -Mbz 5 -Lbz 500 -y 1 \endverbatim 15 16 import sys 17 import time 18 import numpy as np 19 20 # try different netCDF modules 21 try: 22 from netCDF4 import Dataset as CDF 23 except: 24 print("netCDF4 is not installed!") 25 sys.exit(1) 26 27 # set up the grid: 28 Lx = 1e6 29 Ly = 1e6 30 Mx = 51 31 My = 71 32 x = np.linspace(-Lx, Lx, Mx) 33 y = np.linspace(-Ly, Ly, My) 34 35 # create dummy fields 36 [xx, yy] = np.meshgrid(x, y) # if there were "ndgrid" in numpy we would use it 37 acab = np.zeros((Mx, My)) 38 artm = np.zeros((Mx, My)) + 273.15 + 10.0 # 10 degrees Celsius 39 topg = 1000.0 + 200.0 * (xx + yy) / max(Lx, Ly) # change "1000.0" to "0.0" to test 40 # flotation criterion, etc. 41 thk = 3000.0 * (1.0 - 3.0 * (xx ** 2 + yy ** 2) / Lx ** 2) 42 thk[thk < 0.0] = 0.0 43 44 # Output filename 45 ncfile = 'foo.nc' 46 47 # Write the data: 48 nc = CDF(ncfile, "w", format='NETCDF3_CLASSIC') # for netCDF4 module 49 50 # Create dimensions x and y 51 nc.createDimension("x", size=Mx) 52 nc.createDimension("y", size=My) 53 54 x_var = nc.createVariable("x", 'f4', dimensions=("x",)) 55 x_var.units = "m" 56 x_var.long_name = "easting" 57 x_var.standard_name = "projection_x_coordinate" 58 x_var[:] = x 59 60 y_var = nc.createVariable("y", 'f4', dimensions=("y",)) 61 y_var.units = "m" 62 y_var.long_name = "northing" 63 y_var.standard_name = "projection_y_coordinate" 64 y_var[:] = y 65 66 fill_value = np.nan 67 68 69 def def_var(nc, name, units, fillvalue): 70 # dimension transpose is standard: "float thk(y, x)" in NetCDF file 71 var = nc.createVariable(name, 'f', dimensions=("y", "x"), fill_value=fillvalue) 72 var.units = units 73 return var 74 75 76 bed_var = def_var(nc, "topg", "m", fill_value) 77 bed_var.standard_name = "bedrock_altitude" 78 bed_var[:] = topg 79 80 thk_var = def_var(nc, "thk", "m", fill_value) 81 thk_var.standard_name = "land_ice_thickness" 82 thk_var[:] = thk 83 84 acab_var = def_var(nc, "climatic_mass_balance", "m year-1", fill_value) 85 acab_var.standard_name = "land_ice_surface_specific_mass_balance" 86 acab_var[:] = acab 87 88 artm_var = def_var(nc, "ice_surface_temp", "K", fill_value) 89 artm_var[:] = artm 90 91 # set global attributes 92 nc.Conventions = "CF-1.4" 93 historysep = ' ' 94 historystr = time.asctime() + ': ' + historysep.join(sys.argv) + '\n' 95 setattr(nc, 'history', historystr) 96 97 nc.close() 98 print(' PISM-bootable NetCDF file %s written' % ncfile) 99 print(' for example, run:') 100 print(' $ pismr -i foo.nc -bootstrap -Mx 41 -My 41 -Mz 21 -Lz 4000 -Mbz 5 -Lbz 500 -y 1')