[7765] | 1 | |
---|
| 2 | ## |
---|
| 3 | # @brief Convert 3 URS files back to 4 NC files. |
---|
| 4 | # @param basename_in Stem of the input filenames. |
---|
| 5 | # @param basename_outStem of the output filenames. |
---|
| 6 | # @note The name of the urs file names must be: |
---|
| 7 | # [basename_in]-z-mux |
---|
| 8 | # [basename_in]-e-mux |
---|
| 9 | # [basename_in]-n-mux |
---|
| 10 | def urs2nc(basename_in='o', basename_out='urs'): |
---|
| 11 | """Convert the 3 urs files to 4 nc files. |
---|
| 12 | |
---|
| 13 | The name of the urs file names must be; |
---|
| 14 | [basename_in]-z-mux |
---|
| 15 | [basename_in]-e-mux |
---|
| 16 | [basename_in]-n-mux |
---|
| 17 | """ |
---|
| 18 | |
---|
| 19 | files_in = [basename_in + WAVEHEIGHT_MUX_LABEL, |
---|
| 20 | basename_in + EAST_VELOCITY_LABEL, |
---|
| 21 | basename_in + NORTH_VELOCITY_LABEL] |
---|
| 22 | files_out = [basename_out + '_ha.nc', |
---|
| 23 | basename_out + '_ua.nc', |
---|
| 24 | basename_out + '_va.nc'] |
---|
| 25 | quantities = ['HA', 'UA', 'VA'] |
---|
| 26 | |
---|
| 27 | #if os.access(files_in[0]+'.mux', os.F_OK) == 0 : |
---|
| 28 | for i, file_name in enumerate(files_in): |
---|
| 29 | if os.access(file_name, os.F_OK) == 0: |
---|
| 30 | if os.access(file_name + '.mux', os.F_OK) == 0 : |
---|
| 31 | msg = 'File %s does not exist or is not accessible' % file_name |
---|
| 32 | raise IOError, msg |
---|
| 33 | else: |
---|
| 34 | files_in[i] += '.mux' |
---|
| 35 | log.critical("file_name %s" % file_name) |
---|
| 36 | |
---|
| 37 | hashed_elevation = None |
---|
| 38 | for file_in, file_out, quantity in map(None, files_in, |
---|
| 39 | files_out, |
---|
| 40 | quantities): |
---|
| 41 | lonlatdep, lon, lat, depth = _binary_c2nc(file_in, |
---|
| 42 | file_out, |
---|
| 43 | quantity) |
---|
| 44 | if hashed_elevation == None: |
---|
| 45 | elevation_file = basename_out + '_e.nc' |
---|
| 46 | write_elevation_nc(elevation_file, |
---|
| 47 | lon, |
---|
| 48 | lat, |
---|
| 49 | depth) |
---|
| 50 | hashed_elevation = myhash(lonlatdep) |
---|
| 51 | else: |
---|
| 52 | msg = "The elevation information in the mux files is inconsistent" |
---|
| 53 | assert hashed_elevation == myhash(lonlatdep), msg |
---|
| 54 | |
---|
| 55 | files_out.append(elevation_file) |
---|
| 56 | |
---|
| 57 | return files_out |
---|
| 58 | |
---|
| 59 | |
---|
| 60 | ## |
---|
| 61 | # @brief Convert a quantity URS file to a NetCDF file. |
---|
| 62 | # @param file_in Path to input URS file. |
---|
| 63 | # @param file_out Path to the output file. |
---|
| 64 | # @param quantity Name of the quantity to be written to the output file. |
---|
| 65 | # @return A tuple (lonlatdep, lon, lat, depth). |
---|
| 66 | def _binary_c2nc(file_in, file_out, quantity): |
---|
| 67 | """Reads in a quantity urs file and writes a quantity nc file. |
---|
| 68 | Additionally, returns the depth and lat, long info, |
---|
| 69 | so it can be written to a file. |
---|
| 70 | """ |
---|
| 71 | |
---|
| 72 | columns = 3 # long, lat , depth |
---|
| 73 | mux_file = open(file_in, 'rb') |
---|
| 74 | |
---|
| 75 | # Number of points/stations |
---|
| 76 | (points_num,) = unpack('i', mux_file.read(4)) |
---|
| 77 | |
---|
| 78 | # nt, int - Number of time steps |
---|
| 79 | (time_step_count,) = unpack('i', mux_file.read(4)) |
---|
| 80 | |
---|
| 81 | #dt, float - time step, seconds |
---|
| 82 | (time_step,) = unpack('f', mux_file.read(4)) |
---|
| 83 | |
---|
| 84 | msg = "Bad data in the mux file." |
---|
| 85 | if points_num < 0: |
---|
| 86 | mux_file.close() |
---|
| 87 | raise ANUGAError, msg |
---|
| 88 | if time_step_count < 0: |
---|
| 89 | mux_file.close() |
---|
| 90 | raise ANUGAError, msg |
---|
| 91 | if time_step < 0: |
---|
| 92 | mux_file.close() |
---|
| 93 | raise ANUGAError, msg |
---|
| 94 | |
---|
| 95 | lonlatdep = p_array.array('f') |
---|
| 96 | lonlatdep.read(mux_file, columns * points_num) |
---|
| 97 | lonlatdep = num.array(lonlatdep, dtype=num.float) |
---|
| 98 | lonlatdep = num.reshape(lonlatdep, (points_num, columns)) |
---|
| 99 | |
---|
| 100 | lon, lat, depth = lon_lat2grid(lonlatdep) |
---|
| 101 | lon_sorted = list(lon) |
---|
| 102 | lon_sorted.sort() |
---|
| 103 | |
---|
| 104 | if not num.alltrue(lon == lon_sorted): |
---|
| 105 | msg = "Longitudes in mux file are not in ascending order" |
---|
| 106 | raise IOError, msg |
---|
| 107 | |
---|
| 108 | lat_sorted = list(lat) |
---|
| 109 | lat_sorted.sort() |
---|
| 110 | |
---|
| 111 | nc_file = Write_nc(quantity, |
---|
| 112 | file_out, |
---|
| 113 | time_step_count, |
---|
| 114 | time_step, |
---|
| 115 | lon, |
---|
| 116 | lat) |
---|
| 117 | |
---|
| 118 | for i in range(time_step_count): |
---|
| 119 | #Read in a time slice from mux file |
---|
| 120 | hz_p_array = p_array.array('f') |
---|
| 121 | hz_p_array.read(mux_file, points_num) |
---|
| 122 | hz_p = num.array(hz_p_array, dtype=num.float) |
---|
| 123 | hz_p = num.reshape(hz_p, (len(lon), len(lat))) |
---|
| 124 | hz_p = num.transpose(hz_p) # mux has lat varying fastest, nc has long v.f. |
---|
| 125 | |
---|
| 126 | #write time slice to nc file |
---|
| 127 | nc_file.store_timestep(hz_p) |
---|
| 128 | |
---|
| 129 | mux_file.close() |
---|
| 130 | nc_file.close() |
---|
| 131 | |
---|
| 132 | return lonlatdep, lon, lat, depth |
---|
| 133 | |
---|
| 134 | |
---|
[7776] | 135 | |
---|
| 136 | ## |
---|
| 137 | # @brief |
---|
| 138 | # @param long_lat_dep |
---|
| 139 | # @return A tuple (long, lat, quantity). |
---|
| 140 | # @note The latitude is the fastest varying dimension - in mux files. |
---|
| 141 | def lon_lat2grid(long_lat_dep): |
---|
| 142 | """ |
---|
| 143 | given a list of points that are assumed to be an a grid, |
---|
| 144 | return the long's and lat's of the grid. |
---|
| 145 | long_lat_dep is an array where each row is a position. |
---|
| 146 | The first column is longitudes. |
---|
| 147 | The second column is latitudes. |
---|
| 148 | |
---|
| 149 | The latitude is the fastest varying dimension - in mux files |
---|
| 150 | """ |
---|
| 151 | |
---|
| 152 | LONG = 0 |
---|
| 153 | LAT = 1 |
---|
| 154 | QUANTITY = 2 |
---|
| 155 | |
---|
| 156 | long_lat_dep = ensure_numeric(long_lat_dep, num.float) |
---|
| 157 | |
---|
| 158 | num_points = long_lat_dep.shape[0] |
---|
| 159 | this_rows_long = long_lat_dep[0,LONG] |
---|
| 160 | |
---|
| 161 | # Count the length of unique latitudes |
---|
| 162 | i = 0 |
---|
| 163 | while long_lat_dep[i,LONG] == this_rows_long and i < num_points: |
---|
| 164 | i += 1 |
---|
| 165 | |
---|
| 166 | # determine the lats and longsfrom the grid |
---|
| 167 | lat = long_lat_dep[:i, LAT] |
---|
| 168 | long = long_lat_dep[::i, LONG] |
---|
| 169 | |
---|
| 170 | lenlong = len(long) |
---|
| 171 | lenlat = len(lat) |
---|
| 172 | |
---|
| 173 | msg = 'Input data is not gridded' |
---|
| 174 | assert num_points % lenlat == 0, msg |
---|
| 175 | assert num_points % lenlong == 0, msg |
---|
| 176 | |
---|
| 177 | # Test that data is gridded |
---|
| 178 | for i in range(lenlong): |
---|
| 179 | msg = 'Data is not gridded. It must be for this operation' |
---|
| 180 | first = i * lenlat |
---|
| 181 | last = first + lenlat |
---|
| 182 | |
---|
| 183 | assert num.allclose(long_lat_dep[first:last,LAT], lat), msg |
---|
| 184 | assert num.allclose(long_lat_dep[first:last,LONG], long[i]), msg |
---|
| 185 | |
---|
| 186 | msg = 'Out of range latitudes/longitudes' |
---|
| 187 | for l in lat:assert -90 < l < 90 , msg |
---|
| 188 | for l in long:assert -180 < l < 180 , msg |
---|
| 189 | |
---|
| 190 | # Changing quantity from lat being the fastest varying dimension to |
---|
| 191 | # long being the fastest varying dimension |
---|
| 192 | # FIXME - make this faster/do this a better way |
---|
| 193 | # use numeric transpose, after reshaping the quantity vector |
---|
| 194 | quantity = num.zeros(num_points, num.float) |
---|
| 195 | |
---|
| 196 | for lat_i, _ in enumerate(lat): |
---|
| 197 | for long_i, _ in enumerate(long): |
---|
| 198 | q_index = lat_i*lenlong + long_i |
---|
| 199 | lld_index = long_i*lenlat + lat_i |
---|
| 200 | temp = long_lat_dep[lld_index, QUANTITY] |
---|
| 201 | quantity[q_index] = temp |
---|
| 202 | |
---|
| 203 | return long, lat, quantity |
---|
| 204 | |
---|
| 205 | |
---|
| 206 | |
---|