Part 1, importing the DEM
From GeoMod
[edit]
Importing the DEM
After trial and tribulation, I have managed to get my DEM imported. It is a section of Shelby County surrounding the Loosahatchie River. It is done on a 30 meter mesh and the file was HUGE. When plotted by python, it first appears to be a large green blur unitl zoomed in and revealing high res dem model.
I then added the sample points for my chemical data and placed them on the map. Next was added a line conecting the points, to profile the river location. This was elevated to make it easily visible. The color changes along this line based on section of stream as defined in the next section.
Here is my code for that section:
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 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))
print "center", xc, yc
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 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)
'''THIS FUNCTION READS IN X, Y, DATA text
Code to add if using Ryans ASCII conversion (reminder: dx is defaulted to 1)'''
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)," ")
#print l
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:
#print ct
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]
xmax = old_long[-1]
ymin = old_lat[-1]
ymax = old_lat[0]
print xmin, xmax,"Delta x = ", xmax-xmin
print ymin, ymax,"Delta y = ", ymax-ymin
pos1 = array([35.311139,-89.639073],Float)#hwy 70
pos2 = array([35.306869,-89.668077],Float)# Hwy 205
pos3 = array([35.281040,-89.766085],Float)# New Brunswick
pos4 = array([35.275161,-89.887326],Float)#204 Singleton Pkwy
pos5 = array([35.263172,-89.929469],Float)#Raliegh Millington
pos6 = array([35.259577,-89.993882],Float)#Hwy 51
pos7 = array([35.254232,-90.026146],Float)#N. Watkins
pt1 = sphere(pos=(xmin+((90.10764+pos1[1])/.000255)*dx,ymax-((35.31319-pos1[0])/.000255)*dx,100),radius=10,color=color.red)
pt2 = sphere(pos=(xmin+((90.10764+pos2[1])/.000255)*dx,ymax-((35.31319-pos2[0])/.000255)*dx,100),radius=10,color=color.red)
pt3 = sphere(pos=(xmin+((90.10764+pos3[1])/.000255)*dx,ymax-((35.31319-pos3[0])/.000255)*dx,100),radius=10,color=color.red)
pt4 = sphere(pos=(xmin+((90.10764+pos4[1])/.000255)*dx,ymax-((35.31319-pos4[0])/.000255)*dx,100),radius=10,color=color.red)
pt5 = sphere(pos=(xmin+((90.10764+pos5[1])/.000255)*dx,ymax-((35.31319-pos5[0])/.000255)*dx,100),radius=10,color=color.red)
pt6 = sphere(pos=(xmin+((90.10764+pos6[1])/.000255)*dx,ymax-((35.31319-pos6[0])/.000255)*dx,100),radius=10,color=color.red)
pt7 = sphere(pos=(xmin+((90.10764+pos7[1])/.000255)*dx,ymax-((35.31319-pos7[0])/.000255)*dx,100),radius=10,color=color.red)
river = curve(color=color.red)
river.append(pos=(xmin+((90.10764+pos1[1])/.000255)*dx,ymax-((35.31319-pos1[0])/.000255)*dx,300),color=color.blue)
river.append(pos=(xmin+((90.10764+pos2[1])/.000255)*dx,ymax-((35.31319-pos2[0])/.000255)*dx,300),color=color.blue)
river.append(pos=(xmin+((90.10764+pos3[1])/.000255)*dx,ymax-((35.31319-pos3[0])/.000255)*dx,300),color=color.red)
river.append(pos=(xmin+((90.10764+pos4[1])/.000255)*dx,ymax-((35.31319-pos4[0])/.000255)*dx,300),color=color.red)
river.append(pos=(xmin+((90.10764+pos5[1])/.000255)*dx,ymax-((35.31319-pos5[0])/.000255)*dx,300),color=color.red)
river.append(pos=(xmin+((90.10764+pos6[1])/.000255)*dx,ymax-((35.31319-pos6[0])/.000255)*dx,300),color=color.green)
river.append(pos=(xmin+((90.10764+pos7[1])/.000255)*dx,ymax-((35.31319-pos7[0])/.000255)*dx,300),color=color.green)
scene.center = vector(1000,1000,-1000)
print xmin+(0.468567/.000255)
return raster_map(resize(data, (nrows,ncols)), xmin, ymax, dx)
#pos1 = array([35.311139,-89.639073],Float)#hwy 70
#newbyy = -1*(35.31319 - 35.311139)*100000
#newbyx = -1*(-90.10764 - (-89.639073))*100000
topog = raster_import_ryan("dantopo30_2.txt")
#topog.line_3d(scale = 1.0)
topog.draw_boundary()
