Python Dataset.projection方法代码示例

# 需要导入模块: from netCDF4 import Dataset [as 别名]
# 或者: from netCDF4.Dataset import projection [as 别名]
def fix_netcdf(infile,outfile):
    """
    Write a new netcdf but this time do the coordinate vars correctly
    """
    rootgrp = Dataset(outfile,'w', format='NETCDF3_64BIT')

    data, targetAttrs = read_netcdf(infile,vars=('Prec','Wind','Tmax','Tmin','time','nav_lat','nav_lon'))
    res = 0.5
    # set dimensions
    lon = rootgrp.createDimension('lon',data['Prec'].shape[2])
    lat = rootgrp.createDimension('lat',data['Prec'].shape[1])
    time = rootgrp.createDimension('time',data['Prec'].shape[0])

    # do vars
    
    times = rootgrp.createVariable('time','f8',('time',))
    times[:] = np.arange(data['Prec'].shape[0])*86400
    times.units = targetAttrs['time']['units']
    times.long_name = targetAttrs['time']['long_name']
    
    lat = rootgrp.createVariable('lat','f8',('lat',))
    lat[:] = np.arange(data['nav_lat'].min(),data['nav_lat'].max()+res,res)
    lat.units = 'degrees_north'
    lat.long_name = 'Latitude'
    
    lon = rootgrp.createVariable('lon','f8',('lon',))
    lon[:] = np.arange(data['nav_lon'].min(),data['nav_lon'].max()+res,res)
    lon.units = 'degrees_east'
    lon.long_name = 'Longitude'
    
    Precip = rootgrp.createVariable('Precip','f8',('time','lat','lon',),fill_value=data['Prec'].fill_value)
    Precip[:,:,:] = data['Prec']
    Precip.units = targetAttrs['Prec']['units']
    Precip.long_name = targetAttrs['Prec']['long_name']
    
    Tmax = rootgrp.createVariable('Tmax','f8',('time','lat','lon',),fill_value=data['Tmax'].fill_value)
    Tmax[:,:,:] = data['Tmax']
    Tmax.units = targetAttrs['Tmax']['units']
    Tmax.long_name = targetAttrs['Tmax']['long_name']

    Tmin = rootgrp.createVariable('Tmin','f8',('time','lat','lon',),fill_value=data['Tmin'].fill_value)
    Tmin[:,:,:] = data['Tmin']
    Tmin.units = targetAttrs['Tmin']['units']
    Tmin.long_name = targetAttrs['Tmin']['long_name']    

    Wind = rootgrp.createVariable('Wind','f8',('time','lat','lon',),fill_value=data['Wind'].fill_value)
    Wind[:,:,:] = data['Wind']
    Wind.units = targetAttrs['Wind']['units']
    Wind.long_name = targetAttrs['Wind']['long_name']
    
    rootgrp.description = 'Global 1/2 Degree Gridded Meteorological VIC Forcing Data Set '
    rootgrp.history = 'Created: {}\n'.format(tm.ctime(tm.time()))
    rootgrp.source = sys.argv[0] # prints the name of script used
    rootgrp.institution = "University of Washington Dept. of Civil and Environmental Engineering"
    rootgrp.sources = "UDel (Willmott and Matsuura 2007), CRU (Mitchell et al., 2004), NCEP/NCAR (Kalnay et al. 1996)"
    rootgrp.projection = "Geographic"
    rootgrp.surfSng_convention = "Traditional"

    rootgrp.close()

