From GeoMod
from Numeric import *
from string import *
from visual import *
'''A raster map crreated with a 2d array'''
class raster_map:
def __init__(self, data_array, xmin=0.0, ymax=0.0, dx=1): #datarr, xmin, ymax, dx, nrows, ncols):
'''data_array is a 2d array '''
self.data = data_array
(self.nrows, self.ncols) = self.data.shape
self.xmin = xmin
self.ymax = ymax
self.dx = dx
def get_x(self, c):
return self.xmin + c * self.dx
def get_y(self, r):
return self.ymax - r * self.dx
def get_val(self, pos):
r, c = self.get_rc(pos)
return self.data[r,c]
def get_rc(self, pos):
r = int(((self.ymax - pos.y + self.dx/2.0) ) /self.dx )
c = int(((pos.x - self.xmin + self.dx/2.0) ) /self.dx)
return (r, c)
def draw_map(self, imin=0, imax=-99999, jmin=0, jmax=-99999, scale = 1.0, center = 0):
''' draw a specific region of the raster map'''
self.scale = scale
default_val = -99999
if imax == default_val:
imax = self.nrows
if jmax == default_val:
jmax = self.ncols
print "drawing map "
print "imin, imax ", imin, imax
print "jmin, jmax ", jmin, jmax
tdata = self.data[imin:imax,jmin:jmax]
#draw in topography
tmax = max(max(tdata))
tmin = min(min(tdata))
for i in range(imin,imax):
for j in range(jmin,jmax):
ipos = self.xmin + j* self.dx
jpos = self.ymax - i* self.dx
kpos = self.data[i,j] #self.get_val(vector(ipos,jpos,0))
colval = (kpos - tmin) / (tmax - tmin)
box(pos=(ipos,jpos, 0.5*kpos*scale),
length = self.dx, height = self.dx, width = self.data[i,j]*scale,
color = (colval,colval,colval))
#center scene on area
#print self.xmin + jmin* self.dx, self.ymax - imax* self.dx
if center == 1:
xc = self.xmin + (self.dx * ((jmin+jmax)/2))
yc = self.ymax - (self.dx * ((imin+imax)/2))
zc = self.get_val(vector(xc,yc,0))
print "center", xc, yc
scene.center = vector(xc, yc, zc)
def draw_boundary(self):
bound = curve(color=color.blue)
bound.append(pos=(self.xmin, self.ymax,self.data[0,0]))
bound.append(pos=(self.xmin+self.dx*self.ncols, self.ymax,self.data[0,0]))
bound.append(pos=(self.xmin+self.dx*self.ncols,
self.ymax-self.dx*self.nrows,self.data[0,0]))
bound.append(pos=(self.xmin,
self.ymax-self.dx*self.nrows,self.data[0,0]))
bound.append(pos=(self.xmin, self.ymax,self.data[0,0]))
def draw_post(self,posz):
rod = cylinder(pos=vector(posz.x,posz.y,posz.z*self.scale), axis=(0,0,1.5*self.scale), radius=0.2*self.scale)
def draw_curve(self,posx):
self.top_curve= curve(color=color.red, pos=vector(posx.x,posx.y,posx.z*self.scale), radius=0.3*self.scale)
def append_curve(self,posx):
self.top_curve.append(pos= vector(posx.x,posx.y,posx.z*self.scale))
def distance(self, point1, point2):
sum_sq = 0 #summ of squares of the differences
if len(point1) != len(point2):
#check for errors here
raise ValueError("Points have different number of dimensions")
for i in range(len(point1)):
x,y = point1[i], point2[i] #extract the next two ordinates
sum_sq = sum_sq + (y-x)**2 #and add the square of the difference
return sqrt(sum_sq)
def mid_dis (self, point1,point2):
mid= (self.distance (point1,point2))/2
return mid
def midpoint (self, point1, point2):
#FOR MIDPOINT COORDINATES
x_coord = (point1.x + point2.x)/2
y_coord = (point1.y + point2.y)/2
z_coord = (point1.z + point2.z)/2
self.mid_pos = vector (x_coord, y_coord, z_coord)
return (self.mid_pos)
def line_3d(self, imin=0, imax=-99999, jmin=0, jmax=-99999, scale=1.0, center=1):
''' draw a specific region of the raster map as a wire mesh
imin - minimum column value to be drawn
imax - maximum column value
jmin - miinimum row
jmax - maximum row
scale - to add a vertical exageration
center - to center scene on area drawn if center == 1'''
self.scale = scale
default_val = -99999
if imax == default_val:
imax = self.nrows
if jmax == default_val:
jmax = self.ncols
print
tdata = self.data[imin:imax,jmin:jmax]
self.imin=imin
self.imax=imax
self.jmin=jmin
self.jmax=jmax
for i in range(imin, imax):
lin = curve(color=(0,1,0))
for j in range (jmin, jmax):
ipos = self.xmin + j* self.dx
jpos = self.ymax - i* self.dx
kpos = self.data[i,j]
lin.append(pos=(ipos,jpos,kpos*scale))
#print lin.pos[-1]
for j in range (jmin, jmax):
lin = curve(color=(0,1,0))
for i in range(imin, imax):
ipos = self.xmin + j* self.dx
jpos = self.ymax - i* self.dx
kpos = self.data[i,j]
lin.append(pos=(ipos,jpos,kpos*scale))
#center of scene of area
#print self.xmin + jmin* self.dx, self.ymax - imax* self.dx
if center == 1:
xc = self.xmin + (self.dx * ((jmin+jmax)/2))
yc = self.ymax - (self.dx * ((imin+imax)/2))
zc = self.get_val(vector(xc,yc,0))*self.scale
print "center", xc, yc, zc
scene.center = vector(xc, yc, zc)
def local_3x3(self, pos):
#find the local 3x3 matrix of topography around the cell centered at pos
r, c = self.get_rc(pos)
#print "r,c=", r, c
return self.data[r-1:r+2, c-1:c+2]
def get_slope(self, pos, dir):
#the slope given the positon (pos) and direction and the adjacent cells
'''WARNING: THIS IS A PROTOTYPE VERSION AND SHOULD BE USED WITH CAUTION'''
print "WARNING: THIS IS A PROTOTYPE VERSION (sub_raster_array) AND SHOULD BE USED WITH CAUTION"
posz = self.get_val(pos)
pnp = pos + (dir*self.dx)
pnelev = self.get_val(pnp)
slope = (pnelev - posz/self.scale)/(self.dx)
#print "r,c=", r, c
return slope
def sub_raster_array(self, center_pos, cell_range=1):
'''to create a sub raster using only cells within 1 cell_range of the
center cell'''
r, c = self.get_rc(center_pos)
xmin = self.get_x(c - cell_range)
ymax = self.get_y(r - cell_range)
#print center_pos
#print "x, y", self.get_x(c), self.get_y(r)
return raster_map(self.local_3x3(center_pos), xmin, ymax, self.dx)
def get_highest(self, pos):
#local_topo=zeros((self.nrows,self.ncols))
#print local_topo
imax= self.imax
jmax= self.jmax
high_num=0
for i in range(0, self.nrows):
for j in range(0,self.ncols):
val= self.data[i,j]
if val > high_num :
high_num=val
#print 'the highest elevation (so far) is' ,high_num, ' and it came from row' ,i, 'and column' ,j
high_r = i
high_c = j
self.high_r = high_r
self.high_c = high_c
return high_r, high_c, high_num
def get_lowest(self, pos,bound_row):
#local_topo=zeros((self.nrows,self.ncols))
#print local_topo
imax= self.imax
jmax= self.jmax
low_num= 9999
for i in range(0, bound_row):
for j in range(0,self.ncols):
val= self.data[i,j]
if val < low_num :
low_num=val
#print 'the lowest elevation (so far) is' ,low_num, ' and it came from row' ,i, 'and column' ,j
low_r = i
low_c = j
self.low_r = low_r
self.low_c = low_c
return low_r, low_c, low_num
def raster_import(rasterfile):
infile = open(rasterfile,"r")
ct = 0
line = infile.readline()
while line:
ct += 1
l = split(strip(line))
if ct == 1:
ncols = int(l[1])
if ct == 2:
nrows = int(l[1])
if ct == 3:
xmin = float(l[1])
if ct == 4:
ymin = float(l[1])
if ct == 5:
dx = float(l[1])
if ct == 7:
nodata = l[1]
if ct >= 7:
if ct == 7:
topo = zeros((nrows*ncols,), Float)
t_ct = -1
for i in l:
t_ct += 1
if i == nodata: #set no-data values to zero
topo[t_ct] = 0.0
else:
topo[t_ct] = float(i)
line = infile.readline()
infile.close
ymax = ymin + (dx * (nrows-1)) # top left corner
'''datarr is a 2d array '''
#self.xmin = xmin
#self.ymax = ymax
#self.dx = dx
#self.nrows = nrows
#self.ncols = ncols
#self.data = resize(topo, (nrows, ncols))
return raster_map(resize(topo, (nrows, ncols)), xmin, ymax, dx)
def raster_import_ryan(rasterfile, dx=1.0):
infile = open(rasterfile,"r")
old_long = []
old_lat = []
old_data = []
line = infile.readline()
l = splitfields(strip(line),",")
old_long.append(float(l[0]))
old_lat.append(float(l[1]))
old_data.append(float(l[2]))
ct = 0
#data[ct] = float(l[2])
ncols = 1
nrows = 1
while line:
l = splitfields(strip(line),",")
#print float(l[1]), float(l[2])
old_long.append(float(l[0]))
old_lat.append(float(l[1]))
old_data.append(float(l[2]))
#data[ct+1] = float(l[2])
line = infile.readline()
ct += 1
if (old_long[-1] < old_long[-2]):
nrows += 1
ncols = ct/nrows
data = zeros((nrows*ncols,), Float)
print "length = ", len(old_data), nrows*ncols
for i in range (nrows*ncols):
data[i] = old_data[i]
print "rows, columns, number of items ", nrows, ncols, ct
xmin = old_long[0]
ymax = old_lat[0]
'''self.xmin = xmin
self.ymax = ymax
self.dx = dx
self.nrows = nrows
self.ncols = ncols
self.data = resize(data, (nrows, ncols))'''
return raster_map(resize(data, (nrows,ncols)), xmin, ymax, dx)
