<s> import os <EOL> import re <EOL> import struct <EOL> from math import floor <EOL> from graphserver . vincenty import vincenty <EOL> def floatrange ( start , stop , step ) : <EOL> i = start <EOL> while i <= stop : <EOL> yield i <EOL> i += step <EOL> def cons ( ary ) : <EOL> for i in range ( len ( ary ) - <NUM_LIT:1> ) : <EOL> yield ary [ i ] , ary [ i + <NUM_LIT:1> ] <EOL> def split_line_segment ( lng1 , lat1 , lng2 , lat2 , max_section_length ) : <EOL> if lng1 == lng2 and lat1 == lat2 : <EOL> yield [ lng1 , lat1 , <NUM_LIT:0> ] <EOL> yield [ lng2 , lat2 , <NUM_LIT:0> ] <EOL> return <EOL> street_len = vincenty ( lat1 , lng1 , lat2 , lng2 ) <EOL> n_sections = int ( street_len / max_section_length ) + <NUM_LIT:1> <EOL> geolen = ( ( lat2 - lat1 ) ** <NUM_LIT:2> + ( lng2 - lng1 ) ** <NUM_LIT:2> ) ** <NUM_LIT:0.5> <EOL> section_len = geolen / n_sections <EOL> street_vector = ( lng2 - lng1 , lat2 - lat1 ) <EOL> unit_vector = [ x / geolen for x in street_vector ] <EOL> for i in range ( n_sections + <NUM_LIT:1> ) : <EOL> vec = [ x * section_len * i for x in unit_vector ] <EOL> vec = [ lng1 + vec [ <NUM_LIT:0> ] , lat1 + vec [ <NUM_LIT:1> ] , ( street_len / n_sections ) * i ] <EOL> yield vec <EOL> def split_line_string ( points , max_section_length ) : <EOL> split_segs = [ ] <EOL> for ( lng1 , lat1 ) , ( lng2 , lat2 ) in cons ( points ) : <EOL> split_seg = list ( split_line_segment ( lng1 , lat1 , lng2 , lat2 , max_section_length ) ) <EOL> split_segs . append ( split_seg ) <EOL> ret = [ ] <EOL> segstart_s = <NUM_LIT:0> <EOL> for i , split_seg in enumerate ( split_segs ) : <EOL> for x , y , s in split_seg [ : - <NUM_LIT:1> ] : <EOL> ret . append ( ( x , y , s + segstart_s ) ) <EOL> if i == len ( split_segs ) - <NUM_LIT:1> : <EOL> x , y , s = split_seg [ - <NUM_LIT:1> ] <EOL> ret . append ( ( x , y , s + segstart_s ) ) <EOL> segstart_s += split_seg [ - <NUM_LIT:1> ] [ <NUM_LIT:2> ] <EOL> return ret <EOL> class GridFloat : <EOL> def __init__ ( self , basename ) : <EOL> self . _read_header ( basename + "<STR_LIT>" ) <EOL> self . fp = open ( basename + "<STR_LIT>" , "<STR_LIT:rb>" ) <EOL> def _read_header ( self , filename ) : <EOL> fp = open ( filename , "<STR_LIT:r>" ) <EOL> self . ncols = int ( fp . readline ( ) [ <NUM_LIT> : ] . strip ( ) ) <EOL> self . nrows = int ( fp . readline ( ) [ <NUM_LIT> : ] . strip ( ) ) <EOL> self . xllcorner = float ( fp . readline ( ) [ <NUM_LIT> : ] . strip ( ) ) <EOL> self . yllcorner = float ( fp . readline ( ) [ <NUM_LIT> : ] . strip ( ) ) <EOL> self . cellsize = float ( fp . readline ( ) [ <NUM_LIT> : ] . strip ( ) ) <EOL> self . NODATA_value = int ( fp . readline ( ) [ <NUM_LIT> : ] . strip ( ) ) <EOL> self . byteorder = "<STR_LIT:<>" if fp . readline ( ) [ <NUM_LIT> : ] . strip ( ) == "<STR_LIT>" else "<STR_LIT:>>" <EOL> self . left = self . xllcorner <EOL> self . right = self . xllcorner + ( self . ncols - <NUM_LIT:1> ) * self . cellsize <EOL> self . bottom = self . yllcorner <EOL> self . top = self . yllcorner + ( self . nrows - <NUM_LIT:1> ) * self . cellsize <EOL> @ property <EOL> def extent ( self ) : <EOL> return ( self . xllcorner , <EOL> self . yllcorner , <EOL> self . xllcorner + self . cellsize * ( self . ncols - <NUM_LIT:1> ) , <EOL> self . yllcorner + self . cellsize * ( self . nrows - <NUM_LIT:1> ) ) <EOL> def contains ( self , lng , lat ) : <EOL> return not ( lng < self . left or lng >= self . right or lat <= self . bottom or lat > self . top ) <EOL> def allcells ( self ) : <EOL> self . fp . seek ( <NUM_LIT:0> ) <EOL> return struct . unpack ( "<STR_LIT>" % ( self . byteorder , self . nrows * self . ncols ) , self . fp . read ( ) ) <EOL> def extremes ( self ) : <EOL> mem = self . allcells ( ) <EOL> return ( min ( mem ) , max ( mem ) ) <EOL> def cell ( self , x , y ) : <EOL> position = ( y * self . ncols + x ) * <NUM_LIT:4> <EOL> self . fp . seek ( position ) <EOL> return struct . unpack ( "<STR_LIT>" % ( self . byteorder ) , self . fp . read ( <NUM_LIT:4> ) ) [ <NUM_LIT:0> ] <EOL> def elevation ( self , lng , lat , interpolate = True ) : <EOL> if lng < self . left or lng >= self . right or lat <= self . bottom or lat > self . top : <EOL> return None <EOL> x = ( lng - self . left ) / self . cellsize <EOL> y = ( self . top - lat ) / self . cellsize <EOL> ulx = int ( floor ( x ) ) <EOL> uly = int ( floor ( y ) ) <EOL> ul = self . cell ( ulx , uly ) <EOL> if not interpolate : <EOL> return ul <EOL> ur = self . cell ( ulx + <NUM_LIT:1> , uly ) <EOL> ll = self . cell ( ulx , uly + <NUM_LIT:1> ) <EOL> lr = self . cell ( ulx + <NUM_LIT:1> , uly + <NUM_LIT:1> ) <EOL> cellleft = x % <NUM_LIT:1> <EOL> celltop = y % <NUM_LIT:1> <EOL> um = ( ur - ul ) * cellleft + ul <EOL> lm = ( lr - ll ) * cellleft + ll <EOL> return ( lm - um ) * celltop + um <EOL> def profile ( self , points , resolution = <NUM_LIT:10> ) : <EOL> return [ ( s , self . elevation ( lng , lat ) ) for lng , lat , s in split_line_string ( points , resolution ) ] <EOL> class BIL : <EOL> def __init__ ( self , basename ) : <EOL> self . _read_header ( basename + "<STR_LIT>" ) <EOL> self . fp = open ( basename + "<STR_LIT>" , "<STR_LIT:rb>" ) <EOL> def _read_header ( self , filename ) : <EOL> HCW = <NUM_LIT:15> <EOL> fp = open ( filename , "<STR_LIT:r>" ) <EOL> raw_header = dict ( [ x . strip ( ) . split ( ) for x in fp . read ( ) . strip ( ) . split ( "<STR_LIT:\n>" ) ] ) <EOL> self . byteorder = "<STR_LIT:<>" if raw_header [ '<STR_LIT>' ] == "<STR_LIT:I>" else "<STR_LIT:>>" <EOL> self . layout = raw_header [ '<STR_LIT>' ] <EOL> self . ncols = int ( raw_header [ '<STR_LIT>' ] ) <EOL> self . nrows = int ( raw_header [ '<STR_LIT>' ] ) <EOL> self . nbands = int ( raw_header [ '<STR_LIT>' ] ) <EOL> self . nbits = int ( raw_header [ '<STR_LIT>' ] ) <EOL> self . bandrowbytes = int ( raw_header [ '<STR_LIT>' ] ) <EOL> self . totalrowbytes = int ( raw_header [ '<STR_LIT>' ] ) <EOL> self . pixeltype = raw_header [ '<STR_LIT>' ] <EOL> self . ulxmap = float ( raw_header [ '<STR_LIT>' ] ) <EOL> self . ulymap = float ( raw_header [ '<STR_LIT>' ] ) <EOL> self . xdim = float ( raw_header [ '<STR_LIT>' ] ) <EOL> self . ydim = float ( raw_header [ '<STR_LIT>' ] ) <EOL> self . nodata = float ( raw_header [ '<STR_LIT>' ] ) <EOL> self . left = self . ulxmap <EOL> self . right = self . ulxmap + ( self . ncols - <NUM_LIT:1> ) * self . xdim <EOL> self . bottom = self . ulymap - ( self . nrows - <NUM_LIT:1> ) * self . ydim <EOL> self . top = self . ulymap <EOL> @ property <EOL> def extent ( self ) : <EOL> return ( self . left , <EOL> self . bottom , <EOL> self . right , <EOL> self . top ) <EOL> def contains ( self , lng , lat ) : <EOL> return not ( lng < self . left or lng >= self . right or lat <= self . bottom or lat > self . top ) <EOL> def allcells ( self ) : <EOL> self . fp . seek ( <NUM_LIT:0> ) <EOL> return struct . unpack ( "<STR_LIT>" % ( self . byteorder , self . nrows * self . ncols ) , self . fp . read ( ) ) <EOL> def extremes ( self ) : <EOL> mem = self . allcells ( ) <EOL> return ( min ( mem ) , max ( mem ) ) <EOL> def cell ( self , x , y ) : <EOL> position = ( y * self . ncols + x ) * <NUM_LIT:4> <EOL> self . fp . seek ( position ) <EOL> return struct . unpack ( "<STR_LIT>" % ( self . byteorder ) , self . fp . read ( <NUM_LIT:4> ) ) [ <NUM_LIT:0> ] <EOL> def elevation ( self , lng , lat , interpolate = True ) : <EOL> if lng < self . left or lng >= self . right or lat <= self . bottom or lat > self . top : <EOL> return None <EOL> x = ( lng - self . left ) / self . xdim <EOL> y = ( self . top - lat ) / self . ydim <EOL> ulx = int ( floor ( x ) ) <EOL> uly = int ( floor ( y ) ) <EOL> ul = self . cell ( ulx , uly ) <EOL> if not interpolate : <EOL> return ul <EOL> ur = self . cell ( ulx + <NUM_LIT:1> , uly ) <EOL> ll = self . cell ( ulx , uly + <NUM_LIT:1> ) <EOL> lr = self . cell ( ulx + <NUM_LIT:1> , uly + <NUM_LIT:1> ) <EOL> cellleft = x % <NUM_LIT:1> <EOL> celltop = y % <NUM_LIT:1> <EOL> um = ( ur - ul ) * cellleft + ul <EOL> lm = ( lr - ll ) * cellleft + ll <EOL> return ( lm - um ) * celltop + um <EOL> def profile ( self , points , resolution = <NUM_LIT:10> ) : <EOL> return [ ( s , self . elevation ( lng , lat ) ) for lng , lat , s in split_line_string ( points , resolution ) ] <EOL> class ElevationPile : <EOL> def __init__ ( self ) : <EOL> self . tiles = [ ] <EOL> def add ( self , dem_basename ) : <EOL> base_basename = "<STR_LIT>" . join ( dem_basename . split ( "<STR_LIT:.>" ) [ <NUM_LIT:0> : - <NUM_LIT:1> ] ) <EOL> format = dem_basename . split ( "<STR_LIT:.>" ) [ - <NUM_LIT:1> ] <EOL> if format == "<STR_LIT>" : <EOL> dem = GridFloat ( base_basename ) <EOL> elif format == "<STR_LIT>" : <EOL> dem = BIL ( base_basename ) <EOL> else : <EOL> raise Exception ( "<STR_LIT>" % format ) <EOL> self . tiles . append ( dem ) <EOL> def elevation ( self , lng , lat , interpolate = True ) : <EOL> for tile in self . tiles : <EOL> if tile . contains ( lng , lat ) : <EOL> return tile . elevation ( lng , lat , interpolate ) <EOL> def profile ( self , points , resolution = <NUM_LIT:10> ) : <EOL> return [ ( s , self . elevation ( lng , lat ) ) for lng , lat , s in split_line_string ( points , resolution ) ] <EOL> def selftest ( ) : <EOL> BASENAME = "<STR_LIT>" <EOL> HOMEAREA = "<STR_LIT>" + BASENAME <EOL> gf = GridFloat ( HOMEAREA , BASENAME ) <EOL> print gf <EOL> print gf . extent <EOL> toprow = [ gf . cell ( x , <NUM_LIT:0> ) for x in range ( gf . ncols ) ] <EOL> assert gf . elevation ( gf . left , gf . top ) == toprow [ <NUM_LIT:0> ] <EOL> assert round ( gf . elevation ( gf . right - <NUM_LIT> , gf . top ) , <NUM_LIT:5> ) == round ( toprow [ - <NUM_LIT:1> ] , <NUM_LIT:5> ) <EOL> bottomrow = [ gf . cell ( x , gf . nrows - <NUM_LIT:2> ) for x in range ( gf . ncols ) ] <EOL> assert gf . elevation ( gf . left , gf . bottom + <NUM_LIT> ) == bottomrow [ <NUM_LIT:0> ] <EOL> assert gf . elevation ( gf . right - <NUM_LIT> , gf . bottom + <NUM_LIT> ) == bottomrow [ - <NUM_LIT:2> ] <EOL> assert gf . extremes ( ) == ( <NUM_LIT> , <NUM_LIT> ) <EOL> assert round ( gf . elevation ( ( gf . right - gf . left ) / <NUM_LIT:2> + gf . left , ( gf . top - gf . bottom ) / <NUM_LIT:2> + gf . bottom ) , <NUM_LIT:6> ) == round ( <NUM_LIT> , <NUM_LIT:6> ) <EOL> def create_elev_circles ( ) : <EOL> from renderer . processing import MapRenderer <EOL> BASENAME = "<STR_LIT>" <EOL> HOMEAREA = "<STR_LIT>" + BASENAME <EOL> gf = GridFloat ( HOMEAREA , BASENAME ) <EOL> mr = MapRenderer ( "<STR_LIT>" ) <EOL> mr . start ( gf . left , gf . bottom , gf . right , gf . top , <NUM_LIT> ) <EOL> mr . smooth ( ) <EOL> mr . fill ( <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ) <EOL> mr . background ( <NUM_LIT:255> , <NUM_LIT:255> , <NUM_LIT:255> ) <EOL> mr . strokeWeight ( <NUM_LIT> ) <EOL> for y in floatrange ( gf . bottom , gf . top , ( gf . top - gf . bottom ) / <NUM_LIT:50> ) : <EOL> for x in floatrange ( gf . left , gf . right , ( gf . right - gf . left ) / <NUM_LIT:50> ) : <EOL> elev = gf . elevation ( x , y ) <EOL> mr . ellipse ( x , y , elev * <NUM_LIT> , elev * <NUM_LIT> ) <EOL> mr . saveLocal ( "<STR_LIT>" ) <EOL> mr . stop ( ) <EOL> if __name__ == '<STR_LIT:__main__>' : <EOL> BASENAME = "<STR_LIT>" <EOL> HOMEAREA = "<STR_LIT>" + BASENAME <EOL> gf = GridFloat ( HOMEAREA , BASENAME ) <EOL> for x in gf . extent : <EOL> print x </s>
