source: trunk/anuga_core/source/anuga/file/netcdf.py @ 7800

Last change on this file since 7800 was 7800, checked in by hudson, 12 years ago

urs2sww has an extra urs_ungridded2sww function.

File size: 9.1 KB
Line 
1""" This module is responsible for loading and saving NetCDF NC files
2"""
3
4import os, sys
5import csv
6import string
7import shutil
8from struct import unpack
9import array as p_array
10from os import sep, path, remove, mkdir, access, F_OK, W_OK, getcwd
11
12import numpy as num
13
14from Scientific.IO.NetCDF import NetCDFFile
15from os.path import exists, basename, join
16from os import getcwd
17
18from anuga.coordinate_transforms.redfearn import redfearn, \
19     convert_from_latlon_to_utm
20from anuga.coordinate_transforms.geo_reference import Geo_reference, \
21     write_NetCDF_georeference, ensure_geo_reference
22from anuga.geospatial_data.geospatial_data import Geospatial_data,\
23     ensure_absolute
24from anuga.config import minimum_storable_height as \
25     default_minimum_storable_height
26from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
27from anuga.config import netcdf_float, netcdf_float32, netcdf_int
28from anuga.config import max_float
29from anuga.utilities.numerical_tools import ensure_numeric,  mean
30from anuga.caching.caching import myhash
31from anuga.shallow_water.shallow_water_domain import Domain
32from anuga.abstract_2d_finite_volumes.pmesh2domain import \
33     pmesh_to_domain_instance
34from anuga.abstract_2d_finite_volumes.util import get_revision_number, \
35     remove_lone_verts, sww2timeseries, get_centroid_values
36
37from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints
38from anuga.load_mesh.loadASCII import export_mesh_file
39from anuga.geometry.polygon import intersection
40from anuga.file_conversion.sww2dem import sww2dem
41
42from anuga.utilities.system_tools import get_vars_in_expression
43import anuga.utilities.log as log
44
45from anuga.utilities.file_utils import create_filename,\
46                        get_all_swwfiles
47from anuga.file.csv_file import load_csv_as_dict
48from sww import Read_sww, Write_sww
49
50from anuga.anuga_exceptions import DataMissingValuesError, \
51                DataFileNotOpenError, DataTimeError, DataDomainError, \
52                NewQuantity
53
54
55# Definitions of various NetCDF dimension names, etc.
56lon_name = 'LON'
57lat_name = 'LAT'
58time_name = 'TIME'
59precision = netcdf_float # So if we want to change the precision its done here
60
61##
62# @brief Clas for a NetCDF data file writer.
63class Write_nc:
64    """Write an nc file.
65
66    Note, this should be checked to meet cdc netcdf conventions for gridded
67    data. http://www.cdc.noaa.gov/cdc/conventions/cdc_netcdf_standard.shtml
68    """
69
70    ##
71    # @brief Instantiate a Write_nc instance.
72    # @param quantity_name
73    # @param file_name
74    # @param time_step_count The number of time steps.
75    # @param time_step The time_step size.
76    # @param lon
77    # @param lat
78    def __init__(self,
79                 quantity_name,
80                 file_name,
81                 time_step_count,
82                 time_step,
83                 lon,
84                 lat):
85        """Instantiate a Write_nc instance (NetCDF file writer).
86
87        time_step_count is the number of time steps.
88        time_step is the time step size
89
90        pre-condition: quantity_name must be 'HA', 'UA'or 'VA'.
91        """
92
93        self.quantity_name = quantity_name
94        quantity_units = {'HA':'CENTIMETERS',
95                          'UA':'CENTIMETERS/SECOND',
96                          'VA':'CENTIMETERS/SECOND'}
97
98        multiplier_dic = {'HA':100.0,   # To convert from m to cm
99                          'UA':100.0,   #             and m/s to cm/sec
100                          'VA':-100.0}  # MUX files have positive x in the
101                                        # Southern direction.  This corrects
102                                        # for it, when writing nc files.
103
104        self.quantity_multiplier =  multiplier_dic[self.quantity_name]
105
106        #self.file_name = file_name
107        self.time_step_count = time_step_count
108        self.time_step = time_step
109
110        # NetCDF file definition
111        self.outfile = NetCDFFile(file_name, netcdf_mode_w)
112        outfile = self.outfile
113
114        #Create new file
115        nc_lon_lat_header(outfile, lon, lat)
116
117        # TIME
118        outfile.createDimension(time_name, None)
119        outfile.createVariable(time_name, precision, (time_name,))
120
121        #QUANTITY
122        outfile.createVariable(self.quantity_name, precision,
123                               (time_name, lat_name, lon_name))
124        outfile.variables[self.quantity_name].missing_value = -1.e+034
125        outfile.variables[self.quantity_name].units = \
126                                 quantity_units[self.quantity_name]
127        outfile.variables[lon_name][:]= ensure_numeric(lon)
128        outfile.variables[lat_name][:]= ensure_numeric(lat)
129
130        #Assume no one will be wanting to read this, while we are writing
131        #outfile.close()
132
133    ##
134    # @brief Write a time-step of quantity data.
135    # @param quantity_slice The data to be stored for this time-step.
136    def store_timestep(self, quantity_slice):
137        """Write a time slice of quantity info
138
139        quantity_slice is the data to be stored at this time step
140        """
141
142        # Get the variables
143        time = self.outfile.variables[time_name]
144        quantity = self.outfile.variables[self.quantity_name]
145
146        # get index oflice to write
147        i = len(time)
148
149        #Store time
150        time[i] = i * self.time_step    #self.domain.time
151        quantity[i,:] = quantity_slice * self.quantity_multiplier
152
153    ##
154    # @brief Close file underlying the class instance.
155    def close(self):
156        self.outfile.close()
157
158
159
160##
161# @brief Write an NC elevation file.
162# @param file_out Path to the output file.
163# @param lon ??
164# @param lat ??
165# @param depth_vector The elevation data to write.
166def write_elevation_nc(file_out, lon, lat, depth_vector):
167    """Write an nc elevation file."""
168
169    # NetCDF file definition
170    outfile = NetCDFFile(file_out, netcdf_mode_w)
171
172    #Create new file
173    nc_lon_lat_header(outfile, lon, lat)
174
175    # ELEVATION
176    zname = 'ELEVATION'
177    outfile.createVariable(zname, precision, (lat_name, lon_name))
178    outfile.variables[zname].units = 'CENTIMETERS'
179    outfile.variables[zname].missing_value = -1.e+034
180
181    outfile.variables[lon_name][:] = ensure_numeric(lon)
182    outfile.variables[lat_name][:] = ensure_numeric(lat)
183
184    depth = num.reshape(depth_vector, (len(lat), len(lon)))
185    outfile.variables[zname][:] = depth
186
187    outfile.close()
188
189
190##
191# @brief Write lat/lon headers to a NetCDF file.
192# @param outfile Handle to open file to write to.
193# @param lon An iterable of the longitudes.
194# @param lat An iterable of the latitudes.
195# @note Defines lat/long dimensions and variables. Sets various attributes:
196#          .point_spacing  and  .units
197#       and writes lat/lon data.
198
199def nc_lon_lat_header(outfile, lon, lat):
200    """Write lat/lon headers to a NetCDF file.
201
202    outfile is the netcdf file handle.
203    lon - a list/array of the longitudes
204    lat - a list/array of the latitudes
205    """
206
207    outfile.institution = 'Geoscience Australia'
208    outfile.description = 'Converted from URS binary C'
209
210    # Longitude
211    outfile.createDimension(lon_name, len(lon))
212    outfile.createVariable(lon_name, precision, (lon_name,))
213    outfile.variables[lon_name].point_spacing = 'uneven'
214    outfile.variables[lon_name].units = 'degrees_east'
215    outfile.variables[lon_name].assignValue(lon)
216
217    # Latitude
218    outfile.createDimension(lat_name, len(lat))
219    outfile.createVariable(lat_name, precision, (lat_name,))
220    outfile.variables[lat_name].point_spacing = 'uneven'
221    outfile.variables[lat_name].units = 'degrees_north'
222    outfile.variables[lat_name].assignValue(lat)
223
224
225
226
227##
228# @brief Filter data file, selecting timesteps first:step:last.
229# @param filename1 Data file to filter.
230# @param filename2 File to write filtered timesteps to.
231# @param first First timestep.
232# @param last Last timestep.
233# @param step Timestep stride.
234def filter_netcdf(filename1, filename2, first=0, last=None, step=1):
235    """Filter data file, selecting timesteps first:step:last.
236   
237    Read netcdf filename1, pick timesteps first:step:last and save to
238    nettcdf file filename2
239    """
240
241    from Scientific.IO.NetCDF import NetCDFFile
242
243    # Get NetCDF
244    infile = NetCDFFile(filename1, netcdf_mode_r)  #Open existing file for read
245    outfile = NetCDFFile(filename2, netcdf_mode_w)  #Open new file
246
247    # Copy dimensions
248    for d in infile.dimensions:
249        outfile.createDimension(d, infile.dimensions[d])
250
251    # Copy variable definitions
252    for name in infile.variables:
253        var = infile.variables[name]
254        outfile.createVariable(name, var.dtype.char, var.dimensions)
255
256    # Copy the static variables
257    for name in infile.variables:
258        if name == 'time' or name == 'stage':
259            pass
260        else:
261            outfile.variables[name][:] = infile.variables[name][:]
262
263    # Copy selected timesteps
264    time = infile.variables['time']
265    stage = infile.variables['stage']
266
267    newtime = outfile.variables['time']
268    newstage = outfile.variables['stage']
269
270    if last is None:
271        last = len(time)
272
273    selection = range(first, last, step)
274    for i, j in enumerate(selection):
275        log.critical('Copying timestep %d of %d (%f)'
276                     % (j, last-first, time[j]))
277        newtime[i] = time[j]
278        newstage[i,:] = stage[j,:]
279
280    # Close
281    infile.close()
282    outfile.close()
283
Note: See TracBrowser for help on using the repository browser.