[2488] | 1 | """ |
---|
[4451] | 2 | |
---|
[2488] | 3 | """ |
---|
| 4 | |
---|
| 5 | |
---|
| 6 | #FIXME: Ensure that all attributes of a georef are treated everywhere |
---|
| 7 | #and unit test |
---|
| 8 | |
---|
[6438] | 9 | |
---|
[2593] | 10 | import types, sys |
---|
[6438] | 11 | import copy |
---|
[3514] | 12 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
| 13 | from anuga.utilities.anuga_exceptions import ANUGAError, TitleError, ParsingError, \ |
---|
[3134] | 14 | ShapeError |
---|
[2488] | 15 | |
---|
[6149] | 16 | import Numeric as num |
---|
[2941] | 17 | |
---|
[6149] | 18 | |
---|
[2602] | 19 | DEFAULT_ZONE = -1 |
---|
[4663] | 20 | TITLE = '#geo reference' + "\n" #this title is referred to in the test format |
---|
[2488] | 21 | |
---|
[4387] | 22 | DEFAULT_PROJECTION = 'UTM' |
---|
| 23 | DEFAULT_DATUM = 'wgs84' |
---|
| 24 | DEFAULT_UNITS = 'm' |
---|
| 25 | DEFAULT_FALSE_EASTING = 500000 |
---|
| 26 | DEFAULT_FALSE_NORTHING = 10000000 #Default for southern hemisphere |
---|
| 27 | |
---|
[2488] | 28 | class Geo_reference: |
---|
| 29 | """ |
---|
| 30 | """ |
---|
| 31 | |
---|
| 32 | def __init__(self, |
---|
| 33 | zone = DEFAULT_ZONE, |
---|
| 34 | xllcorner = 0.0, |
---|
| 35 | yllcorner = 0.0, |
---|
[4387] | 36 | datum = DEFAULT_DATUM, |
---|
| 37 | projection = DEFAULT_PROJECTION, |
---|
| 38 | units = DEFAULT_UNITS, |
---|
| 39 | false_easting = DEFAULT_FALSE_EASTING, |
---|
| 40 | false_northing = DEFAULT_FALSE_NORTHING, |
---|
[2488] | 41 | NetCDFObject=None, |
---|
| 42 | ASCIIFile=None, |
---|
| 43 | read_title=None): |
---|
| 44 | """ |
---|
| 45 | input: |
---|
| 46 | NetCDFObject - a handle to the netCDF file to be written to |
---|
| 47 | ASCIIFile - a handle to the text file |
---|
| 48 | read_title - the title of the georeference text, if it was read in. |
---|
[2941] | 49 | If the function that calls this has already read the title line, |
---|
| 50 | it can't unread it, so this info has to be passed. |
---|
| 51 | If you know of a way to unread this info, then tell us. |
---|
[2488] | 52 | |
---|
[4387] | 53 | Note, the text file only saves a sub set of the info the |
---|
| 54 | points file does. Currently the info not written in text |
---|
| 55 | must be the default info, since ANUGA assumes it isn't |
---|
| 56 | changing. |
---|
| 57 | """ |
---|
[4451] | 58 | if zone is None: |
---|
| 59 | zone = DEFAULT_ZONE |
---|
[2488] | 60 | self.false_easting = false_easting |
---|
| 61 | self.false_northing = false_northing |
---|
| 62 | self.datum = datum |
---|
| 63 | self.projection = projection |
---|
| 64 | self.zone = zone |
---|
| 65 | self.units = units |
---|
| 66 | self.xllcorner = xllcorner |
---|
| 67 | self.yllcorner = yllcorner |
---|
| 68 | |
---|
| 69 | if NetCDFObject is not None: |
---|
| 70 | self.read_NetCDF(NetCDFObject) |
---|
| 71 | |
---|
| 72 | if ASCIIFile is not None: |
---|
| 73 | self.read_ASCII(ASCIIFile, read_title=read_title) |
---|
| 74 | |
---|
| 75 | def get_xllcorner(self): |
---|
| 76 | return self.xllcorner |
---|
| 77 | |
---|
| 78 | def get_yllcorner(self): |
---|
| 79 | return self.yllcorner |
---|
| 80 | |
---|
| 81 | def get_zone(self): |
---|
| 82 | return self.zone |
---|
| 83 | |
---|
| 84 | def write_NetCDF(self, outfile): |
---|
| 85 | outfile.xllcorner = self.xllcorner |
---|
| 86 | outfile.yllcorner = self.yllcorner |
---|
| 87 | outfile.zone = self.zone |
---|
| 88 | |
---|
| 89 | outfile.false_easting = self.false_easting |
---|
| 90 | outfile.false_northing = self.false_northing |
---|
| 91 | |
---|
| 92 | outfile.datum = self.datum |
---|
| 93 | outfile.projection = self.projection |
---|
| 94 | outfile.units = self.units |
---|
| 95 | |
---|
| 96 | def read_NetCDF(self, infile): |
---|
| 97 | self.xllcorner = infile.xllcorner[0] |
---|
| 98 | self.yllcorner = infile.yllcorner[0] |
---|
| 99 | self.zone = infile.zone[0] |
---|
| 100 | |
---|
[4387] | 101 | |
---|
[2488] | 102 | # Fix some assertion failures |
---|
[6149] | 103 | if type(self.zone) == num.ArrayType and self.zone.shape == (): |
---|
[2488] | 104 | self.zone = self.zone[0] |
---|
[6149] | 105 | if type(self.xllcorner) == num.ArrayType and self.xllcorner.shape == (): |
---|
[2488] | 106 | self.xllcorner = self.xllcorner[0] |
---|
[6149] | 107 | if type(self.yllcorner) == num.ArrayType and self.yllcorner.shape == (): |
---|
[2488] | 108 | self.yllcorner = self.yllcorner[0] |
---|
| 109 | |
---|
| 110 | assert (type(self.xllcorner) == types.FloatType or\ |
---|
| 111 | type(self.xllcorner) == types.IntType) |
---|
| 112 | assert (type(self.yllcorner) == types.FloatType or\ |
---|
| 113 | type(self.yllcorner) == types.IntType) |
---|
| 114 | assert (type(self.zone) == types.IntType) |
---|
[4387] | 115 | |
---|
[4389] | 116 | try: |
---|
| 117 | self.false_easting = infile.false_easting[0] |
---|
| 118 | self.false_northing = infile.false_northing[0] |
---|
| 119 | |
---|
| 120 | self.datum = infile.datum |
---|
| 121 | self.projection = infile.projection |
---|
| 122 | self.units = infile.units |
---|
| 123 | except: |
---|
| 124 | pass |
---|
[4451] | 125 | if (self.false_easting != DEFAULT_FALSE_EASTING): |
---|
| 126 | print "WARNING: False easting of %f specified." %self.false_easting |
---|
| 127 | print "Default false easting is %f." %DEFAULT_FALSE_EASTING |
---|
| 128 | print "ANUGA does not correct for differences in False Eastings." |
---|
| 129 | |
---|
| 130 | if (self.false_northing != DEFAULT_FALSE_NORTHING): |
---|
| 131 | print "WARNING: False northing of %f specified." \ |
---|
| 132 | %self.false_northing |
---|
| 133 | print "Default false northing is %f." %DEFAULT_FALSE_NORTHING |
---|
| 134 | print "ANUGA does not correct for differences in False Northings." |
---|
| 135 | |
---|
| 136 | if (self.datum.upper() != DEFAULT_DATUM.upper()): |
---|
| 137 | print "WARNING: Datum of %s specified." \ |
---|
| 138 | %self.datum |
---|
| 139 | print "Default Datum is %s." %DEFAULT_DATUM |
---|
| 140 | print "ANUGA does not correct for differences in datums." |
---|
| 141 | |
---|
| 142 | if (self.projection.upper() != DEFAULT_PROJECTION.upper()): |
---|
| 143 | print "WARNING: Projection of %s specified." \ |
---|
| 144 | %self.projection |
---|
| 145 | print "Default Projection is %s." %DEFAULT_PROJECTION |
---|
| 146 | print "ANUGA does not correct for differences in Projection." |
---|
| 147 | |
---|
| 148 | if (self.units.upper() != DEFAULT_UNITS.upper()): |
---|
| 149 | print "WARNING: Units of %s specified." \ |
---|
| 150 | %self.units |
---|
| 151 | print "Default units is %s." %DEFAULT_UNITS |
---|
| 152 | print "ANUGA does not correct for differences in units." |
---|
[2488] | 153 | |
---|
[4663] | 154 | ### ASCII files with geo-refs are currently not used |
---|
[2488] | 155 | def write_ASCII(self, fd): |
---|
| 156 | fd.write(TITLE) |
---|
| 157 | fd.write(str(self.zone) + "\n") |
---|
| 158 | fd.write(str(self.xllcorner) + "\n") |
---|
| 159 | fd.write(str(self.yllcorner) + "\n") |
---|
| 160 | |
---|
| 161 | def read_ASCII(self, fd,read_title=None): |
---|
[2942] | 162 | try: |
---|
| 163 | if read_title == None: |
---|
| 164 | read_title = fd.readline() # remove the title line |
---|
| 165 | if read_title[0:2].upper() != TITLE[0:2].upper(): |
---|
| 166 | msg = 'File error. Expecting line: %s. Got this line: %s' \ |
---|
| 167 | %(TITLE, read_title) |
---|
| 168 | raise TitleError, msg |
---|
| 169 | self.zone = int(fd.readline()) |
---|
| 170 | self.xllcorner = float(fd.readline()) |
---|
| 171 | self.yllcorner = float(fd.readline()) |
---|
| 172 | except SyntaxError: |
---|
| 173 | msg = 'File error. Got syntax error while parsing geo reference' |
---|
| 174 | raise ParsingError, msg |
---|
[2488] | 175 | |
---|
| 176 | # Fix some assertion failures |
---|
[6149] | 177 | if(type(self.zone) == num.ArrayType and self.zone.shape == ()): |
---|
[2488] | 178 | self.zone = self.zone[0] |
---|
[6149] | 179 | if type(self.xllcorner) == num.ArrayType and self.xllcorner.shape == (): |
---|
[2488] | 180 | self.xllcorner = self.xllcorner[0] |
---|
[6149] | 181 | if type(self.yllcorner) == num.ArrayType and self.yllcorner.shape == (): |
---|
[2488] | 182 | self.yllcorner = self.yllcorner[0] |
---|
| 183 | |
---|
| 184 | assert (type(self.xllcorner) == types.FloatType) |
---|
| 185 | assert (type(self.yllcorner) == types.FloatType) |
---|
| 186 | assert (type(self.zone) == types.IntType) |
---|
| 187 | |
---|
| 188 | |
---|
| 189 | def change_points_geo_ref(self, points, |
---|
| 190 | points_geo_ref=None): |
---|
| 191 | """ |
---|
| 192 | Change the geo reference of a list or Numeric array of points to |
---|
| 193 | be this reference.(The reference used for this object) |
---|
| 194 | If the points do not have a geo ref, assume 'absolute' values |
---|
| 195 | """ |
---|
[6438] | 196 | import copy |
---|
| 197 | |
---|
[2488] | 198 | is_list = False |
---|
| 199 | if type(points) == types.ListType: |
---|
| 200 | is_list = True |
---|
[2594] | 201 | |
---|
[6438] | 202 | points = ensure_numeric(copy.copy(points), num.Float) |
---|
[2595] | 203 | |
---|
[2594] | 204 | if len(points.shape) == 1: |
---|
| 205 | #One point has been passed |
---|
| 206 | msg = 'Single point must have two elements' |
---|
| 207 | assert len(points) == 2, msg |
---|
[6149] | 208 | points = num.reshape(points, (1,2)) |
---|
[2594] | 209 | |
---|
[5665] | 210 | msg = 'Points array must be two dimensional.\n' |
---|
| 211 | msg += 'I got %d dimensions' %len(points.shape) |
---|
| 212 | assert len(points.shape) == 2, msg |
---|
| 213 | |
---|
[3071] | 214 | msg = 'Input must be an N x 2 array or list of (x,y) values. ' |
---|
| 215 | msg += 'I got an %d x %d array' %points.shape |
---|
| 216 | assert points.shape[1] == 2, msg |
---|
| 217 | |
---|
[5736] | 218 | |
---|
| 219 | # FIXME (Ole): Could also check if zone, xllcorner, yllcorner |
---|
| 220 | # are identical in the two geo refs. |
---|
[2594] | 221 | if points_geo_ref is not self: |
---|
[5736] | 222 | # If georeferences are different |
---|
| 223 | |
---|
[2594] | 224 | if not points_geo_ref is None: |
---|
[5736] | 225 | # Convert points to absolute coordinates |
---|
[2594] | 226 | points[:,0] += points_geo_ref.xllcorner |
---|
| 227 | points[:,1] += points_geo_ref.yllcorner |
---|
[2488] | 228 | |
---|
[5736] | 229 | # Make points relative to primary geo reference |
---|
[2594] | 230 | points[:,0] -= self.xllcorner |
---|
| 231 | points[:,1] -= self.yllcorner |
---|
| 232 | |
---|
[2488] | 233 | |
---|
| 234 | if is_list: |
---|
| 235 | points = points.tolist() |
---|
[2594] | 236 | |
---|
[2488] | 237 | return points |
---|
[3129] | 238 | |
---|
[2488] | 239 | |
---|
[3129] | 240 | def is_absolute(self): |
---|
| 241 | """Return True if xllcorner==yllcorner==0 indicating that points |
---|
| 242 | in question are absolute. |
---|
| 243 | """ |
---|
| 244 | |
---|
[6149] | 245 | return num.allclose([self.xllcorner, self.yllcorner], 0) |
---|
[3129] | 246 | |
---|
| 247 | |
---|
[2594] | 248 | |
---|
[2488] | 249 | def get_absolute(self, points): |
---|
| 250 | """ |
---|
| 251 | Given a set of points geo referenced to this instance, |
---|
| 252 | return the points as absolute values. |
---|
| 253 | """ |
---|
[3071] | 254 | |
---|
[3134] | 255 | #if self.is_absolute(): |
---|
| 256 | # return points |
---|
[2488] | 257 | is_list = False |
---|
| 258 | if type(points) == types.ListType: |
---|
| 259 | is_list = True |
---|
[3071] | 260 | |
---|
[6438] | 261 | points = ensure_numeric(copy.copy(points), num.Float) |
---|
[3071] | 262 | if len(points.shape) == 1: |
---|
| 263 | #One point has been passed |
---|
| 264 | msg = 'Single point must have two elements' |
---|
[3134] | 265 | if not len(points) == 2: |
---|
| 266 | raise ShapeError, msg |
---|
| 267 | #points = reshape(points, (1,2)) |
---|
[3071] | 268 | |
---|
| 269 | |
---|
| 270 | msg = 'Input must be an N x 2 array or list of (x,y) values. ' |
---|
| 271 | msg += 'I got an %d x %d array' %points.shape |
---|
[3134] | 272 | if not points.shape[1] == 2: |
---|
| 273 | raise ShapeError, msg |
---|
[2488] | 274 | |
---|
[3071] | 275 | |
---|
[3134] | 276 | # Add geo ref to points |
---|
| 277 | if not self.is_absolute(): |
---|
| 278 | points[:,0] += self.xllcorner |
---|
| 279 | points[:,1] += self.yllcorner |
---|
[3071] | 280 | |
---|
| 281 | |
---|
[2488] | 282 | if is_list: |
---|
| 283 | points = points.tolist() |
---|
| 284 | |
---|
| 285 | return points |
---|
| 286 | |
---|
[2594] | 287 | |
---|
[5288] | 288 | def get_relative(self, points): |
---|
| 289 | """ |
---|
| 290 | Given a set of points in absolute UTM coordinates, |
---|
| 291 | make them relative to this geo_reference instance, |
---|
| 292 | return the points as relative values. |
---|
| 293 | |
---|
| 294 | This is the inverse of get_absolute. |
---|
| 295 | """ |
---|
| 296 | |
---|
| 297 | is_list = False |
---|
| 298 | if type(points) == types.ListType: |
---|
| 299 | is_list = True |
---|
| 300 | |
---|
[6149] | 301 | points = ensure_numeric(points, num.Float) |
---|
[5288] | 302 | if len(points.shape) == 1: |
---|
| 303 | #One point has been passed |
---|
| 304 | msg = 'Single point must have two elements' |
---|
| 305 | if not len(points) == 2: |
---|
| 306 | raise ShapeError, msg |
---|
| 307 | #points = reshape(points, (1,2)) |
---|
| 308 | |
---|
| 309 | |
---|
| 310 | msg = 'Input must be an N x 2 array or list of (x,y) values. ' |
---|
| 311 | msg += 'I got an %d x %d array' %points.shape |
---|
| 312 | if not points.shape[1] == 2: |
---|
| 313 | raise ShapeError, msg |
---|
| 314 | |
---|
| 315 | |
---|
| 316 | # Subtract geo ref from points |
---|
| 317 | if not self.is_absolute(): |
---|
| 318 | points[:,0] -= self.xllcorner |
---|
| 319 | points[:,1] -= self.yllcorner |
---|
| 320 | |
---|
| 321 | |
---|
| 322 | if is_list: |
---|
| 323 | points = points.tolist() |
---|
| 324 | |
---|
| 325 | return points |
---|
| 326 | |
---|
| 327 | |
---|
| 328 | |
---|
[2594] | 329 | def reconcile_zones(self, other): |
---|
[4180] | 330 | if other is None: |
---|
| 331 | other = Geo_reference() |
---|
[3071] | 332 | if self.zone == other.zone or\ |
---|
| 333 | self.zone == DEFAULT_ZONE and other.zone == DEFAULT_ZONE: |
---|
[2594] | 334 | pass |
---|
| 335 | elif self.zone == DEFAULT_ZONE: |
---|
| 336 | self.zone = other.zone |
---|
| 337 | elif other.zone == DEFAULT_ZONE: |
---|
| 338 | other.zone = self.zone |
---|
| 339 | else: |
---|
[4470] | 340 | msg = 'Geospatial data must be in the same '+\ |
---|
[3614] | 341 | 'ZONE to allow reconciliation. I got zone %d and %d'\ |
---|
| 342 | %(self.zone, other.zone) |
---|
[2608] | 343 | raise ANUGAError, msg |
---|
[2488] | 344 | |
---|
| 345 | #def easting_northing2geo_reffed_point(self, x, y): |
---|
| 346 | # return [x-self.xllcorner, y - self.xllcorner] |
---|
| 347 | |
---|
| 348 | #def easting_northing2geo_reffed_points(self, x, y): |
---|
| 349 | # return [x-self.xllcorner, y - self.xllcorner] |
---|
| 350 | |
---|
| 351 | def get_origin(self): |
---|
| 352 | return (self.zone, self.xllcorner, self.yllcorner) |
---|
| 353 | |
---|
| 354 | def __repr__(self): |
---|
| 355 | return "(zone=%i easting=%f, northing=%f)" %(self.zone, self.xllcorner, self.yllcorner) |
---|
| 356 | |
---|
| 357 | def __cmp__(self,other): |
---|
| 358 | # FIXME (DSG) add a tolerence |
---|
| 359 | if other is None: |
---|
| 360 | return 1 |
---|
| 361 | cmp = 0 |
---|
| 362 | if not (self.xllcorner == self.xllcorner): |
---|
| 363 | cmp = 1 |
---|
| 364 | if not (self.yllcorner == self.yllcorner): |
---|
| 365 | cmp = 1 |
---|
| 366 | if not (self.zone == self.zone): |
---|
| 367 | cmp = 1 |
---|
| 368 | return cmp |
---|
[4387] | 369 | |
---|
| 370 | def write_NetCDF_georeference(origin, outfile): |
---|
| 371 | """ |
---|
| 372 | Write georeferrence info to a netcdf file, usually sww. |
---|
| 373 | |
---|
| 374 | The origin can be a georef instance or parrameters for a geo_ref instance |
---|
| 375 | |
---|
| 376 | outfile is the name of the file to be written to. |
---|
| 377 | """ |
---|
[4451] | 378 | geo_ref = ensure_geo_reference(origin) |
---|
| 379 | geo_ref.write_NetCDF(outfile) |
---|
| 380 | return geo_ref |
---|
| 381 | |
---|
| 382 | def ensure_geo_reference(origin): |
---|
| 383 | """ |
---|
| 384 | Given a list/tuple of zone, xllcorner and yllcorner of a geo-ref object, |
---|
| 385 | return a geo ref object. |
---|
| 386 | |
---|
| 387 | If the origin is None, return None, so calling this function doesn't |
---|
| 388 | effect code logic |
---|
| 389 | """ |
---|
[4387] | 390 | if isinstance(origin, Geo_reference): |
---|
| 391 | geo_ref = origin |
---|
[4451] | 392 | elif origin is None: |
---|
| 393 | geo_ref = None |
---|
[4387] | 394 | else: |
---|
[5730] | 395 | geo_ref = apply(Geo_reference, origin) |
---|
[4387] | 396 | return geo_ref |
---|
[4451] | 397 | |
---|
[4387] | 398 | |
---|
[2488] | 399 | #----------------------------------------------------------------------- |
---|
| 400 | |
---|
| 401 | if __name__ == "__main__": |
---|
| 402 | pass |
---|