# 需要导入模块: from netCDF4 import Dataset [as 别名]
# 或者: from netCDF4.Dataset import projection [as 别名]
def write_netcdf(ncfile, format, domain, data):
    nc = Dataset(ncfile, 'w', format=format)
    nc.createDimension('x', domain['nx'])
    nc.createDimension('y', domain['ny'])

    nav_lon = nc.createVariable('nav_lon', 'f4', ('y','x',))
    nav_lat = nc.createVariable('nav_lat', 'f4', ('y','x',))
    soilsat_inst = nc.createVariable('soilsat_inst', 'f4', ('y','x',))
    soilsat_mean = nc.createVariable('soilsat_mean', 'f4', ('y','x',))
    soilsat_anom = nc.createVariable('soilsat_anom', 'f4', ('y','x',))
    soilsat_perc = nc.createVariable('soilsat_perc', 'f4', ('y','x',))

    nav_lon[:] = data['nav_lon']
    nav_lat[:] = data['nav_lat']
    soilsat_inst[:] = data['soilsat_inst']
    soilsat_mean[:] = data['soilsat_mean']
    soilsat_anom[:] = data['soilsat_anom']
    soilsat_perc[:] = data['soilsat_perc']

    nav_lon.units = 'degrees_east'
    nav_lat.units = 'degrees_north'
    soilsat_inst.units = '-'
    soilsat_mean.units = '-'
    soilsat_anom.units = '-'
    soilsat_perc.units = '-'

    nav_lon.long_name = 'Longitude'
    nav_lat.long_name = 'Latitude'
    soilsat_inst.long_name = 'Soil saturation'
    soilsat_mean.long_name = 'Soil saturation mean'
    soilsat_anom.long_name = 'Soil saturation anomaly'
    soilsat_perc.long_name = 'Soil saturation percentile'

    nav_lon.FillValue = default_fillvals['f4']
    nav_lat.FillValue = default_fillvals['f4']
    soilsat_inst.FillValue = default_fillvals['f4']
    soilsat_mean.FillValue = default_fillvals['f4']
    soilsat_anom.FillValue = default_fillvals['f4']
    soilsat_perc.FillValue = default_fillvals['f4']

    nav_lon.missing_value = 1.e+20
    nav_lat.missing_value = 1.e+20
    soilsat_inst.missing_value = 1.e+20
    soilsat_mean.missing_value = 1.e+20
    soilsat_anom.missing_value = 1.e+20
    soilsat_perc.missing_value = 1.e+20

    nav_lon.axis = 'YX'
    nav_lat.axis = 'YX'
    soilsat_inst.axis = 'YX'
    soilsat_mean.axis = 'YX'
    soilsat_anom.axis = 'YX'
    soilsat_perc.axis = 'YX'

    nav_lon.description = 'Longitude of grid cell center'
    nav_lat.description = 'Latitude of grid cell center'
    soilsat_inst.description = ('Simulated total column soil saturation for '
                              'a specific date')
    soilsat_mean.description = ('Long-term mean simulated total colunm soil '
                              'saturation for a specific day of the year. '
                              'The long-term mean is calculated as the 5 day '
                              'moving average centered on the current day. '
                              'The averaging pariod is 1916-2004')
    soilsat_anom.description = ('Total column soil saturation anomaly. Calculated '
                              'as soilsat_inst - soilsat_anom')
    soilsat_perc.description = ('Total column soil saturation percentile. This '
                              'value shows how often during the 1916-2004 '
                              'reference period the soil saturation on this '
                              'day of the year (using a 5 day centered '
                              'window) was less than soilsat_inst')

    nav_lon.valid_min = -180.
    nav_lat.valid_min = -90.
    soilsat_inst.valid_min = 0.
    soilsat_mean.valid_min = 0.
    soilsat_anom.valid_min = -1.
    soilsat_perc.valid_min = 0.

    nav_lon.valid_max = 180.
    nav_lat.valid_max = 90.
    soilsat_inst.valid_max = 1.
    soilsat_mean.valid_max = 1.
    soilsat_anom.valid_max = 1.
    soilsat_perc.valid_max = 1.

    nav_lon.modulo = 360.

    nc.history = 'Created ' + time.ctime(time.time())
    nc.description = ('Soil saturation data from the Variable Infiltration Model '
                      'as part of the operational surface water monitor. Note '
                      'that since this data is from the operational surface '
                      'water monitor, there may be occassional data problems. '
                      'Please check the data carefully')
    nc.source = ('Surface Water Monitor, Surface Water Hydrology Group, '
                 'University of Washington, Seattle, Washington, USA')
    nc.website = 'http://www.hydro.washington.edu/forecast/monitor'
    nc.contact = 'Bart Nijssen, email: [email protected]'
    nc.history = ' '.join(sys.argv)
    nc.projection = 'Geographic'
    nc.resolution = ('Spatial resolution: Longitude ({} degrees), latitude ({} '
#.........这里部分代码省略.........

# 需要导入模块: from netCDF4 import Dataset [as 别名]
# 或者: from netCDF4.Dataset import projection [as 别名]
mb = np.zeros_like(bed)


## acab_var = def_var(nc, "acab", "m year-1", fill_value)
## acab_var.standard_name = "land_ice_surface_specific_mass_balance"
## acab_var[:] = mb

# Set boundary conditions
# ------------------------------------------------------------------------------
#
# (A) Surface temperature for temperature equation bc
T0    = 273.15 # K
Tma   =  -6.0  # degC, mean annual air temperature at Tarfala
zcts  = 1300   # m a.s.l.; altitude where CTS is at the surface, projected to topg
zbts  = 1250   # m a.s.l.; altitude where CTS is at the base; just a wild guess

artm  = np.zeros((M,N),float) + T0
artm[bed<zcts] = T0 + Tma # Scandinavian-type polythermal glacier

artm_var = def_var(nc, "artm", "K", fill_value)
artm_var[:] = artm

# set global attributes
nc.Conventions = "CF-1.4"
historysep = ' '
historystr = time.asctime() + ': ' + historysep.join(sys.argv) + '\n'
setattr(nc, 'history', historystr)
nc.projection = projRT90
nc.close()
write('Done writing NetCDF file %s!\n' % ncfile)        

# 需要导入模块: from netCDF4 import Dataset [as 别名]
# 或者: from netCDF4.Dataset import projection [as 别名]
acab_var[:] = acab * ice_density

# Set boundary conditions for Scandinavian-type polythermal glacier
# ------------------------------------------------------------------------------
#
# (A) Surface temperature for temperature equation bc
T0 = 273.15  # K
Tma = -6.0  # degC, mean annual air temperature at Tarfala
zcts = 1300   # m a.s.l.; altitude where CTS is at the surface, projected to topg
slope = 100    # m; range around which surface temp transition happens


# smoothed version; FIXME:  can't we at least have it depend on initial DEM?
#   additional lapse rate?
artm = T0 + Tma * (zcts + slope - m_bed) / (2.0 * slope)
artm[m_bed < zcts-slope] = T0 + Tma
artm[m_bed > zcts+slope] = T0

artm_var = def_var(nc, "ice_surface_temp", "K", fill_value)
artm_var[:] = artm

# set global attributes
nc.Conventions = "CF-1.4"
historysep = ' '
historystr = time.asctime() + ': ' + historysep.join(sys.argv) + '\n'
setattr(nc, 'history', historystr)
# proj4 string equivalent to EPSG 3021
nc.projection = "+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"
nc.close()
write('Done writing NetCDF file %s!\n' % ncfile)