""" """ #FIXME: Ensure that all attributes of a georef are treated everywhere #and unit test import types, sys from Numeric import array, Float, ArrayType, reshape from utilities.numerical_tools import ensure_numeric from utilities.anuga_exceptions import ANUGAError DEFAULT_ZONE = -1 TITLE = '#geo reference' + "\n" #this title is referred to in the .xya format class Geo_reference: """ """ def __init__(self, zone = DEFAULT_ZONE, xllcorner = 0.0, yllcorner = 0.0, datum = 'wgs84', projection = 'UTM', units = 'm', false_easting = 500000, false_northing = 10000000, #Default for southern hemisphere NetCDFObject=None, ASCIIFile=None, read_title=None): """ input: NetCDFObject - a handle to the netCDF file to be written to ASCIIFile - a handle to the text file read_title - the title of the georeference text, if it was read in. """ self.false_easting = false_easting self.false_northing = false_northing self.datum = datum self.projection = projection self.zone = zone self.units = units self.xllcorner = xllcorner self.yllcorner = yllcorner if NetCDFObject is not None: self.read_NetCDF(NetCDFObject) if ASCIIFile is not None: self.read_ASCII(ASCIIFile, read_title=read_title) def get_xllcorner(self): return self.xllcorner def get_yllcorner(self): return self.yllcorner def get_zone(self): return self.zone def write_NetCDF(self, outfile): outfile.xllcorner = self.xllcorner outfile.yllcorner = self.yllcorner outfile.zone = self.zone outfile.false_easting = self.false_easting outfile.false_northing = self.false_northing outfile.datum = self.datum outfile.projection = self.projection outfile.units = self.units def read_NetCDF(self, infile): self.xllcorner = infile.xllcorner[0] self.yllcorner = infile.yllcorner[0] self.zone = infile.zone[0] # Fix some assertion failures if type(self.zone) == ArrayType and self.zone.shape == (): self.zone = self.zone[0] if type(self.xllcorner) == ArrayType and self.xllcorner.shape == (): self.xllcorner = self.xllcorner[0] if type(self.yllcorner) == ArrayType and self.yllcorner.shape == (): self.yllcorner = self.yllcorner[0] assert (type(self.xllcorner) == types.FloatType or\ type(self.xllcorner) == types.IntType) assert (type(self.yllcorner) == types.FloatType or\ type(self.yllcorner) == types.IntType) assert (type(self.zone) == types.IntType) #FIXME (Ole) Read in the rest... def write_ASCII(self, fd): fd.write(TITLE) fd.write(str(self.zone) + "\n") fd.write(str(self.xllcorner) + "\n") fd.write(str(self.yllcorner) + "\n") def read_ASCII(self, fd,read_title=None): if read_title == None: read_title = fd.readline() # remove the title line #print "gr read_title[0:2].upper()",read_title[0:2].upper() #print "gr TITLE[0:2].upper()",TITLE[0:2].upper() if read_title[0:2].upper() != TITLE[0:2].upper(): raise SyntaxError self.zone = int(fd.readline()) self.xllcorner = float(fd.readline()) self.yllcorner = float(fd.readline()) # Fix some assertion failures if(type(self.zone) == ArrayType and self.zone.shape == ()): self.zone = self.zone[0] if type(self.xllcorner) == ArrayType and self.xllcorner.shape == (): self.xllcorner = self.xllcorner[0] if type(self.yllcorner) == ArrayType and self.yllcorner.shape == (): self.yllcorner = self.yllcorner[0] assert (type(self.xllcorner) == types.FloatType) assert (type(self.yllcorner) == types.FloatType) assert (type(self.zone) == types.IntType) #false_easting = infile.false_easting[0] #false_northing = infile.false_northing[0] def change_points_geo_ref(self, points, points_geo_ref=None): """ Change the geo reference of a list or Numeric array of points to be this reference.(The reference used for this object) If the points do not have a geo ref, assume 'absolute' values """ # FIXME: Use ensure_numeric is_list = False if type(points) == types.ListType: is_list = True points = ensure_numeric(points, Float) if len(points.shape) == 1: #One point has been passed msg = 'Single point must have two elements' assert len(points) == 2, msg points = reshape(points, (1,2)) if points_geo_ref is not self: #add point geo ref to points if not points_geo_ref is None: points[:,0] += points_geo_ref.xllcorner points[:,1] += points_geo_ref.yllcorner #subtract primary geo ref from points points[:,0] -= self.xllcorner points[:,1] -= self.yllcorner if is_list: points = points.tolist() return points def get_absolute(self, points): """ Given a set of points geo referenced to this instance, return the points as absolute values. """ is_list = False if type(points) == types.ListType: is_list = True if len(points)>0 and type(points[0]) \ not in [types.ListType,types.TupleType]: #a single point is being passed. make it a list of lists points = [points] elif type(points) == ArrayType: if len(points.shape) == 1: points = [points] # convert into array points = array(points).astype(Float) # add primary geo ref from points points[:,0] += self.xllcorner points[:,1] += self.yllcorner if is_list: points = points.tolist() return points def reconcile_zones(self, other): if self.zone == other.zone or self.zone == DEFAULT_ZONE and other.zone == DEFAULT_ZONE: pass elif self.zone == DEFAULT_ZONE: self.zone = other.zone elif other.zone == DEFAULT_ZONE: other.zone = self.zone else: msg = 'Both geospatial_data objects must be in the same \ ZONE to allow reconciliation. I got zone %d and %d' %(self.zone, other.zone) raise ANUGAError, msg #def easting_northing2geo_reffed_point(self, x, y): # return [x-self.xllcorner, y - self.xllcorner] #def easting_northing2geo_reffed_points(self, x, y): # return [x-self.xllcorner, y - self.xllcorner] def get_origin(self): return (self.zone, self.xllcorner, self.yllcorner) def __repr__(self): return "(zone=%i easting=%f, northing=%f)" %(self.zone, self.xllcorner, self.yllcorner) def __cmp__(self,other): # FIXME (DSG) add a tolerence if other is None: return 1 cmp = 0 if not (self.xllcorner == self.xllcorner): cmp = 1 if not (self.yllcorner == self.yllcorner): cmp = 1 if not (self.zone == self.zone): cmp = 1 return cmp #----------------------------------------------------------------------- if __name__ == "__main__": pass