| 1 | """ |
|---|
| 2 | """ |
|---|
| 3 | |
|---|
| 4 | |
|---|
| 5 | #FIXME: Ensure that all attributes of a georef are treated everywhere |
|---|
| 6 | #and unit test |
|---|
| 7 | |
|---|
| 8 | import types |
|---|
| 9 | from Numeric import array,Float,ArrayType |
|---|
| 10 | |
|---|
| 11 | DEFAULT_ZONE = 56 |
|---|
| 12 | TITLE = '#geo reference' + "\n" #this title is referred to in the .xya format |
|---|
| 13 | |
|---|
| 14 | class Geo_reference: |
|---|
| 15 | """ |
|---|
| 16 | """ |
|---|
| 17 | |
|---|
| 18 | def __init__(self, |
|---|
| 19 | zone = DEFAULT_ZONE, |
|---|
| 20 | xllcorner = 0.0, |
|---|
| 21 | yllcorner = 0.0, |
|---|
| 22 | datum = 'wgs84', |
|---|
| 23 | projection = 'UTM', |
|---|
| 24 | units = 'm', |
|---|
| 25 | false_easting = 500000, |
|---|
| 26 | false_northing = 10000000, #Default for southern hemisphere |
|---|
| 27 | NetCDFObject=None, |
|---|
| 28 | ASCIIFile=None, |
|---|
| 29 | read_title=None): |
|---|
| 30 | """ |
|---|
| 31 | input: |
|---|
| 32 | NetCDFObject - a handle to the netCDF file |
|---|
| 33 | ASCIIFile - a handle to the text file |
|---|
| 34 | read_title - the title of the georeference text, if it was read in. |
|---|
| 35 | """ |
|---|
| 36 | |
|---|
| 37 | self.false_easting = false_easting |
|---|
| 38 | self.false_northing = false_northing |
|---|
| 39 | self.datum = datum |
|---|
| 40 | self.projection = projection |
|---|
| 41 | self.zone = zone |
|---|
| 42 | self.units = units |
|---|
| 43 | self.xllcorner = xllcorner |
|---|
| 44 | self.yllcorner = yllcorner |
|---|
| 45 | |
|---|
| 46 | if NetCDFObject is not None: |
|---|
| 47 | self.read_NetCDF(NetCDFObject) |
|---|
| 48 | |
|---|
| 49 | if ASCIIFile is not None: |
|---|
| 50 | self.read_ASCII(ASCIIFile, read_title=read_title) |
|---|
| 51 | |
|---|
| 52 | def get_xllcorner(self): |
|---|
| 53 | return self.xllcorner |
|---|
| 54 | |
|---|
| 55 | def get_yllcorner(self): |
|---|
| 56 | return self.yllcorner |
|---|
| 57 | |
|---|
| 58 | def get_zone(self): |
|---|
| 59 | return self.zone |
|---|
| 60 | |
|---|
| 61 | def write_NetCDF(self, outfile): |
|---|
| 62 | outfile.xllcorner = self.xllcorner |
|---|
| 63 | outfile.yllcorner = self.yllcorner |
|---|
| 64 | outfile.zone = self.zone |
|---|
| 65 | |
|---|
| 66 | outfile.false_easting = self.false_easting |
|---|
| 67 | outfile.false_northing = self.false_northing |
|---|
| 68 | |
|---|
| 69 | outfile.datum = self.datum |
|---|
| 70 | outfile.projection = self.projection |
|---|
| 71 | outfile.units = self.units |
|---|
| 72 | |
|---|
| 73 | def read_NetCDF(self, infile): |
|---|
| 74 | self.xllcorner = infile.xllcorner[0] |
|---|
| 75 | self.yllcorner = infile.yllcorner[0] |
|---|
| 76 | self.zone = infile.zone[0] |
|---|
| 77 | |
|---|
| 78 | |
|---|
| 79 | assert (type(self.xllcorner) == types.FloatType or\ |
|---|
| 80 | type(self.xllcorner) == types.IntType) |
|---|
| 81 | assert (type(self.yllcorner) == types.FloatType or\ |
|---|
| 82 | type(self.yllcorner) == types.IntType) |
|---|
| 83 | assert (type(self.zone) == types.IntType) |
|---|
| 84 | |
|---|
| 85 | #FIXME (Ole) Read in the rest... |
|---|
| 86 | |
|---|
| 87 | |
|---|
| 88 | def write_ASCII(self, fd): |
|---|
| 89 | fd.write(TITLE) |
|---|
| 90 | fd.write(str(self.zone) + "\n") |
|---|
| 91 | fd.write(str(self.xllcorner) + "\n") |
|---|
| 92 | fd.write(str(self.yllcorner) + "\n") |
|---|
| 93 | |
|---|
| 94 | def read_ASCII(self, fd,read_title=None): |
|---|
| 95 | if read_title == None: |
|---|
| 96 | read_title = fd.readline() # remove the title line |
|---|
| 97 | #print "gr read_title[0:2].upper()",read_title[0:2].upper() |
|---|
| 98 | #print "gr TITLE[0:2].upper()",TITLE[0:2].upper() |
|---|
| 99 | if read_title[0:2].upper() != TITLE[0:2].upper(): |
|---|
| 100 | raise SyntaxError |
|---|
| 101 | self.zone = int(fd.readline()) |
|---|
| 102 | self.xllcorner = float(fd.readline()) |
|---|
| 103 | self.yllcorner = float(fd.readline()) |
|---|
| 104 | |
|---|
| 105 | |
|---|
| 106 | assert (type(self.xllcorner) == types.FloatType) |
|---|
| 107 | assert (type(self.yllcorner) == types.FloatType) |
|---|
| 108 | assert (type(self.zone) == types.IntType) |
|---|
| 109 | |
|---|
| 110 | #false_easting = infile.false_easting[0] |
|---|
| 111 | #false_northing = infile.false_northing[0] |
|---|
| 112 | |
|---|
| 113 | def change_points_geo_ref(self, points, |
|---|
| 114 | points_geo_ref=None): |
|---|
| 115 | """ |
|---|
| 116 | Change the reference of a list or Numeric array of points. |
|---|
| 117 | If the points do not have ageo ref, assume 'absolute' values |
|---|
| 118 | """ |
|---|
| 119 | is_list = False |
|---|
| 120 | if type(points) == types.ListType: |
|---|
| 121 | is_list = True |
|---|
| 122 | if len(points)>0 and type(points[0]) not in [types.ListType,types.TupleType]: |
|---|
| 123 | #a single pont is being passed. make it a list of lists |
|---|
| 124 | points = [points] |
|---|
| 125 | elif type(points) == ArrayType: |
|---|
| 126 | if len(points.shape) == 1: |
|---|
| 127 | points = [points] |
|---|
| 128 | |
|---|
| 129 | |
|---|
| 130 | #convert into array |
|---|
| 131 | points = array(points).astype(Float) |
|---|
| 132 | #add point geo ref to points |
|---|
| 133 | if not points_geo_ref is None: |
|---|
| 134 | points[:,0] += points_geo_ref.xllcorner |
|---|
| 135 | points[:,1] += points_geo_ref.yllcorner |
|---|
| 136 | |
|---|
| 137 | #subtract primary geo ref from points |
|---|
| 138 | points[:,0] -= self.xllcorner |
|---|
| 139 | points[:,1] -= self.yllcorner |
|---|
| 140 | if is_list: |
|---|
| 141 | points = points.tolist() |
|---|
| 142 | # return points |
|---|
| 143 | return points |
|---|
| 144 | |
|---|
| 145 | |
|---|
| 146 | #def easting_northing2geo_reffed_point(self, x, y): |
|---|
| 147 | # return [x-self.xllcorner, y - self.xllcorner] |
|---|
| 148 | |
|---|
| 149 | #def easting_northing2geo_reffed_points(self, x, y): |
|---|
| 150 | # return [x-self.xllcorner, y - self.xllcorner] |
|---|
| 151 | |
|---|
| 152 | def get_origin(self): |
|---|
| 153 | return (self.zone, self.xllcorner, self.yllcorner) |
|---|
| 154 | |
|---|
| 155 | def __repr__(self): |
|---|
| 156 | return "(%i %f,%f)" % (self.zone, self.xllcorner, self.yllcorner) |
|---|
| 157 | |
|---|
| 158 | def __cmp__(self,other): |
|---|
| 159 | # FIXME (DSG) add a tolerence |
|---|
| 160 | if other is None: |
|---|
| 161 | return 1 |
|---|
| 162 | cmp = 0 |
|---|
| 163 | if not (self.xllcorner == self.xllcorner): |
|---|
| 164 | cmp = 1 |
|---|
| 165 | if not (self.yllcorner == self.yllcorner): |
|---|
| 166 | cmp = 1 |
|---|
| 167 | if not (self.zone == self.zone): |
|---|
| 168 | cmp = 1 |
|---|
| 169 | return cmp |
|---|
| 170 | |
|---|
| 171 | #----------------------------------------------------------------------- |
|---|
| 172 | |
|---|
| 173 | if __name__ == "__main__": |
|---|
| 174 | pass |
|---|
