In [1]:
clf()
from skimage import draw
In [2]:
import gdal
from gdalconst import *
filename = '/Users/Anya/Data/HiRISE_dems/DTEPC_022699_0985_022607_0985_A01.IMG'
dataset = gdal.Open( filename, GA_ReadOnly )
In [3]:
if dataset is None: print 'none'
In [4]:
print 'Driver: ', dataset.GetDriver().ShortName,'/', \
dataset.GetDriver().LongName
print 'Size is ',dataset.RasterXSize,'x',dataset.RasterYSize, \
'x',dataset.RasterCount
print 'Projection is ',dataset.GetProjection()
geotransform = dataset.GetGeoTransform()
if not geotransform is None:
print 'Origin = (',geotransform[0], ',',geotransform[3],')'
print 'Pixel Size = (',geotransform[1], ',',geotransform[5],')'
In [5]:
band = dataset.GetRasterBand(1)
print 'Band Type=',gdal.GetDataTypeName(band.DataType)
print 'DEM size:', band.XSize, 'x', band.YSize
min = band.GetMinimum()
max = band.GetMaximum()
if min is None or max is None:
(min,max) = band.ComputeRasterMinMax(1)
print 'Min=%.3f, Max=%.3f' % (min,max)
if band.GetOverviewCount() > 0:
print 'Band has ', band.GetOverviewCount(), ' overviews.'
if not band.GetRasterColorTable() is None:
print 'Band has a color table with ', \
band.GetRasterColorTable().GetCount(), ' entries.'
In [6]:
x_start = 10
y_start = 3000
x_points = 12000
y_points = 10
In [7]:
data = band.ReadAsArray(x_start, y_start, x_points, y_points)
In [8]:
data[data == band.GetNoDataValue()] = nan
data[data < -6000] = nan
print data.shape
print 'Min=', data.min(), 'Max=', data.max()
In [9]:
clf()
plot(data[5,:], '*')
#imshow(data, cmap='Paired', vmax = 1870)
#colorbar()
Out[9]:
In [10]:
def height_profile(x1,y1,x2,y2):
line_x, line_y = draw.line(x1,y1,x2,y2)
height = zeros(line_x.size)
dist = zeros(line_x.size)
for ind in range(line_x.size):
height[ind] = band.ReadAsArray(int(line_x[ind]), int(line_y[ind]), 1, 1)
dist[ind] = sqrt((line_x[ind] - line_x[0])**2 + (line_y[ind] - line_y[0])**2)
return height, dist
In [11]:
clf()
In [12]:
#RIM = [10940, 4280, 12195, 3290]
#RIM = [8845, 4925, 9555, 4405]
RIM = [5770, 1060, 6640, 2166] # E-W
#RIM = [7101, 3460, 8190, 5280] # E-W
#RIM = [5675, 4143, 7014, 3121]
#RIM = [8079, 2330, 8880,1655]
height, dist = height_profile(RIM[0], RIM[1], RIM[2], RIM[3])
for nr in range(3):
height2, dist2 = height_profile(RIM[0]+nr,RIM[1]+nr,RIM[2]+nr,RIM[3]+nr)
height3, dist3 = height_profile(RIM[0]-nr,RIM[1]-nr,RIM[2]-nr,RIM[3]-nr)
for i in range(dist.shape[0]):
height[i] = ( height[i] + height2[i] + height3[i] )/3
c1 = 'b' #'yellow' #'black' #'blue' #'green'#'red' #'orange'
plot(500+dist, height , 'x', color = c1)
plot(500+dist, height , color = c1)
xlabel('distance, m')
ylabel('height, m')
grid(True)
In [13]:
spiders = band.ReadAsArray(9970, 3230, 300,300)#9723, 3435, 300, 300)
from scipy import ndimage
mean_dem = ndimage.filters.gaussian_filter(spiders, 20, mode='nearest')
spiders2= spiders-mean_dem
In [14]:
clf()
imshow(spiders2[10:-10,10:-10], vmax = 0)#1890)
bar=colorbar()
bar.set_label('height, m')
xlabel('distance, m')
ylabel('distance, m')
rcParams.update({'font.size': 22})
In [15]:
x, y = np.mgrid[0:data.shape[0], 0:data.shape[1]]
#from mayavi import mlab
#s = mlab.mesh(x, y, data)
In [16]:
def gaussianFilter(sizex,sizey=None,scale=0.333):
'''
Generate and return a 2D Gaussian function
of dimensions (sizex,sizey)
If sizey is not set, it defaults to sizex
A scale can be defined to widen the function (default = 0.333)
'''
sizey = sizey or sizex
x, y = np.mgrid[-sizex:sizex+1, -sizey:sizey+1]
g = np.exp(-scale*(x**2/float(sizex)+y**2/float(sizey)))
return g/g.sum()
In [17]:
#alternative aspect calculation
def grad2d(dem):
'''
Calculate the slope and gradient of a DEM
'''
from scipy import signal
f0 = gaussianFilter(3)
I = signal.convolve(dem,f0,mode='valid')
f1 = np.array([[-1,0,1],[-2,0,2],[-1,0,1]])
f2 = f1.transpose()
g1 = signal.convolve(I,f1,mode='valid')
g2 = signal.convolve(I,f2,mode='valid')
slope = sqrt(g1**2 + g2**2)
aspect = np.arctan2(g2,g1)
return slope, aspect
In [18]:
slope, aspect = grad2d(data)
In [19]:
degrees(aspect.max()+pi), degrees(aspect.min()+pi)
Out[19]:
In [20]:
#figure(1)
#clf()
#imshow(degrees(aspect+pi), cmap='gray')
#colorbar()
In [21]:
#open slope map created by gdaldem
slope_file = "/Users/Anya/Data/HiRISE_dems/DTEPC_022699_0985_022607_0985_A01_SLOPE.TIFF"
slope_dataset = gdal.Open( slope_file, GA_ReadOnly )
if slope_dataset is None: print 'none'
slope_band = slope_dataset.GetRasterBand(1)
slope_data = slope_band.ReadAsArray(x_start, y_start, x_points, y_points)
slope_data[slope_data == slope_band.GetNoDataValue()] = nan
print slope_data.shape
print 'Min=', slope_data.min(), 'Max=', slope_data.max()
In [22]:
slope_spiders = slope_band.ReadAsArray(9970, 3230, 300,300)
In [23]:
figure(2)
clf()
imshow((slope_spiders))
colorbar()
Out[23]:
In [24]:
#open aspect map created by gdaldem
aspect_file = "/Users/Anya/Data/HiRISE_dems/DTEPC_022699_0985_022607_0985_A01_ASPECT.TIFF"
aspect_dataset = gdal.Open( aspect_file, GA_ReadOnly )
if aspect_dataset is None: print 'none'
aspect_band = aspect_dataset.GetRasterBand(1)
aspect_data = aspect_band.ReadAsArray(x_start, y_start, x_points, y_points)
aspect_data[aspect_data == aspect_band.GetNoDataValue()] = nan
print aspect_data.shape
print 'Min=', aspect_data.min(), 'Max=', aspect_data.max()
In [25]:
#figure(3)
#clf()
#imshow(aspect_data, cmap='gray')
#colorbar()
In [26]:
cone_dem = zeros(10000).reshape(100,100)
x, y = np.mgrid[0:cone_dem.shape[0], 0:cone_dem.shape[1]]
center_x=20
center_y=50
top = sqrt((center_x - (center_x+1))**2+(center_y - (center_y+1))**2)+1
diam = 10
for x2 in range(center_x-diam, center_x+diam):
for y2 in range(center_y-diam, center_y+diam): cone_dem[x2, y2] = top - sqrt((center_x - x2)**2+(center_y - y2)**2)
In [27]:
#open complete DEM
x = band.XSize
y = band.YSize
#define nr of 100*100 sub-frames
subframe_x_size = 100
subframe_y_size = 100
overlap = 20
nr_x = (x - overlap)//subframe_x_size
nr_y = (y - overlap)//subframe_y_size
ix = range(overlap, x, subframe_x_size)
iy = range(overlap, y, subframe_y_size)
In [46]:
clf()
for i in range(nr_x): #nrx):
for j in range(nr_y): #nry):
out_file = 'DEM' + '_LH_atPointX'+ str(ix[i]) +'Y'+ str(iy[j])+'-size-'+str(subframe_x_size)+'by'+str(subframe_y_size)
if (ix[i]+subframe_x_size > x): # - check f the last frame not getting out of complete DEM array
subframe_x_size_last = int(x - ix[i]-1)
else:
subframe_x_size_last = subframe_x_size
if (iy[j]+subframe_y_size > y):
subframe_y_size_last = int(y - iy[j]-1)
else:
subframe_y_size_last = subframe_y_size
DEM_subframe = band.ReadAsArray(ix[i], iy[j], subframe_x_size_last, subframe_y_size_last)
DEM_subframe[where(DEM_subframe < -10e20)] = nan
if (-isnan(DEM_subframe)).all():
#print DEM_subframe.shape
imshow(DEM_subframe)
colorbar()
break
In [43]:
-isnan(DEM_subframe).any()
Out[43]:
In [34]:
reshape?
In [ ]: