22'''
33Module to read DTM files published by Geoscience Australia
44Written by Stephen Dade ([email protected] 5- '''
65
7- import os
8- import sys
6+ AP_FLAKE8_CLEAN
7+ '''
98
109import numpy
10+ import os
1111
1212
1313class ERMap :
@@ -37,10 +37,10 @@ def read_ermapper(self, ifile):
3737 nroflines = int (self .header ['nroflines' ])
3838 nrofcellsperlines = int (self .header ['nrofcellsperline' ])
3939 self .data = self .read_ermapper_data (data_file , offset = int (self .header ['headeroffset' ]))
40- self .data = numpy .reshape (self .data ,(nroflines ,nrofcellsperlines ))
40+ self .data = numpy .reshape (self .data , (nroflines , nrofcellsperlines ))
4141
42- longy = numpy .fromstring (self .getHeaderParam ('longitude' ), sep = ':' )
43- latty = numpy .fromstring (self .getHeaderParam ('latitude' ), sep = ':' )
42+ longy = numpy .fromstring (self .getHeaderParam ('longitude' ), sep = ':' )
43+ latty = numpy .fromstring (self .getHeaderParam ('latitude' ), sep = ':' )
4444
4545 self .deltalatitude = float (self .header ['ydimension' ])
4646 self .deltalongitude = float (self .header ['xdimension' ])
@@ -62,32 +62,32 @@ def read_ermapper_header(self, ifile):
6262 # function for reading an ERMapper header from file
6363 header = {}
6464
65- fid = open (ifile ,'rt' )
65+ fid = open (ifile , 'rt' )
6666 header_string = fid .readlines ()
6767 fid .close ()
6868
6969 for line in header_string :
7070 if line .find ('=' ) > 0 :
7171 tmp_string = line .strip ().split ('=' )
72- header [tmp_string [0 ].strip ().lower ()]= tmp_string [1 ].strip ()
72+ header [tmp_string [0 ].strip ().lower ()] = tmp_string [1 ].strip ()
7373
7474 return header
7575
76- def read_ermapper_data (self , ifile , data_format = numpy .float32 , offset = 0 ):
76+ def read_ermapper_data (self , ifile , data_format = numpy .float32 , offset = 0 ):
7777 # open input file in a binary format and read the input string
78- fid = open (ifile ,'rb' )
78+ fid = open (ifile , 'rb' )
7979 if offset != 0 :
8080 fid .seek (offset )
8181 input_string = fid .read ()
8282 fid .close ()
8383
8484 # convert input string to required format (Note default format is numpy.float32)
85- grid_as_float = numpy .fromstring (input_string ,data_format )
85+ grid_as_float = numpy .fromstring (input_string , data_format )
8686 return grid_as_float
8787
8888 def getHeaderParam (self , key ):
89- '''Find a parameter in the associated .ers file'''
90- return self .header [key ]
89+ '''Find a parameter in the associated .ers file'''
90+ return self .header [key ]
9191
9292 def printBoundingBox (self ):
9393 '''Print the bounding box that this DEM covers'''
@@ -113,7 +113,7 @@ def getPercentBlank(self):
113113
114114 def getAltitudeAtPoint (self , latty , longy ):
115115 '''Return the altitude at a particular long/lat'''
116- #check the bounds
116+ # check the bounds
117117 if self .startlongitude > 0 and (longy < self .startlongitude or longy > self .endlongitude ):
118118 return - 1
119119 if self .startlongitude < 0 and (longy > self .startlongitude or longy < self .endlongitude ):
@@ -126,20 +126,20 @@ def getAltitudeAtPoint(self, latty, longy):
126126 x = numpy .abs ((latty - self .startlatitude )/ self .deltalatitude )
127127 y = numpy .abs ((longy - self .startlongitude )/ self .deltalongitude )
128128
129- #do some interpolation
129+ # do some interpolation
130130 # print("x,y", x, y)
131131 x_int = int (x )
132132 x_frac = x - int (x )
133133 y_int = int (y )
134134 y_frac = y - int (y )
135- #print("frac", x_int, x_frac, y_int, y_frac)
135+ # print("frac", x_int, x_frac, y_int, y_frac)
136136 value00 = self .data [x_int , y_int ]
137137 value10 = self .data [x_int + 1 , y_int ]
138138 value01 = self .data [x_int , y_int + 1 ]
139139 value11 = self .data [x_int + 1 , y_int + 1 ]
140- #print("values ", value00, value10, value01, value11)
140+ # print("values ", value00, value10, value01, value11)
141141
142- #check for null values
142+ # check for null values
143143 if value00 == - 99999 :
144144 value00 = 0
145145 if value10 == - 99999 :
@@ -151,7 +151,7 @@ def getAltitudeAtPoint(self, latty, longy):
151151
152152 value1 = self ._avg (value00 , value10 , x_frac )
153153 value2 = self ._avg (value01 , value11 , x_frac )
154- value = self ._avg (value1 , value2 , y_frac )
154+ value = self ._avg (value1 , value2 , y_frac )
155155
156156 return value
157157
@@ -167,20 +167,19 @@ def _avg(value1, value2, weight):
167167 return value2 * weight + value1 * (1 - weight )
168168
169169
170-
171170if __name__ == '__main__' :
172171
173172 print ("./Canberra/GSNSW_P756demg" )
174173 mappy = ERMap ()
175174 mappy .read_ermapper (os .path .join (os .environ ['HOME' ], './Documents/Elevation/Canberra/GSNSW_P756demg' ))
176175
177- #print some header data
176+ # print some header data
178177 mappy .printBoundingBox ()
179178
180- #get a measure of data quality
181- #mappy.getPercentBlank()
179+ # get a measure of data quality
180+ # mappy.getPercentBlank()
182181
183- #test the altitude (around Canberra):
182+ # test the altitude (around Canberra):
184183 alt = mappy .getAltitudeAtPoint (- 35.274411 , 149.097504 )
185184 print ("Alt at (-35.274411, 149.097504) is 807m (Google) or " + str (alt ))
186185 alt = mappy .getAltitudeAtPoint (- 35.239648 , 149.126118 )
@@ -197,7 +196,3 @@ def _avg(value1, value2, weight):
197196 print ("Alt at (-35.045126, 149.257482) is 577m (Google) or " + str (alt ))
198197 alt = mappy .getAltitudeAtPoint (- 35.564044 , 149.177657 )
199198 print ("Alt at (-35.564044, 149.177657) is 1113m (Google) or " + str (alt ))
200-
201-
202-
203-
0 commit comments