treadgeom.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
       ---
       treadgeom.py (2611B)
       ---
            1 #!/usr/bin/env python3
            2 
            3 # readgeom   reads bed elevation and ice thickness from BEDMAP2 files
            4 #            and either shows as figures or writes into NetCDF file
            5 
            6 # edit these:
            7 BM2PATH = '/home/ed/Desktop/bedmap2_bin/'
            8 showgeom = False
            9 
           10 import numpy as np
           11 import numpy.ma as ma
           12 
           13 import matplotlib.pyplot as plt
           14 import sys
           15 
           16 try:
           17     from netCDF4 import Dataset as NC
           18 except:
           19     print("netCDF4 is not installed!")
           20     sys.exit(1)
           21 
           22 from PISMNC import PISMDataset as PNC
           23 
           24 
           25 N = 6667
           26 
           27 print("reading bedmap2 binary files from %s ..." % (BM2PATH))
           28 
           29 fname = BM2PATH + 'bedmap2_bed.flt'
           30 bed = ma.masked_equal(np.reshape(np.fromfile(fname, dtype=np.float32), (N, N)), -9999.0)
           31 fname = BM2PATH + 'bedmap2_thickness.flt'
           32 thk = ma.masked_equal(np.reshape(np.fromfile(fname, dtype=np.float32), (N, N)), -9999.0)
           33 fname = BM2PATH + 'bedmap2_icemask_grounded_and_shelves.flt'
           34 msk = ma.masked_equal(np.reshape(np.fromfile(fname, dtype=np.float32), (N, N)), -9999.0)
           35 
           36 print("  range of bed = [%.2f, %.2f]" % (bed.min(), bed.max()))
           37 print("  range of thk = [%.2f, %.2f]" % (thk.min(), thk.max()))
           38 print("  range of msk = [%.2f, %.2f]" % (msk.min(), msk.max()))
           39 
           40 if showgeom:
           41     print("showing fields ...")
           42 
           43     fig = plt.figure(1)
           44     ax = plt.imshow(msk)
           45     fig.colorbar(ax)
           46 
           47     fig = plt.figure(2)
           48     ax = plt.imshow(bed)
           49     fig.colorbar(ax)
           50 
           51     fig = plt.figure(3)
           52     ax = plt.imshow(thk)
           53     fig.colorbar(ax)
           54 
           55     plt.show()
           56 
           57     sys.exit(0)
           58 
           59 #err = abs(bed[msk<0.5] + thk[msk<0.5] - srf[msk<0.5])
           60 # print err.max()
           61 
           62 outname = 'ant1kmgeom.nc'
           63 
           64 print("writing NetCDF file '%s' ..." % outname)
           65 try:
           66     nc = PNC(outname, 'w', format='NETCDF3_CLASSIC')
           67 except:
           68     print("can't open file %s for writing" % outname)
           69     exit(1)
           70 
           71 print("  writing x,y ...")
           72 dx = 1000.0
           73 dy = 1000.0
           74 x = np.linspace(0.0, (N - 1) * dx, N)
           75 y = np.linspace(0.0, (N - 1) * dy, N)
           76 nc.create_dimensions(x, y, time_dependent=False)
           77 
           78 print("  writing topg ...")
           79 nc.define_2d_field("topg", time_dependent=False,
           80                    attrs={"long_name": "elevation of bedrock",
           81                           "valid_range": (-9000.0, 9000.0),
           82                           "standard_name": "bedrock_altitude",
           83                           "units": "meters"})
           84 nc.write_2d_field("topg", bed)
           85 
           86 print("  writing thk ...")
           87 nc.define_2d_field("thk", time_dependent=False,
           88                    attrs={"long_name": "thickness of ice sheet or ice shelf",
           89                           "valid_range": (0.0, 9000.0),
           90                           "standard_name": "land_ice_thickness",
           91                           "units": "meters"})
           92 nc.write_2d_field("thk", thk)
           93 
           94 nc.close()
           95 print("done")