[2465] | 1 | """Class Geospatial_data - Manipulation of locations on the planet and |
---|
| 2 | associated attributes. |
---|
[2303] | 3 | |
---|
| 4 | """ |
---|
| 5 | |
---|
[3280] | 6 | from os import access, F_OK, R_OK |
---|
[3738] | 7 | from types import DictType |
---|
| 8 | |
---|
[3702] | 9 | from Numeric import concatenate, array, Float, shape, reshape, ravel, take |
---|
[2303] | 10 | |
---|
[3514] | 11 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
| 12 | from anuga.coordinate_transforms.geo_reference import Geo_reference, TitleError |
---|
[3739] | 13 | from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm |
---|
[2570] | 14 | |
---|
[2303] | 15 | |
---|
| 16 | class Geospatial_data: |
---|
| 17 | |
---|
| 18 | def __init__(self, |
---|
[2570] | 19 | data_points = None, |
---|
[2303] | 20 | attributes = None, |
---|
| 21 | geo_reference = None, |
---|
[2570] | 22 | default_attribute_name = None, |
---|
[2928] | 23 | file_name = None, |
---|
[2935] | 24 | delimiter = None, |
---|
[3280] | 25 | latitudes = None, |
---|
| 26 | longitudes = None, |
---|
[3297] | 27 | points_are_lats_longs = False, |
---|
[2935] | 28 | verbose = False): |
---|
[2303] | 29 | |
---|
[2643] | 30 | |
---|
[2928] | 31 | """ |
---|
[2643] | 32 | Create instance from data points and associated attributes |
---|
[2303] | 33 | |
---|
| 34 | data_points: x,y coordinates in meters. Type must be either a |
---|
| 35 | sequence of 2-tuples or an Mx2 Numeric array of floats. |
---|
| 36 | |
---|
| 37 | attributes: Associated values for each data point. The type |
---|
| 38 | must be either a list or an array of length M or a dictionary |
---|
| 39 | of lists (or arrays) of length M. In the latter case the keys |
---|
| 40 | in the dictionary represent the attribute names, in the former |
---|
[2889] | 41 | the attribute will get the default name "attribute". |
---|
[2303] | 42 | |
---|
| 43 | geo_reference: Object representing the origin of the data |
---|
| 44 | points. It contains UTM zone, easting and northing and data |
---|
| 45 | points are assumed to be relative to this origin. |
---|
| 46 | If geo_reference is None, the default geo ref object is used |
---|
| 47 | |
---|
| 48 | default_attribute_name: Name of default attribute to be used with |
---|
| 49 | get_attribute_values. The idea is that the dataset can be |
---|
| 50 | equipped with information about which attribute to return. |
---|
[2889] | 51 | If None, the default is the "first" |
---|
[2303] | 52 | |
---|
[2643] | 53 | file_name: Name of input netCDF file or xya file. netCDF file must |
---|
| 54 | have dimensions "points" etc. |
---|
[3451] | 55 | xya file is a comma seperated file with x, y and attribute data. |
---|
| 56 | the first line must be the attribute names eg elevation |
---|
[2570] | 57 | |
---|
[2698] | 58 | The format for a .xya file is: |
---|
| 59 | 1st line: [attribute names] |
---|
| 60 | other lines: x y [attributes] |
---|
| 61 | |
---|
| 62 | for example: |
---|
| 63 | elevation, friction |
---|
| 64 | 0.6, 0.7, 4.9, 0.3 |
---|
| 65 | 1.9, 2.8, 5, 0.3 |
---|
| 66 | 2.7, 2.4, 5.2, 0.3 |
---|
| 67 | |
---|
| 68 | The first two columns are always implicitly assumed to be x, y coordinates. |
---|
| 69 | Use the same delimiter for the attribute names and the data |
---|
| 70 | |
---|
| 71 | An xya file can optionally end with |
---|
| 72 | #geo reference |
---|
| 73 | 56 |
---|
| 74 | 466600.0 |
---|
| 75 | 8644444.0 |
---|
| 76 | |
---|
| 77 | When the 1st # is the zone, |
---|
| 78 | 2nd # the xllcorner and |
---|
| 79 | 3rd # the yllcorner |
---|
[3451] | 80 | |
---|
| 81 | An issue with the xya format is that the attribute column order |
---|
| 82 | is not be controlled. The info is stored in a dictionary and it's |
---|
| 83 | written however |
---|
[2624] | 84 | |
---|
[2698] | 85 | The format for a Points dictionary is: |
---|
| 86 | |
---|
| 87 | ['pointlist'] a 2 column array describing points. 1st column x, 2nd column y. |
---|
| 88 | ['attributelist'], a dictionary of 1D arrays, representing attribute values |
---|
| 89 | at the point. The dictionary key is the attribute header. |
---|
| 90 | ['geo_reference'] a Geo_refernece object. Use if the point information |
---|
| 91 | is relative. This is optional. |
---|
| 92 | eg |
---|
| 93 | dic['pointlist'] = [[1.0,2.0],[3.0,5.0]] |
---|
| 94 | dic['attributelist']['elevation'] = [[7.0,5.0] |
---|
[2928] | 95 | |
---|
[2940] | 96 | delimiter: is the file delimiter that will be used when |
---|
| 97 | importing the file |
---|
[2928] | 98 | |
---|
| 99 | verbose: |
---|
[2698] | 100 | |
---|
[2928] | 101 | """ |
---|
[2303] | 102 | |
---|
[2711] | 103 | if isinstance(data_points, basestring): |
---|
| 104 | # assume data point is really a file name |
---|
| 105 | file_name = data_points |
---|
| 106 | |
---|
[2940] | 107 | self.set_verbose(verbose) |
---|
[3149] | 108 | self.geo_reference=None #create the attribute |
---|
[2570] | 109 | if file_name is None: |
---|
[2928] | 110 | if delimiter is not None: |
---|
| 111 | msg = 'No file specified yet a delimiter is provided!' |
---|
| 112 | raise ValueError, msg |
---|
[2940] | 113 | file_name = None |
---|
[3297] | 114 | if latitudes is not None or longitudes is not None or \ |
---|
| 115 | points_are_lats_longs: |
---|
[3280] | 116 | data_points, geo_reference = \ |
---|
| 117 | self._set_using_lat_long(latitudes=latitudes, |
---|
[3297] | 118 | longitudes=longitudes, |
---|
| 119 | geo_reference=geo_reference, |
---|
| 120 | data_points=data_points, |
---|
| 121 | points_are_lats_longs=points_are_lats_longs) |
---|
[2570] | 122 | self.check_data_points(data_points) |
---|
| 123 | self.set_attributes(attributes) |
---|
| 124 | self.set_geo_reference(geo_reference) |
---|
| 125 | self.set_default_attribute_name(default_attribute_name) |
---|
[2928] | 126 | |
---|
[2570] | 127 | else: |
---|
[3738] | 128 | # watch for case where file name and points, attributes etc are provided!! |
---|
| 129 | # if file name then all provided info will be removed! |
---|
| 130 | self.import_points_file(file_name, delimiter, verbose) |
---|
[2570] | 131 | |
---|
[3738] | 132 | self.check_data_points(self.data_points) |
---|
| 133 | self.set_attributes(self.attributes) |
---|
| 134 | self.set_geo_reference(self.geo_reference) |
---|
| 135 | self.set_default_attribute_name(default_attribute_name) |
---|
[2303] | 136 | |
---|
[3738] | 137 | |
---|
| 138 | assert self.attributes is None or isinstance(self.attributes, DictType) |
---|
| 139 | |
---|
| 140 | |
---|
[3302] | 141 | def __len__(self): |
---|
| 142 | return len(self.data_points) |
---|
[3702] | 143 | |
---|
| 144 | def __repr__(self): |
---|
| 145 | return str(self.get_data_points(absolute=True)) |
---|
[3302] | 146 | |
---|
[3702] | 147 | |
---|
[2570] | 148 | def check_data_points(self, data_points): |
---|
| 149 | """Checks data points |
---|
| 150 | """ |
---|
| 151 | |
---|
| 152 | if data_points is None: |
---|
| 153 | self.data_points = None |
---|
[2928] | 154 | msg = 'There is no data or file provided!' |
---|
| 155 | raise ValueError, msg |
---|
[2570] | 156 | |
---|
| 157 | else: |
---|
| 158 | self.data_points = ensure_numeric(data_points) |
---|
[2711] | 159 | assert len(self.data_points.shape) == 2 |
---|
| 160 | assert self.data_points.shape[1] == 2 |
---|
[2303] | 161 | |
---|
| 162 | def set_attributes(self, attributes): |
---|
| 163 | """Check and assign attributes dictionary |
---|
| 164 | """ |
---|
| 165 | |
---|
| 166 | if attributes is None: |
---|
| 167 | self.attributes = None |
---|
| 168 | return |
---|
| 169 | |
---|
| 170 | if not isinstance(attributes, DictType): |
---|
| 171 | #Convert single attribute into dictionary |
---|
| 172 | attributes = {'attribute': attributes} |
---|
| 173 | |
---|
| 174 | #Check input attributes |
---|
| 175 | for key in attributes.keys(): |
---|
| 176 | try: |
---|
| 177 | attributes[key] = ensure_numeric(attributes[key]) |
---|
| 178 | except: |
---|
| 179 | msg = 'Attribute %s could not be converted' %key |
---|
| 180 | msg += 'to a numeric vector' |
---|
| 181 | raise msg |
---|
| 182 | |
---|
| 183 | self.attributes = attributes |
---|
| 184 | |
---|
| 185 | |
---|
| 186 | def set_geo_reference(self, geo_reference): |
---|
| 187 | |
---|
[3514] | 188 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
[2303] | 189 | |
---|
| 190 | if geo_reference is None: |
---|
[3149] | 191 | geo_reference = Geo_reference() # Use default |
---|
| 192 | if not isinstance(geo_reference, Geo_reference): |
---|
[2489] | 193 | msg = 'Argument geo_reference must be a valid Geo_reference \n' |
---|
| 194 | msg += 'object or None.' |
---|
[2303] | 195 | raise msg |
---|
| 196 | |
---|
[3149] | 197 | # if a geo ref already exists, change the point data to |
---|
| 198 | # represent the new geo-ref |
---|
| 199 | if self.geo_reference is not None: |
---|
| 200 | #FIXME: Maybe put out a warning here... |
---|
[3154] | 201 | self.data_points = self.get_data_points \ |
---|
| 202 | (geo_reference=geo_reference) |
---|
[3149] | 203 | |
---|
| 204 | self.geo_reference = geo_reference |
---|
[2303] | 205 | |
---|
[3149] | 206 | |
---|
[2309] | 207 | def set_default_attribute_name(self, default_attribute_name): |
---|
| 208 | self.default_attribute_name = default_attribute_name |
---|
[2303] | 209 | |
---|
[2935] | 210 | def set_verbose(self, verbose = False): |
---|
| 211 | if verbose is not False: |
---|
| 212 | verbose = True |
---|
| 213 | else: |
---|
| 214 | verbose = False |
---|
[3702] | 215 | |
---|
| 216 | def clip(self, polygon, closed=True): |
---|
[3712] | 217 | """Clip geospatial data by a polygon |
---|
[3702] | 218 | |
---|
| 219 | Input |
---|
| 220 | polygon - Either a list of points, an Nx2 array or |
---|
| 221 | a Geospatial data object. |
---|
| 222 | closed - (optional) determine whether points on boundary should be |
---|
| 223 | regarded as belonging to the polygon (closed = True) |
---|
| 224 | or not (closed = False). Default is True. |
---|
| 225 | |
---|
| 226 | Output |
---|
[3712] | 227 | New geospatial data object representing points inside |
---|
| 228 | specified polygon. |
---|
[3702] | 229 | |
---|
| 230 | """ |
---|
| 231 | |
---|
| 232 | from anuga.utilities.polygon import inside_polygon |
---|
| 233 | |
---|
| 234 | if isinstance(polygon, Geospatial_data): |
---|
| 235 | # Polygon is an object - extract points |
---|
| 236 | polygon = polygon.get_data_points() |
---|
| 237 | |
---|
| 238 | points = self.get_data_points() |
---|
| 239 | inside_indices = inside_polygon(points, polygon, closed) |
---|
[3738] | 240 | |
---|
| 241 | clipped_points = take(points, inside_indices) |
---|
| 242 | |
---|
| 243 | # Clip all attributes |
---|
| 244 | attributes = self.get_all_attributes() |
---|
| 245 | |
---|
| 246 | clipped_attributes = {} |
---|
| 247 | if attributes is not None: |
---|
| 248 | for key, att in attributes.items(): |
---|
| 249 | clipped_attributes[key] = take(att, inside_indices) |
---|
| 250 | |
---|
| 251 | return Geospatial_data(clipped_points, clipped_attributes) |
---|
[3702] | 252 | |
---|
| 253 | |
---|
[3714] | 254 | def clip_outside(self, polygon, closed=True): |
---|
| 255 | """Clip geospatial date by a polygon, keeping data OUTSIDE of polygon |
---|
| 256 | |
---|
| 257 | Input |
---|
| 258 | polygon - Either a list of points, an Nx2 array or |
---|
| 259 | a Geospatial data object. |
---|
| 260 | closed - (optional) determine whether points on boundary should be |
---|
| 261 | regarded as belonging to the polygon (closed = True) |
---|
| 262 | or not (closed = False). Default is True. |
---|
| 263 | |
---|
| 264 | Output |
---|
| 265 | Geospatial data object representing point OUTSIDE specified polygon |
---|
| 266 | |
---|
| 267 | """ |
---|
| 268 | |
---|
| 269 | from anuga.utilities.polygon import outside_polygon |
---|
| 270 | |
---|
| 271 | if isinstance(polygon, Geospatial_data): |
---|
| 272 | # Polygon is an object - extract points |
---|
| 273 | polygon = polygon.get_data_points() |
---|
| 274 | |
---|
| 275 | points = self.get_data_points() |
---|
| 276 | outside_indices = outside_polygon(points, polygon, closed) |
---|
[3738] | 277 | |
---|
| 278 | clipped_points = take(points, outside_indices) |
---|
| 279 | |
---|
| 280 | # Clip all attributes |
---|
| 281 | attributes = self.get_all_attributes() |
---|
| 282 | |
---|
| 283 | clipped_attributes = {} |
---|
| 284 | if attributes is not None: |
---|
| 285 | for key, att in attributes.items(): |
---|
| 286 | clipped_attributes[key] = take(att, outside_indices) |
---|
| 287 | |
---|
| 288 | return Geospatial_data(clipped_points, clipped_attributes) |
---|
| 289 | |
---|
[3714] | 290 | |
---|
[3280] | 291 | def _set_using_lat_long(self, |
---|
| 292 | latitudes, |
---|
| 293 | longitudes, |
---|
| 294 | geo_reference, |
---|
[3297] | 295 | data_points, |
---|
| 296 | points_are_lats_longs): |
---|
[3745] | 297 | |
---|
[3280] | 298 | if geo_reference is not None: |
---|
[3320] | 299 | msg = """A georeference is specified yet latitude and longitude |
---|
| 300 | are also specified!""" |
---|
[3280] | 301 | raise ValueError, msg |
---|
[3745] | 302 | |
---|
[3297] | 303 | if data_points is not None and not points_are_lats_longs: |
---|
[3320] | 304 | msg = """Data points are specified yet latitude and |
---|
| 305 | longitude are also specified!""" |
---|
[3280] | 306 | raise ValueError, msg |
---|
| 307 | |
---|
[3297] | 308 | if points_are_lats_longs: |
---|
| 309 | if data_points is None: |
---|
| 310 | msg = """Data points are not specified !""" |
---|
| 311 | raise ValueError, msg |
---|
| 312 | lats_longs = ensure_numeric(data_points) |
---|
[3320] | 313 | latitudes = ravel(lats_longs[:,0:1]) |
---|
| 314 | longitudes = ravel(lats_longs[:,1:]) |
---|
[3297] | 315 | |
---|
| 316 | if latitudes is None and longitudes is None: |
---|
| 317 | msg = """Latitudes and Longitudes are not.""" |
---|
| 318 | raise ValueError, msg |
---|
| 319 | |
---|
[3280] | 320 | if latitudes is None: |
---|
| 321 | msg = """Longitudes are specified yet latitudes aren't!""" |
---|
| 322 | raise ValueError, msg |
---|
| 323 | |
---|
| 324 | if longitudes is None: |
---|
[3297] | 325 | msg = """Latitudes are specified yet longitudes aren't!""" |
---|
[3280] | 326 | raise ValueError, msg |
---|
| 327 | |
---|
[3745] | 328 | data_points, zone = convert_from_latlon_to_utm(latitudes=latitudes, |
---|
| 329 | longitudes=longitudes) |
---|
[3280] | 330 | |
---|
| 331 | return data_points, Geo_reference(zone=zone) |
---|
| 332 | |
---|
[2309] | 333 | def get_geo_reference(self): |
---|
| 334 | return self.geo_reference |
---|
[2472] | 335 | |
---|
[3154] | 336 | def get_data_points(self, absolute = True, geo_reference=None): |
---|
[2590] | 337 | """Get coordinates for all data points as an Nx2 array |
---|
| 338 | |
---|
| 339 | If absolute is True absolute UTM coordinates are returned otherwise |
---|
| 340 | returned coordinates are relative to the internal georeference's |
---|
| 341 | xll and yll corners. |
---|
| 342 | |
---|
[3177] | 343 | If a geo_reference is passed the points are returned relative |
---|
| 344 | to that geo_reference. |
---|
| 345 | |
---|
[2590] | 346 | Default: absolute is True. |
---|
| 347 | """ |
---|
[2594] | 348 | |
---|
[3154] | 349 | if absolute is True and geo_reference is None: |
---|
[2594] | 350 | return self.geo_reference.get_absolute(self.data_points) |
---|
[3154] | 351 | elif geo_reference is not None: |
---|
| 352 | return geo_reference.change_points_geo_ref \ |
---|
| 353 | (self.data_points, |
---|
| 354 | self.geo_reference) |
---|
| 355 | else: |
---|
| 356 | return self.data_points |
---|
[2594] | 357 | |
---|
[2472] | 358 | |
---|
[2309] | 359 | def get_attributes(self, attribute_name = None): |
---|
[2303] | 360 | """Return values for one named attribute. |
---|
| 361 | |
---|
| 362 | If attribute_name is None, default_attribute_name is used |
---|
| 363 | """ |
---|
| 364 | |
---|
| 365 | if attribute_name is None: |
---|
| 366 | if self.default_attribute_name is not None: |
---|
| 367 | attribute_name = self.default_attribute_name |
---|
| 368 | else: |
---|
[2465] | 369 | attribute_name = self.attributes.keys()[0] |
---|
| 370 | # above line takes the first one from keys |
---|
[2303] | 371 | |
---|
| 372 | |
---|
| 373 | msg = 'Attribute name %s does not exist in data set' %attribute_name |
---|
| 374 | assert self.attributes.has_key(attribute_name), msg |
---|
| 375 | |
---|
| 376 | return self.attributes[attribute_name] |
---|
[2465] | 377 | |
---|
[2624] | 378 | def get_all_attributes(self): |
---|
| 379 | """ |
---|
| 380 | Return values for all attributes. |
---|
[3738] | 381 | The return value is either None or a dictionary (possibly empty). |
---|
[2624] | 382 | """ |
---|
[2594] | 383 | |
---|
[2624] | 384 | return self.attributes |
---|
| 385 | |
---|
[2592] | 386 | def __add__(self, other): |
---|
[2889] | 387 | """ |
---|
[2928] | 388 | Returns the addition of 2 geospatical objects, |
---|
[2465] | 389 | objects are concatencated to the end of each other |
---|
| 390 | |
---|
| 391 | NOTE: doesn't add if objects contain different |
---|
[2471] | 392 | attributes |
---|
[2624] | 393 | |
---|
| 394 | Always return relative points! |
---|
[2889] | 395 | """ |
---|
[2928] | 396 | |
---|
[2594] | 397 | # find objects zone and checks if the same |
---|
| 398 | geo_ref1 = self.get_geo_reference() |
---|
| 399 | zone1 = geo_ref1.get_zone() |
---|
| 400 | |
---|
| 401 | geo_ref2 = other.get_geo_reference() |
---|
| 402 | zone2 = geo_ref2.get_zone() |
---|
| 403 | |
---|
| 404 | geo_ref1.reconcile_zones(geo_ref2) |
---|
| 405 | |
---|
| 406 | |
---|
[2489] | 407 | # sets xll and yll as the smallest from self and other |
---|
[3723] | 408 | # FIXME (Duncan and Ole): use lower left corner derived from |
---|
| 409 | # absolute coordinates |
---|
[2489] | 410 | if self.geo_reference.xllcorner <= other.geo_reference.xllcorner: |
---|
| 411 | xll = self.geo_reference.xllcorner |
---|
| 412 | else: |
---|
| 413 | xll = other.geo_reference.xllcorner |
---|
| 414 | |
---|
| 415 | if self.geo_reference.yllcorner <= other.geo_reference.yllcorner: |
---|
| 416 | yll = self.geo_reference.yllcorner |
---|
| 417 | else: |
---|
| 418 | yll = other.geo_reference.yllcorner |
---|
[2594] | 419 | new_geo_ref = Geo_reference(geo_ref1.get_zone(), xll, yll) |
---|
[2489] | 420 | |
---|
[2594] | 421 | xll = yll = 0. |
---|
[2624] | 422 | |
---|
| 423 | relative_points1 = self.get_data_points(absolute = False) |
---|
| 424 | relative_points2 = other.get_data_points(absolute = False) |
---|
[2594] | 425 | |
---|
| 426 | |
---|
[3723] | 427 | new_relative_points1 = new_geo_ref.\ |
---|
| 428 | change_points_geo_ref(relative_points1, |
---|
| 429 | geo_ref1) |
---|
| 430 | new_relative_points2 = new_geo_ref.\ |
---|
| 431 | change_points_geo_ref(relative_points2, |
---|
| 432 | geo_ref2) |
---|
[2594] | 433 | |
---|
[3723] | 434 | # Now both point sets are relative to new_geo_ref and |
---|
| 435 | # zones have been reconciled |
---|
[2594] | 436 | |
---|
| 437 | # Concatenate points |
---|
| 438 | new_points = concatenate((new_relative_points1, |
---|
| 439 | new_relative_points2), |
---|
[2624] | 440 | axis = 0) |
---|
[2594] | 441 | |
---|
[3723] | 442 | # Concatenate attributes if any |
---|
| 443 | if self.attributes is None: |
---|
| 444 | if other.attributes is not None: |
---|
[2594] | 445 | msg = 'Both geospatial_data objects must have the same \n' |
---|
| 446 | msg += 'attributes to allow addition.' |
---|
[3723] | 447 | raise Exception, msg |
---|
| 448 | |
---|
| 449 | new_attributes = None |
---|
| 450 | else: |
---|
| 451 | new_attributes = {} |
---|
| 452 | for x in self.attributes.keys(): |
---|
| 453 | if other.attributes.has_key(x): |
---|
[2594] | 454 | |
---|
[3723] | 455 | attrib1 = self.attributes[x] |
---|
| 456 | attrib2 = other.attributes[x] |
---|
| 457 | new_attributes[x] = concatenate((attrib1, attrib2)) |
---|
| 458 | |
---|
| 459 | else: |
---|
| 460 | msg = 'Both geospatial_data objects must have the same \n' |
---|
| 461 | msg += 'attributes to allow addition.' |
---|
| 462 | raise Exception, msg |
---|
| 463 | |
---|
[2594] | 464 | # Instantiate new data object and return |
---|
| 465 | return Geospatial_data(new_points, |
---|
| 466 | new_attributes, |
---|
[2698] | 467 | new_geo_ref) |
---|
[2570] | 468 | |
---|
[2624] | 469 | ### |
---|
| 470 | # IMPORT/EXPORT POINTS FILES |
---|
| 471 | ### |
---|
[2928] | 472 | |
---|
[2940] | 473 | def import_points_file(self, file_name, delimiter = None, verbose = False): |
---|
[2624] | 474 | """ load an .xya or .pts file |
---|
| 475 | Note: will throw an IOError if it can't load the file. |
---|
| 476 | Catch these! |
---|
[3738] | 477 | |
---|
| 478 | Post condition: self.attributes dictionary has been set |
---|
[2624] | 479 | """ |
---|
| 480 | |
---|
[2940] | 481 | if access(file_name, F_OK) == 0 : |
---|
| 482 | msg = 'File %s does not exist or is not accessible' %file_name |
---|
| 483 | raise IOError, msg |
---|
| 484 | |
---|
[2643] | 485 | attributes = {} |
---|
[2940] | 486 | if file_name[-4:]== ".xya": |
---|
[2624] | 487 | try: |
---|
| 488 | if delimiter == None: |
---|
| 489 | try: |
---|
[2940] | 490 | fd = open(file_name) |
---|
[2928] | 491 | data_points, attributes, geo_reference = _read_xya_file(fd, ',') |
---|
[2941] | 492 | except TitleError: |
---|
[2624] | 493 | fd.close() |
---|
[2940] | 494 | fd = open(file_name) |
---|
[2928] | 495 | data_points, attributes, geo_reference = _read_xya_file(fd, ' ') |
---|
[2624] | 496 | else: |
---|
[2940] | 497 | fd = open(file_name) |
---|
[2928] | 498 | data_points, attributes, geo_reference = _read_xya_file(fd, delimiter) |
---|
[2624] | 499 | fd.close() |
---|
| 500 | except (IndexError,ValueError,SyntaxError): |
---|
| 501 | fd.close() |
---|
[2940] | 502 | msg = 'Could not open file %s ' %file_name |
---|
[2624] | 503 | raise IOError, msg |
---|
| 504 | except IOError, e: |
---|
[2955] | 505 | fd.close() |
---|
[2624] | 506 | # Catch this to add an error message |
---|
[2940] | 507 | msg = 'Could not open file or incorrect file format %s:%s' %(file_name, e) |
---|
[2624] | 508 | raise IOError, msg |
---|
| 509 | |
---|
[2940] | 510 | elif file_name[-4:]== ".pts": |
---|
[2624] | 511 | try: |
---|
[2940] | 512 | data_points, attributes, geo_reference = _read_pts_file(file_name, verbose) |
---|
[2624] | 513 | except IOError, e: |
---|
[2940] | 514 | msg = 'Could not open file %s ' %file_name |
---|
[2624] | 515 | raise IOError, msg |
---|
| 516 | else: |
---|
[2940] | 517 | msg = 'Extension %s is unknown' %file_name[-4:] |
---|
[2624] | 518 | raise IOError, msg |
---|
[2643] | 519 | |
---|
| 520 | # print'in import data_points', data_points |
---|
| 521 | # print'in import attributes', attributes |
---|
| 522 | # print'in import data_points', geo_reference |
---|
| 523 | self.data_points = data_points |
---|
| 524 | self.attributes = attributes |
---|
| 525 | self.geo_reference = geo_reference |
---|
[2624] | 526 | |
---|
[2643] | 527 | # return all_data |
---|
[2624] | 528 | |
---|
[2940] | 529 | def export_points_file(self, file_name, absolute=True): |
---|
[2698] | 530 | |
---|
[2624] | 531 | """ |
---|
[2940] | 532 | write a points file, file_name, as a text (.xya) or binary (.pts) file |
---|
| 533 | file_name is the file name, including the extension |
---|
[2624] | 534 | The point_dict is defined at the top of this file. |
---|
[2772] | 535 | |
---|
| 536 | If absolute is True data points at returned added to the xll and yll |
---|
| 537 | and geo_reference as None |
---|
| 538 | |
---|
| 539 | If absolute is False data points at returned as relative to the xll |
---|
| 540 | and yll and geo_reference remains uneffected |
---|
[2624] | 541 | """ |
---|
[2698] | 542 | |
---|
[2940] | 543 | if (file_name[-4:] == ".xya"): |
---|
[2698] | 544 | if absolute is True: |
---|
[3735] | 545 | _write_xya_file(file_name, |
---|
| 546 | self.get_data_points(absolute=True), |
---|
| 547 | self.get_all_attributes()) |
---|
[2698] | 548 | else: |
---|
[3735] | 549 | _write_xya_file(file_name, |
---|
| 550 | self.get_data_points(absolute=False), |
---|
| 551 | self.get_all_attributes(), |
---|
| 552 | self.get_geo_reference()) |
---|
[2698] | 553 | |
---|
[2940] | 554 | elif (file_name[-4:] == ".pts"): |
---|
[2698] | 555 | if absolute is True: |
---|
[3735] | 556 | _write_pts_file(file_name, |
---|
| 557 | self.get_data_points(absolute), |
---|
| 558 | self.get_all_attributes()) |
---|
[2698] | 559 | else: |
---|
[3735] | 560 | _write_pts_file(file_name, |
---|
| 561 | self.get_data_points(absolute), |
---|
| 562 | self.get_all_attributes(), |
---|
| 563 | self.get_geo_reference()) |
---|
[2624] | 564 | else: |
---|
[2940] | 565 | msg = 'Unknown file type %s ' %file_name |
---|
[2624] | 566 | raise IOError, msg |
---|
| 567 | |
---|
[2465] | 568 | |
---|
[2570] | 569 | |
---|
| 570 | |
---|
[2928] | 571 | def _read_pts_file(file_name, verbose = False): |
---|
[2643] | 572 | """Read .pts NetCDF file |
---|
| 573 | |
---|
| 574 | Return a dic of array of points, and dic of array of attribute |
---|
| 575 | eg |
---|
| 576 | dic['points'] = [[1.0,2.0],[3.0,5.0]] |
---|
| 577 | dic['attributelist']['elevation'] = [[7.0,5.0] |
---|
| 578 | """ |
---|
| 579 | |
---|
| 580 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 581 | |
---|
| 582 | if verbose: print 'Reading ', file_name |
---|
| 583 | |
---|
| 584 | |
---|
| 585 | # see if the file is there. Throw a QUIET IO error if it isn't |
---|
| 586 | fd = open(file_name,'r') |
---|
| 587 | fd.close() |
---|
| 588 | |
---|
| 589 | #throws prints to screen if file not present |
---|
| 590 | fid = NetCDFFile(file_name, 'r') |
---|
| 591 | |
---|
| 592 | # point_atts = {} |
---|
| 593 | # Get the variables |
---|
| 594 | # point_atts['pointlist'] = array(fid.variables['points']) |
---|
| 595 | pointlist = array(fid.variables['points']) |
---|
| 596 | keys = fid.variables.keys() |
---|
| 597 | if verbose: print 'Got %d variables: %s' %(len(keys), keys) |
---|
| 598 | try: |
---|
| 599 | keys.remove('points') |
---|
| 600 | except IOError, e: |
---|
| 601 | fid.close() |
---|
| 602 | msg = 'Expected keyword "points" but could not find it' |
---|
| 603 | raise IOError, msg |
---|
| 604 | |
---|
| 605 | attributes = {} |
---|
| 606 | for key in keys: |
---|
| 607 | if verbose: print "reading attribute '%s'" %key |
---|
| 608 | |
---|
| 609 | attributes[key] = array(fid.variables[key]) |
---|
| 610 | |
---|
| 611 | # point_atts['attributelist'] = attributes |
---|
| 612 | |
---|
| 613 | try: |
---|
| 614 | geo_reference = Geo_reference(NetCDFObject=fid) |
---|
| 615 | # point_atts['geo_reference'] = geo_reference |
---|
| 616 | except AttributeError, e: |
---|
| 617 | #geo_ref not compulsory |
---|
| 618 | # point_atts['geo_reference'] = None |
---|
| 619 | geo_reference = None |
---|
| 620 | |
---|
| 621 | fid.close() |
---|
| 622 | |
---|
| 623 | return pointlist, attributes, geo_reference |
---|
| 624 | |
---|
| 625 | |
---|
[2928] | 626 | def _read_xya_file( fd, delimiter): |
---|
[2643] | 627 | points = [] |
---|
| 628 | pointattributes = [] |
---|
| 629 | title = fd.readline() |
---|
| 630 | att_names = clean_line(title,delimiter) |
---|
| 631 | att_dict = {} |
---|
| 632 | line = fd.readline() |
---|
| 633 | numbers = clean_line(line,delimiter) |
---|
[2942] | 634 | |
---|
| 635 | while len(numbers) > 1 and line[0] <> '#': |
---|
[2643] | 636 | if numbers != []: |
---|
| 637 | try: |
---|
| 638 | x = float(numbers[0]) |
---|
| 639 | y = float(numbers[1]) |
---|
| 640 | points.append([x,y]) |
---|
| 641 | numbers.pop(0) |
---|
| 642 | numbers.pop(0) |
---|
| 643 | if len(att_names) != len(numbers): |
---|
| 644 | fd.close() |
---|
| 645 | # It might not be a problem with the title |
---|
| 646 | #raise TitleAmountError |
---|
| 647 | raise IOError |
---|
| 648 | for i,num in enumerate(numbers): |
---|
| 649 | num.strip() |
---|
| 650 | if num != '\n' and num != '': |
---|
| 651 | #attributes.append(float(num)) |
---|
| 652 | att_dict.setdefault(att_names[i],[]).append(float(num)) |
---|
| 653 | except ValueError: |
---|
| 654 | raise SyntaxError |
---|
| 655 | line = fd.readline() |
---|
| 656 | numbers = clean_line(line,delimiter) |
---|
[2942] | 657 | |
---|
[2643] | 658 | if line == '': |
---|
| 659 | geo_reference = None |
---|
| 660 | else: |
---|
| 661 | geo_reference = Geo_reference(ASCIIFile=fd,read_title=line) |
---|
[2940] | 662 | |
---|
[2643] | 663 | |
---|
| 664 | pointlist = array(points).astype(Float) |
---|
| 665 | for key in att_dict.keys(): |
---|
| 666 | att_dict[key] = array(att_dict[key]).astype(Float) |
---|
| 667 | |
---|
| 668 | return pointlist, att_dict, geo_reference |
---|
| 669 | |
---|
[3735] | 670 | def _write_pts_file(file_name, |
---|
| 671 | write_data_points, |
---|
| 672 | write_attributes=None, |
---|
| 673 | write_geo_reference=None): |
---|
[2928] | 674 | """ |
---|
| 675 | Write .pts NetCDF file |
---|
[2570] | 676 | |
---|
[2928] | 677 | NOTE: Below might not be valid ask Duncan : NB 5/2006 |
---|
| 678 | |
---|
[2624] | 679 | WARNING: This function mangles the point_atts data structure |
---|
[2928] | 680 | #F??ME: (DSG)This format has issues. |
---|
[2624] | 681 | # There can't be an attribute called points |
---|
| 682 | # consider format change |
---|
| 683 | # method changed by NB not sure if above statement is correct |
---|
[2928] | 684 | |
---|
| 685 | should create new test for this |
---|
[2624] | 686 | legal_keys = ['pointlist', 'attributelist', 'geo_reference'] |
---|
| 687 | for key in point_atts.keys(): |
---|
| 688 | msg = 'Key %s is illegal. Valid keys are %s' %(key, legal_keys) |
---|
| 689 | assert key in legal_keys, msg |
---|
[2928] | 690 | """ |
---|
[2624] | 691 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 692 | # NetCDF file definition |
---|
| 693 | outfile = NetCDFFile(file_name, 'w') |
---|
| 694 | |
---|
| 695 | #Create new file |
---|
| 696 | outfile.institution = 'Geoscience Australia' |
---|
| 697 | outfile.description = 'NetCDF format for compact and portable storage ' +\ |
---|
[3735] | 698 | 'of spatial point data' |
---|
[2624] | 699 | |
---|
| 700 | # dimension definitions |
---|
| 701 | shape = write_data_points.shape[0] |
---|
| 702 | outfile.createDimension('number_of_points', shape) |
---|
| 703 | outfile.createDimension('number_of_dimensions', 2) #This is 2d data |
---|
| 704 | |
---|
| 705 | # variable definition |
---|
| 706 | outfile.createVariable('points', Float, ('number_of_points', |
---|
| 707 | 'number_of_dimensions')) |
---|
| 708 | |
---|
| 709 | #create variables |
---|
| 710 | outfile.variables['points'][:] = write_data_points #.astype(Float32) |
---|
| 711 | |
---|
[3735] | 712 | if write_attributes is not None: |
---|
| 713 | for key in write_attributes.keys(): |
---|
| 714 | outfile.createVariable(key, Float, ('number_of_points',)) |
---|
| 715 | outfile.variables[key][:] = write_attributes[key] #.astype(Float32) |
---|
[2624] | 716 | |
---|
| 717 | if write_geo_reference is not None: |
---|
| 718 | write_geo_reference.write_NetCDF(outfile) |
---|
| 719 | |
---|
| 720 | outfile.close() |
---|
| 721 | |
---|
| 722 | |
---|
| 723 | |
---|
[3735] | 724 | def _write_xya_file(file_name, |
---|
| 725 | write_data_points, |
---|
| 726 | write_attributes=None, |
---|
| 727 | write_geo_reference=None, |
---|
| 728 | delimiter = ','): |
---|
[2624] | 729 | """ |
---|
[2940] | 730 | export a file, file_name, with the xya format |
---|
[2624] | 731 | |
---|
| 732 | """ |
---|
| 733 | points = write_data_points |
---|
| 734 | pointattributes = write_attributes |
---|
| 735 | |
---|
| 736 | fd = open(file_name,'w') |
---|
| 737 | titlelist = "" |
---|
[3735] | 738 | if pointattributes is not None: |
---|
| 739 | for title in pointattributes.keys(): |
---|
| 740 | titlelist = titlelist + title + delimiter |
---|
| 741 | titlelist = titlelist[0:-len(delimiter)] # remove the last delimiter |
---|
[2624] | 742 | fd.write(titlelist+"\n") |
---|
[3735] | 743 | |
---|
[2624] | 744 | #<vertex #> <x> <y> [attributes] |
---|
[3735] | 745 | for i, vert in enumerate( points): |
---|
[2928] | 746 | |
---|
[3735] | 747 | |
---|
| 748 | if pointattributes is not None: |
---|
| 749 | attlist = "," |
---|
| 750 | for att in pointattributes.keys(): |
---|
| 751 | attlist = attlist + str(pointattributes[att][i])+ delimiter |
---|
| 752 | attlist = attlist[0:-len(delimiter)] # remove the last delimiter |
---|
| 753 | attlist.strip() |
---|
| 754 | else: |
---|
| 755 | attlist = '' |
---|
| 756 | |
---|
| 757 | fd.write(str(vert[0]) + delimiter + |
---|
| 758 | str(vert[1]) + attlist + "\n") |
---|
| 759 | |
---|
[2624] | 760 | if write_geo_reference is not None: |
---|
| 761 | write_geo_reference.write_ASCII(fd) |
---|
| 762 | fd.close() |
---|
[2570] | 763 | |
---|
| 764 | |
---|
| 765 | |
---|
| 766 | def _point_atts2array(point_atts): |
---|
| 767 | point_atts['pointlist'] = array(point_atts['pointlist']).astype(Float) |
---|
| 768 | |
---|
| 769 | for key in point_atts['attributelist'].keys(): |
---|
| 770 | point_atts['attributelist'][key]= array(point_atts['attributelist'][key]).astype(Float) |
---|
| 771 | return point_atts |
---|
| 772 | |
---|
| 773 | |
---|
| 774 | |
---|
| 775 | |
---|
[2306] | 776 | def geospatial_data2points_dictionary(geospatial_data): |
---|
[2303] | 777 | """Convert geospatial data to points_dictionary |
---|
| 778 | """ |
---|
| 779 | |
---|
[2306] | 780 | points_dictionary = {} |
---|
| 781 | points_dictionary['pointlist'] = geospatial_data.data_points |
---|
| 782 | |
---|
| 783 | points_dictionary['attributelist'] = {} |
---|
| 784 | |
---|
| 785 | for attribute_name in geospatial_data.attributes.keys(): |
---|
| 786 | val = geospatial_data.attributes[attribute_name] |
---|
| 787 | points_dictionary['attributelist'][attribute_name] = val |
---|
| 788 | |
---|
| 789 | points_dictionary['geo_reference'] = geospatial_data.geo_reference |
---|
| 790 | |
---|
| 791 | return points_dictionary |
---|
| 792 | |
---|
| 793 | |
---|
[2303] | 794 | def points_dictionary2geospatial_data(points_dictionary): |
---|
| 795 | """Convert points_dictionary to geospatial data object |
---|
| 796 | """ |
---|
| 797 | |
---|
[2306] | 798 | msg = 'Points dictionary must have key pointlist' |
---|
| 799 | assert points_dictionary.has_key('pointlist'), msg |
---|
[2303] | 800 | |
---|
[2306] | 801 | msg = 'Points dictionary must have key attributelist' |
---|
| 802 | assert points_dictionary.has_key('attributelist'), msg |
---|
| 803 | |
---|
| 804 | if points_dictionary.has_key('geo_reference'): |
---|
| 805 | geo = points_dictionary['geo_reference'] |
---|
| 806 | else: |
---|
| 807 | geo = None |
---|
| 808 | |
---|
| 809 | return Geospatial_data(points_dictionary['pointlist'], |
---|
| 810 | points_dictionary['attributelist'], |
---|
| 811 | geo_reference = geo) |
---|
| 812 | |
---|
[2570] | 813 | def clean_line(line,delimiter): |
---|
| 814 | """Remove whitespace |
---|
| 815 | """ |
---|
| 816 | #print ">%s" %line |
---|
| 817 | line = line.strip() |
---|
| 818 | #print "stripped>%s" %line |
---|
| 819 | numbers = line.split(delimiter) |
---|
| 820 | i = len(numbers) - 1 |
---|
| 821 | while i >= 0: |
---|
| 822 | if numbers[i] == '': |
---|
| 823 | numbers.pop(i) |
---|
| 824 | else: |
---|
| 825 | numbers[i] = numbers[i].strip() |
---|
| 826 | |
---|
| 827 | i += -1 |
---|
| 828 | #for num in numbers: |
---|
| 829 | # print "num>%s<" %num |
---|
| 830 | return numbers |
---|
[2303] | 831 | |
---|
[2940] | 832 | def xxx_add_points_files(add_file1, add_file2, results_file): |
---|
[2928] | 833 | """ adds the points and attruibutes of 2 pts or xya files and |
---|
[2584] | 834 | writes it to a pts file |
---|
| 835 | |
---|
| 836 | NOTE will add the function to check and remove points from one set |
---|
| 837 | that are shared. This will require some work and maybe __subtract__ function |
---|
[2928] | 838 | """ |
---|
[2591] | 839 | |
---|
[2584] | 840 | G1 = Geospatial_data(file_name = add_file1) |
---|
| 841 | G2 = Geospatial_data(file_name = add_file2) |
---|
[2590] | 842 | new_add_file2 = add_file2[:-4] + '.pts' |
---|
[2584] | 843 | |
---|
[2590] | 844 | G = G1 + G2 |
---|
[2584] | 845 | |
---|
| 846 | #FIXME remove dependance on points to dict in export only! |
---|
[2624] | 847 | # G_points_dict = geospatial_data2points_dictionary(G) |
---|
| 848 | # export_points_file(results_file, G_points_dict) |
---|
| 849 | |
---|
| 850 | # G_points_dict = geospatial_data2points_dictionary(G) |
---|
| 851 | |
---|
[2698] | 852 | G.export_points_file(results_file) |
---|
[2889] | 853 | |
---|
[2890] | 854 | # ' |
---|
| 855 | def ensure_absolute(points, geo_reference = None): |
---|
[2889] | 856 | """ |
---|
| 857 | This function inputs several formats and |
---|
[2890] | 858 | outputs one format. - a numeric array of absolute points. |
---|
| 859 | |
---|
| 860 | Inputed formats are; |
---|
| 861 | points: List or numeric array of coordinate pairs [xi, eta] of |
---|
[3266] | 862 | points or geospatial object or points file name |
---|
[2890] | 863 | |
---|
| 864 | mesh_origin: A geo_reference object or 3-tuples consisting of |
---|
| 865 | UTM zone, easting and northing. |
---|
| 866 | If specified vertex coordinates are assumed to be |
---|
| 867 | relative to their respective origins. |
---|
[2889] | 868 | """ |
---|
[3262] | 869 | if isinstance(points,type('')): |
---|
| 870 | #It's a string |
---|
| 871 | #assume it is a point file |
---|
| 872 | points = Geospatial_data(file_name = points) |
---|
| 873 | |
---|
[2890] | 874 | if isinstance(points,Geospatial_data): |
---|
| 875 | points = points.get_data_points( \ |
---|
| 876 | absolute = True) |
---|
| 877 | msg = "Use a Geospatial_data object or a mesh origin. Not both." |
---|
| 878 | assert geo_reference == None, msg |
---|
| 879 | |
---|
| 880 | else: |
---|
| 881 | points = ensure_numeric(points, Float) |
---|
| 882 | if geo_reference is None: |
---|
| 883 | geo = None #Geo_reference() |
---|
| 884 | else: |
---|
| 885 | if isinstance(geo_reference, Geo_reference): |
---|
| 886 | geo = geo_reference |
---|
| 887 | else: |
---|
| 888 | geo = Geo_reference(geo_reference[0], |
---|
| 889 | geo_reference[1], |
---|
| 890 | geo_reference[2]) |
---|
| 891 | points = geo.get_absolute(points) |
---|
| 892 | return points |
---|
[2584] | 893 | |
---|
[3056] | 894 | |
---|
| 895 | def ensure_geospatial(points, geo_reference = None): |
---|
| 896 | """ |
---|
| 897 | This function inputs several formats and |
---|
| 898 | outputs one format. - a geospatial_data instance. |
---|
| 899 | |
---|
| 900 | Inputed formats are; |
---|
| 901 | points: List or numeric array of coordinate pairs [xi, eta] of |
---|
| 902 | points or geospatial object |
---|
| 903 | |
---|
| 904 | mesh_origin: A geo_reference object or 3-tuples consisting of |
---|
| 905 | UTM zone, easting and northing. |
---|
| 906 | If specified vertex coordinates are assumed to be |
---|
| 907 | relative to their respective origins. |
---|
| 908 | """ |
---|
| 909 | if isinstance(points,Geospatial_data): |
---|
| 910 | msg = "Use a Geospatial_data object or a mesh origin. Not both." |
---|
| 911 | assert geo_reference == None, msg |
---|
[3351] | 912 | return points |
---|
[3056] | 913 | else: |
---|
| 914 | points = ensure_numeric(points, Float) |
---|
| 915 | if geo_reference is None: |
---|
[3262] | 916 | geo = None |
---|
[3056] | 917 | else: |
---|
| 918 | if isinstance(geo_reference, Geo_reference): |
---|
| 919 | geo = geo_reference |
---|
| 920 | else: |
---|
| 921 | geo = Geo_reference(geo_reference[0], |
---|
| 922 | geo_reference[1], |
---|
| 923 | geo_reference[2]) |
---|
[3351] | 924 | points = Geospatial_data(data_points=points, geo_reference=geo) |
---|
[3056] | 925 | return points |
---|
| 926 | |
---|