Changeset 5534
- Timestamp:
- Jul 18, 2008, 3:14:35 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/data_manager.py
r5485 r5534 4841 4841 NORTH_VELOCITY_MUX2_LABEL = '-n-mux2' 4842 4842 4843 def read_mux2_py(filenames,weights,verbose=False): 4843 def read_mux2_py(filenames, 4844 weights, 4845 permutation=None, 4846 verbose=False): 4847 """Access the mux files specified in the filenames list. Combine the 4848 data found therin as a weighted linear sum as specifed by the weights. 4849 If permutation is None extract timeseries data for all gauges within the 4850 files. 4851 4852 Input: 4853 filenames: List of filenames specifiying the file containing the 4854 timeseries data (mux2 format) for each source 4855 weights: Weights associated with each source 4856 permutation: Specifies the gauge numbers that for which data is to be 4857 extracted 4858 """ 4844 4859 4845 4860 from Numeric import ones,Float,compress,zeros,arange … … 4856 4871 verbose=0 4857 4872 4873 # Call underlying C implementation urs2sts_ext.c 4858 4874 data=read_mux2(numSrc,filenames,weights,file_params,verbose) 4859 4875 … … 4910 4926 4911 4927 4912 def urs2sts(basename_in, basename_out = None, weights=None, 4913 verbose = False, origin = None, 4914 mean_stage=0.0, zscale=1.0, 4915 ordering_filename = None, 4916 minlat = None, maxlat = None, 4917 minlon = None, maxlon = None): 4928 def urs2sts(basename_in, basename_out=None, 4929 weights=None, 4930 verbose=False, 4931 origin=None, 4932 mean_stage=0.0, 4933 zscale=1.0, 4934 ordering_filename=None): 4918 4935 """Convert URS mux2 format for wave propagation to sts format 4919 4936 … … 4924 4941 UTM coordinates (zone, easting, northing) 4925 4942 4926 input: 4943 inputs: 4944 4927 4945 basename_in: list of source file prefixes 4928 4946 4929 These are combined with the extensions: 4930 WAVEHEIGHT_MUX2_LABEL = '-z-mux2' for stage 4931 EAST_VELOCITY_MUX2_LABEL = '-e-mux2' xmomentum 4932 NORTH_VELOCITY_MUX2_LABEL = '-n-mux2' and ymomentum 4933 4934 to create a 2D list of mux2 file. The rows are associated with each quantity and must have the above extensions 4935 the columns are the list of file prefixes. 4936 4937 ordering: a .txt file name specifying which mux2 gauge points are sto be stored. This is indicated by the index of 4938 the gauge in the ordering file. 4939 4940 ordering file format: 4941 header 4942 index,longitude,latitude 4943 4944 header='index,longitude,latitude\n' 4945 If ordering is None all points are taken in the order they appear in the mux2 file subject to rectangular clipping box (see below). 4946 4947 4948 minlat, minlon, maxlat, maxlon define a rectangular clipping box and can be used to extract gauge information in a rectangular box. This is useful for extracting 4949 information at gauge points not on the boundary. 4950 if ordering_filename is specified clipping box cannot be used 4947 These are combined with the extensions: 4948 WAVEHEIGHT_MUX2_LABEL = '-z-mux2' for stage 4949 EAST_VELOCITY_MUX2_LABEL = '-e-mux2' xmomentum 4950 NORTH_VELOCITY_MUX2_LABEL = '-n-mux2' and ymomentum 4951 4952 to create a 2D list of mux2 file. The rows are associated with each 4953 quantity and must have the above extensions 4954 the columns are the list of file prefixes. 4955 4956 ordering: a .txt file name specifying which mux2 gauge points are 4957 to be stored. This is indicated by the index of the gauge 4958 in the ordering file. 4959 4960 ordering file format: 4961 1st line: 'index,longitude,latitude\n' 4962 other lines: index,longitude,latitude 4963 4964 4965 If ordering is None all points are taken in the order they 4966 appear in the mux2 file. 4951 4967 4952 4968 4953 4969 output: 4954 baename_out: name of sts file in which mux2 data is stored.4970 basename_out: name of sts file in which mux2 data is stored. 4955 4971 4956 4972 … … 4968 4984 4969 4985 # Check that basename is a list of strings 4970 if not reduce(__and__, map(lambda z:isinstance(z,StringType), basename_in)):4986 if not reduce(__and__, map(lambda z:isinstance(z,StringType), basename_in)): 4971 4987 msg= 'basename_in must be a string or list of strings' 4972 4988 raise Exception, msg … … 4977 4993 # A weight must be specified for each source 4978 4994 if weights is None: 4995 # Default is equal weighting 4979 4996 weights=ones(numSrc,Float)/numSrc 4980 4997 else: … … 4984 5001 assert len(weights)== numSrc, msg 4985 5002 4986 precision = Float4987 4988 if ordering_filename is not None:4989 msg = 'If ordering_filename is specified,'4990 msg += ' rectangular clipping box cannot be specified as well. '4991 assert minlat is None and maxlat is None and\4992 minlon is None and maxlon is None, msg4993 4994 4995 msg = 'Must use latitudes and longitudes for minlat, maxlon etc'4996 if minlat != None:4997 assert -90 < minlat < 90 , msg4998 if maxlat != None:4999 assert -90 < maxlat < 90 , msg5000 if minlat != None:5001 assert maxlat > minlat5002 if minlon != None:5003 assert -180 < minlon < 180 , msg5004 if maxlon != None:5005 assert -180 < maxlon < 180 , msg5006 if minlon != None:5007 assert maxlon > minlon5008 5009 5010 5003 # Check output filename 5011 5004 if basename_out is None: 5012 5005 msg = 'STS filename must be specified' 5013 5006 raise Exception, msg 5014 5015 5007 stsname = basename_out + '.sts' 5016 5008 5009 # Create input filenames from basenames and check their existence 5017 5010 files_in=[[],[],[]] 5018 5011 for files in basename_in: … … 5020 5013 files_in[1].append(files + EAST_VELOCITY_MUX2_LABEL) 5021 5014 files_in[2].append(files + NORTH_VELOCITY_MUX2_LABEL) 5022 5023 #files_in = [basename_in + WAVEHEIGHT_MUX2_LABEL, 5024 # basename_in + EAST_VELOCITY_MUX2_LABEL, 5025 # basename_in + NORTH_VELOCITY_MUX2_LABEL] 5026 5027 quantities = ['HA','UA','VA'] 5028 5029 # For each source file check that there exists three files ending with: 5030 # WAVEHEIGHT_MUX2_LABEL, 5031 # EAST_VELOCITY_MUX2_LABEL, and 5032 # NORTH_VELOCITY_MUX2_LABEL 5015 5016 quantities = ['HA','UA','VA'] # Quantity names used in the MUX2 format 5033 5017 for i in range(len(quantities)): 5034 5018 for file_in in files_in[i]: … … 5037 5021 raise IOError, msg 5038 5022 5039 #need to do this for velocity-e-mux2 and velocity-n-mux2 files as well 5040 #for now set x_momentum and y_moentum quantities to zero 5023 # Read MUX2 files 5041 5024 if (verbose): print 'reading mux2 file' 5042 5025 mux={} 5043 for i,quantity in enumerate(quantities): 5044 # For each quantity read the associated list of source mux2 file with extenstion associated with that quantity 5045 times, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights,verbose=verbose) 5026 for i, quantity in enumerate(quantities): 5027 5028 # For each quantity read the associated list of source mux2 file with 5029 # extention associated with that quantity 5030 times,\ 5031 latitudes_urs,\ 5032 longitudes_urs,\ 5033 elevation,\ 5034 mux[quantity],\ 5035 starttime = read_mux2_py(files_in[i], weights, verbose=verbose) 5036 5037 # Check that all quantities have consistent time and space information 5046 5038 if quantity!=quantities[0]: 5047 msg='%s, %s and %s have inconsitent gauge data'%(files_in[0],files_in[1],files_in[2]) 5048 assert allclose(times,times_old),msg 5049 assert allclose(latitudes_urs,latitudes_urs_old),msg 5050 assert allclose(longitudes_urs,longitudes_urs_old),msg 5051 assert allclose(elevation,elevation_old),msg 5052 assert allclose(starttime,starttime_old) 5039 msg='%s, %s and %s have inconsistent gauge data'\ 5040 %(files_in[0], files_in[1], files_in[2]) 5041 assert allclose(times, times_old), msg 5042 assert allclose(latitudes_urs, latitudes_urs_old), msg 5043 assert allclose(longitudes_urs, longitudes_urs_old), msg 5044 assert allclose(elevation, elevation_old), msg 5045 assert allclose(starttime, starttime_old), msg 5053 5046 times_old=times 5054 5047 latitudes_urs_old=latitudes_urs … … 5104 5097 assert allclose(latitudes, ref_latitudes), msg 5105 5098 5106 elif (minlat is not None) and (minlon is not None) and (maxlat is not None) and (maxlon is not None):5107 #FIXME - this branch is probably not working5108 if verbose: print 'Cliping urs data'5109 latitudes = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),latitudes_urs)5110 longitudes = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),longitudes_urs)5111 5112 permutation = range(latitudes.shape[0])5113 5114 msg = 'Clipping is not done yet'5115 raise Exception, msg5116 5117 5099 else: 5118 5100 latitudes=latitudes_urs … … 5120 5102 permutation = range(latitudes.shape[0]) 5121 5103 5122 5104 5105 5106 # Store timeseries in NetCDF STS file 5123 5107 msg='File is empty and or clipped region not in file region' 5124 5108 assert len(latitudes>0),msg
Note: See TracChangeset
for help on using the changeset viewer.