[2488] | 1 | """ |
---|
| 2 | """ |
---|
| 3 | |
---|
| 4 | |
---|
| 5 | #FIXME: Ensure that all attributes of a georef are treated everywhere |
---|
| 6 | #and unit test |
---|
| 7 | |
---|
[2593] | 8 | import types, sys |
---|
[3129] | 9 | from Numeric import array, Float, ArrayType, reshape, allclose |
---|
[3514] | 10 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
| 11 | from anuga.utilities.anuga_exceptions import ANUGAError, TitleError, ParsingError, \ |
---|
[3134] | 12 | ShapeError |
---|
[2488] | 13 | |
---|
[2941] | 14 | |
---|
[2602] | 15 | DEFAULT_ZONE = -1 |
---|
[2488] | 16 | TITLE = '#geo reference' + "\n" #this title is referred to in the .xya format |
---|
| 17 | |
---|
| 18 | class Geo_reference: |
---|
| 19 | """ |
---|
| 20 | """ |
---|
| 21 | |
---|
| 22 | def __init__(self, |
---|
| 23 | zone = DEFAULT_ZONE, |
---|
| 24 | xllcorner = 0.0, |
---|
| 25 | yllcorner = 0.0, |
---|
| 26 | datum = 'wgs84', |
---|
| 27 | projection = 'UTM', units = 'm', |
---|
| 28 | false_easting = 500000, |
---|
| 29 | false_northing = 10000000, #Default for southern hemisphere |
---|
| 30 | NetCDFObject=None, |
---|
| 31 | ASCIIFile=None, |
---|
| 32 | read_title=None): |
---|
| 33 | """ |
---|
| 34 | input: |
---|
| 35 | NetCDFObject - a handle to the netCDF file to be written to |
---|
| 36 | ASCIIFile - a handle to the text file |
---|
| 37 | read_title - the title of the georeference text, if it was read in. |
---|
[2941] | 38 | If the function that calls this has already read the title line, |
---|
| 39 | it can't unread it, so this info has to be passed. |
---|
| 40 | If you know of a way to unread this info, then tell us. |
---|
[2488] | 41 | """ |
---|
| 42 | |
---|
| 43 | self.false_easting = false_easting |
---|
| 44 | self.false_northing = false_northing |
---|
| 45 | self.datum = datum |
---|
| 46 | self.projection = projection |
---|
| 47 | self.zone = zone |
---|
| 48 | self.units = units |
---|
| 49 | self.xllcorner = xllcorner |
---|
| 50 | self.yllcorner = yllcorner |
---|
| 51 | |
---|
| 52 | if NetCDFObject is not None: |
---|
| 53 | self.read_NetCDF(NetCDFObject) |
---|
| 54 | |
---|
| 55 | if ASCIIFile is not None: |
---|
| 56 | self.read_ASCII(ASCIIFile, read_title=read_title) |
---|
| 57 | |
---|
| 58 | def get_xllcorner(self): |
---|
| 59 | return self.xllcorner |
---|
| 60 | |
---|
| 61 | def get_yllcorner(self): |
---|
| 62 | return self.yllcorner |
---|
| 63 | |
---|
| 64 | def get_zone(self): |
---|
| 65 | return self.zone |
---|
| 66 | |
---|
| 67 | def write_NetCDF(self, outfile): |
---|
| 68 | outfile.xllcorner = self.xllcorner |
---|
| 69 | outfile.yllcorner = self.yllcorner |
---|
| 70 | outfile.zone = self.zone |
---|
| 71 | |
---|
| 72 | outfile.false_easting = self.false_easting |
---|
| 73 | outfile.false_northing = self.false_northing |
---|
| 74 | |
---|
| 75 | outfile.datum = self.datum |
---|
| 76 | outfile.projection = self.projection |
---|
| 77 | outfile.units = self.units |
---|
| 78 | |
---|
| 79 | def read_NetCDF(self, infile): |
---|
| 80 | self.xllcorner = infile.xllcorner[0] |
---|
| 81 | self.yllcorner = infile.yllcorner[0] |
---|
| 82 | self.zone = infile.zone[0] |
---|
| 83 | |
---|
| 84 | # Fix some assertion failures |
---|
| 85 | if type(self.zone) == ArrayType and self.zone.shape == (): |
---|
| 86 | self.zone = self.zone[0] |
---|
| 87 | if type(self.xllcorner) == ArrayType and self.xllcorner.shape == (): |
---|
| 88 | self.xllcorner = self.xllcorner[0] |
---|
| 89 | if type(self.yllcorner) == ArrayType and self.yllcorner.shape == (): |
---|
| 90 | self.yllcorner = self.yllcorner[0] |
---|
| 91 | |
---|
| 92 | assert (type(self.xllcorner) == types.FloatType or\ |
---|
| 93 | type(self.xllcorner) == types.IntType) |
---|
| 94 | assert (type(self.yllcorner) == types.FloatType or\ |
---|
| 95 | type(self.yllcorner) == types.IntType) |
---|
| 96 | assert (type(self.zone) == types.IntType) |
---|
| 97 | |
---|
| 98 | #FIXME (Ole) Read in the rest... |
---|
| 99 | |
---|
| 100 | |
---|
| 101 | def write_ASCII(self, fd): |
---|
| 102 | fd.write(TITLE) |
---|
| 103 | fd.write(str(self.zone) + "\n") |
---|
| 104 | fd.write(str(self.xllcorner) + "\n") |
---|
| 105 | fd.write(str(self.yllcorner) + "\n") |
---|
| 106 | |
---|
| 107 | def read_ASCII(self, fd,read_title=None): |
---|
[2942] | 108 | try: |
---|
| 109 | if read_title == None: |
---|
| 110 | read_title = fd.readline() # remove the title line |
---|
| 111 | if read_title[0:2].upper() != TITLE[0:2].upper(): |
---|
| 112 | msg = 'File error. Expecting line: %s. Got this line: %s' \ |
---|
| 113 | %(TITLE, read_title) |
---|
| 114 | raise TitleError, msg |
---|
| 115 | self.zone = int(fd.readline()) |
---|
| 116 | self.xllcorner = float(fd.readline()) |
---|
| 117 | self.yllcorner = float(fd.readline()) |
---|
| 118 | except SyntaxError: |
---|
| 119 | msg = 'File error. Got syntax error while parsing geo reference' |
---|
| 120 | raise ParsingError, msg |
---|
[2488] | 121 | |
---|
| 122 | # Fix some assertion failures |
---|
| 123 | if(type(self.zone) == ArrayType and self.zone.shape == ()): |
---|
| 124 | self.zone = self.zone[0] |
---|
| 125 | if type(self.xllcorner) == ArrayType and self.xllcorner.shape == (): |
---|
| 126 | self.xllcorner = self.xllcorner[0] |
---|
| 127 | if type(self.yllcorner) == ArrayType and self.yllcorner.shape == (): |
---|
| 128 | self.yllcorner = self.yllcorner[0] |
---|
| 129 | |
---|
| 130 | assert (type(self.xllcorner) == types.FloatType) |
---|
| 131 | assert (type(self.yllcorner) == types.FloatType) |
---|
| 132 | assert (type(self.zone) == types.IntType) |
---|
| 133 | |
---|
| 134 | #false_easting = infile.false_easting[0] |
---|
| 135 | #false_northing = infile.false_northing[0] |
---|
| 136 | |
---|
| 137 | def change_points_geo_ref(self, points, |
---|
| 138 | points_geo_ref=None): |
---|
| 139 | """ |
---|
| 140 | Change the geo reference of a list or Numeric array of points to |
---|
| 141 | be this reference.(The reference used for this object) |
---|
| 142 | If the points do not have a geo ref, assume 'absolute' values |
---|
| 143 | """ |
---|
[2594] | 144 | |
---|
[2488] | 145 | is_list = False |
---|
| 146 | if type(points) == types.ListType: |
---|
| 147 | is_list = True |
---|
[2594] | 148 | |
---|
[2595] | 149 | points = ensure_numeric(points, Float) |
---|
| 150 | |
---|
[2594] | 151 | if len(points.shape) == 1: |
---|
| 152 | #One point has been passed |
---|
| 153 | msg = 'Single point must have two elements' |
---|
| 154 | assert len(points) == 2, msg |
---|
| 155 | points = reshape(points, (1,2)) |
---|
| 156 | |
---|
[3071] | 157 | msg = 'Input must be an N x 2 array or list of (x,y) values. ' |
---|
| 158 | msg += 'I got an %d x %d array' %points.shape |
---|
| 159 | assert points.shape[1] == 2, msg |
---|
| 160 | |
---|
[2594] | 161 | |
---|
| 162 | if points_geo_ref is not self: |
---|
| 163 | #add point geo ref to points |
---|
| 164 | if not points_geo_ref is None: |
---|
| 165 | points[:,0] += points_geo_ref.xllcorner |
---|
| 166 | points[:,1] += points_geo_ref.yllcorner |
---|
[2488] | 167 | |
---|
[2594] | 168 | #subtract primary geo ref from points |
---|
| 169 | points[:,0] -= self.xllcorner |
---|
| 170 | points[:,1] -= self.yllcorner |
---|
| 171 | |
---|
[2488] | 172 | |
---|
| 173 | if is_list: |
---|
| 174 | points = points.tolist() |
---|
[2594] | 175 | |
---|
[2488] | 176 | return points |
---|
[3129] | 177 | |
---|
[2488] | 178 | |
---|
[3129] | 179 | def is_absolute(self): |
---|
| 180 | """Return True if xllcorner==yllcorner==0 indicating that points |
---|
| 181 | in question are absolute. |
---|
| 182 | """ |
---|
| 183 | |
---|
| 184 | return allclose([self.xllcorner, self.yllcorner], 0) |
---|
| 185 | |
---|
| 186 | |
---|
[2594] | 187 | |
---|
[2488] | 188 | def get_absolute(self, points): |
---|
| 189 | """ |
---|
| 190 | Given a set of points geo referenced to this instance, |
---|
| 191 | return the points as absolute values. |
---|
| 192 | """ |
---|
[3071] | 193 | |
---|
[3134] | 194 | #if self.is_absolute(): |
---|
| 195 | # return points |
---|
[3129] | 196 | |
---|
[3071] | 197 | |
---|
[2488] | 198 | is_list = False |
---|
| 199 | if type(points) == types.ListType: |
---|
| 200 | is_list = True |
---|
[3071] | 201 | |
---|
| 202 | points = ensure_numeric(points, Float) |
---|
| 203 | if len(points.shape) == 1: |
---|
| 204 | #One point has been passed |
---|
| 205 | msg = 'Single point must have two elements' |
---|
[3134] | 206 | if not len(points) == 2: |
---|
| 207 | raise ShapeError, msg |
---|
| 208 | #points = reshape(points, (1,2)) |
---|
[3071] | 209 | |
---|
| 210 | |
---|
| 211 | msg = 'Input must be an N x 2 array or list of (x,y) values. ' |
---|
| 212 | msg += 'I got an %d x %d array' %points.shape |
---|
[3134] | 213 | if not points.shape[1] == 2: |
---|
| 214 | raise ShapeError, msg |
---|
[2488] | 215 | |
---|
[3071] | 216 | |
---|
[3134] | 217 | # Add geo ref to points |
---|
| 218 | if not self.is_absolute(): |
---|
| 219 | points[:,0] += self.xllcorner |
---|
| 220 | points[:,1] += self.yllcorner |
---|
[3071] | 221 | |
---|
| 222 | |
---|
[2488] | 223 | if is_list: |
---|
| 224 | points = points.tolist() |
---|
| 225 | |
---|
| 226 | return points |
---|
| 227 | |
---|
[2594] | 228 | |
---|
| 229 | def reconcile_zones(self, other): |
---|
| 230 | |
---|
[3071] | 231 | if self.zone == other.zone or\ |
---|
| 232 | self.zone == DEFAULT_ZONE and other.zone == DEFAULT_ZONE: |
---|
[2594] | 233 | pass |
---|
| 234 | elif self.zone == DEFAULT_ZONE: |
---|
| 235 | self.zone = other.zone |
---|
| 236 | elif other.zone == DEFAULT_ZONE: |
---|
| 237 | other.zone = self.zone |
---|
| 238 | else: |
---|
| 239 | msg = 'Both geospatial_data objects must be in the same \ |
---|
[3071] | 240 | ZONE to allow reconciliation. I got zone %d and %d'\ |
---|
| 241 | %(self.zone, other.zone) |
---|
[2608] | 242 | raise ANUGAError, msg |
---|
[2488] | 243 | |
---|
| 244 | #def easting_northing2geo_reffed_point(self, x, y): |
---|
| 245 | # return [x-self.xllcorner, y - self.xllcorner] |
---|
| 246 | |
---|
| 247 | #def easting_northing2geo_reffed_points(self, x, y): |
---|
| 248 | # return [x-self.xllcorner, y - self.xllcorner] |
---|
| 249 | |
---|
| 250 | def get_origin(self): |
---|
| 251 | return (self.zone, self.xllcorner, self.yllcorner) |
---|
| 252 | |
---|
| 253 | def __repr__(self): |
---|
| 254 | return "(zone=%i easting=%f, northing=%f)" %(self.zone, self.xllcorner, self.yllcorner) |
---|
| 255 | |
---|
| 256 | def __cmp__(self,other): |
---|
| 257 | # FIXME (DSG) add a tolerence |
---|
| 258 | if other is None: |
---|
| 259 | return 1 |
---|
| 260 | cmp = 0 |
---|
| 261 | if not (self.xllcorner == self.xllcorner): |
---|
| 262 | cmp = 1 |
---|
| 263 | if not (self.yllcorner == self.yllcorner): |
---|
| 264 | cmp = 1 |
---|
| 265 | if not (self.zone == self.zone): |
---|
| 266 | cmp = 1 |
---|
| 267 | return cmp |
---|
| 268 | |
---|
| 269 | #----------------------------------------------------------------------- |
---|
| 270 | |
---|
| 271 | if __name__ == "__main__": |
---|
| 272 | pass |
---|