## LSDMappingTools.py
##=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
## These functions are tools to deal with rasters
##=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
## SMM
## 26/07/2014
##=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
import osgeo.gdal as gdal
import numpy as np
import numpy.ma as ma
from osgeo import osr
from os.path import exists
from osgeo.gdalconst import GA_ReadOnly
from numpy import uint8
from matplotlib import rcParams
#==============================================================================
def GetUTMMaxMin(FileName):
if exists(FileName) is False:
raise Exception('[Errno 2] No such file or directory: \'' + FileName + '\'')
NDV, xsize, ysize, GeoT, Projection, DataType = GetGeoInfo(FileName)
CellSize = GeoT[1]
XMin = GeoT[0]
XMax = XMin+CellSize*xsize
YMax = GeoT[3]
YMin = YMax-CellSize*ysize
return CellSize,XMin,XMax,YMin,YMax
#==============================================================================
#==============================================================================
# this takes rows and columns of minium and maximum values and converts them
# to UTM
def GetUTMMaxMinFromRowsCol(FileName,x_max_col,x_min_col,y_max_row,y_min_row):
if exists(FileName) is False:
raise Exception('[Errno 2] No such file or directory: \'' + FileName + '\'')
NDV, xsize, ysize, GeoT, Projection, DataType = GetGeoInfo(FileName)
CellSize = GeoT[1]
XMin = GeoT[0]
YMax = GeoT[3]
YMin = YMax-CellSize*ysize
xmax_UTM = XMin+x_max_col*CellSize
xmin_UTM = XMin+x_min_col*CellSize
# need to be careful with the ymax_UTM since the rows go from the top
# but the header index is to bottom corner
print "yll: "+str(YMin)+" and nrows: " +str(ysize) + " dx: "+str(CellSize)
ymax_from_bottom = ysize-y_min_row
ymin_from_bottom = ysize-y_max_row
ymax_UTM = YMin+ymax_from_bottom*CellSize
ymin_UTM = YMin+ymin_from_bottom*CellSize
return xmax_UTM,xmin_UTM,ymax_UTM,ymin_UTM
#==============================================================================
#==============================================================================
# This gets the extent of the raster
def GetRasterExtent(FileName):
CellSize,XMin,XMax,YMin,YMax = GetUTMMaxMin(FileName)
extent = [XMin,XMax,YMin,YMax]
return extent
#==============================================================================
# Function to read the original file's projection:
def GetGeoInfo(FileName):
if exists(FileName) is False:
raise Exception('[Errno 2] No such file or directory: \'' + FileName + '\'')
SourceDS = gdal.Open(FileName, gdal.GA_ReadOnly)
if SourceDS == None:
raise Exception("Unable to read the data file")
NDV = SourceDS.GetRasterBand(1).GetNoDataValue()
xsize = SourceDS.RasterXSize
ysize = SourceDS.RasterYSize
GeoT = SourceDS.GetGeoTransform()
Projection = osr.SpatialReference()
Projection.ImportFromWkt(SourceDS.GetProjectionRef())
DataType = SourceDS.GetRasterBand(1).DataType
DataType = gdal.GetDataTypeName(DataType)
return NDV, xsize, ysize, GeoT, Projection, DataType
#==============================================================================
#==============================================================================
# Function to read the original file's projection:
def GetNPixelsInRaster(FileName):
NDV, xsize, ysize, GeoT, Projection, DataType = GetGeoInfo(FileName)
return xsize*ysize
#==============================================================================
#==============================================================================
def ReadRasterArrayBlocks(raster_file,raster_band=1):
if exists(raster_file) is False:
raise Exception('[Errno 2] No such file or directory: \'' + raster_file + '\'')
dataset = gdal.Open(raster_file, GA_ReadOnly )
if dataset == None:
raise Exception("Unable to read the data file")
band = dataset.GetRasterBand(raster_band)
block_sizes = band.GetBlockSize()
x_block_size = block_sizes[0]
y_block_size = block_sizes[1]
#If the block y size is 1, as in a GeoTIFF image, the gradient can't be calculated,
#so more than one block is used. In this case, using8 lines gives a similar
#result as taking the whole array.
if y_block_size < 8:
y_block_size = 8
xsize = band.XSize
ysize = band.YSize
print "xsize: " +str(xsize)+" and y size: " + str(ysize)
max_value = band.GetMaximum()
min_value = band.GetMinimum()
# now initiate the array
data_array = np.zeros((ysize,xsize))
#print "data shape is: "
#print data_array.shape
if max_value == None or min_value == None:
stats = band.GetStatistics(0, 1)
max_value = stats[1]
min_value = stats[0]
for i in range(0, ysize, y_block_size):
if i + y_block_size < ysize:
rows = y_block_size
else:
rows = ysize - i
for j in range(0, xsize, x_block_size):
if j + x_block_size < xsize:
cols = x_block_size
else:
cols = xsize - j
# get the values for this block
values = band.ReadAsArray(j, i, cols, rows)
# move these values to the data array
data_array[i:i+rows,j:j+cols] = values
return data_array
#==============================================================================
#==============================================================================
# Formats ticks for an imshow plot in UTM
# Filename is the name of the file with full path
# x_max, x_min, y_max, y_min are the extent of the plotting area (NOT the DEM)
# n_target ticks are the number of ticks for plotting
#------------------------------------------------------------------------------
def GetTicksForUTM(FileName,x_max,x_min,y_max,y_min,n_target_tics):
CellSize,XMin,XMax,YMin,YMax = GetUTMMaxMin(FileName)
NDV, xsize, ysize, GeoT, Projection, DataType = GetGeoInfo(FileName)
xmax_UTM = XMax
xmin_UTM = XMin
ymax_UTM = YMax
ymin_UTM = YMin
print "now UTM, xmax: " +str(xmax_UTM)+" x_min: " +str(xmin_UTM)+" y_maxb: " +str(ymax_UTM)+" y_minb: " +str(ymin_UTM)
dy_fig = ymax_UTM-ymin_UTM
dx_fig = xmax_UTM-xmin_UTM
dx_spacing = dx_fig/n_target_tics
dy_spacing = dy_fig/n_target_tics
if (dx_spacing>dy_spacing):
dy_spacing = dx_spacing
str_dy = str(dy_spacing)
str_dy = str_dy.split('.')[0]
n_digits = str_dy.__len__()
nd = int(n_digits)
first_digit = float(str_dy[0])
print "str_dy: " +str_dy+ " n_digits: " +str(nd)+" first_digit: " + str(first_digit)
dy_spacing_rounded = first_digit*pow(10,(nd-1))
print "n_digits: "+str(n_digits)+" dy_spacing: " +str(dy_spacing) + " and rounded: "+str(dy_spacing_rounded)
str_xmin = str(xmin_UTM)
str_ymin = str(ymin_UTM)
print "before split str_xmin: "+ str_xmin + " str ymin: " + str_ymin
str_xmin = str_xmin.split('.')[0]
str_ymin = str_ymin.split('.')[0]
print "after split str_xmin: "+ str_xmin + " str ymin: " + str_ymin
xmin_UTM = float(str_xmin)
ymin_UTM = float(str_ymin)
print "UTM: "+ str(xmin_UTM) + " str ymin: " + str(ymin_UTM)
n_digx = str_xmin.__len__()
n_digy = str_ymin.__len__()
print "n_dig_x: " + str(n_digx)+ " nd: " + str(nd)
if (n_digx-nd+1) >= 1:
front_x = str_xmin[:(n_digx-nd+1)]
else:
front_x = str_xmin
if (n_digy-nd+1) >= 1:
front_y = str_ymin[:(n_digy-nd+1)]
else:
front_y = str_ymin
print "xmin: " + str_xmin + " ymin: " + str_ymin + " n_digx: " + str(n_digx)+ " n_digy: " + str(n_digy)
print "frontx: " +front_x+" and fronty: "+ front_y
round_xmin = float(front_x)*pow(10,nd-1)
round_ymin = float(front_y)*pow(10,nd-1)
print "x_min: " +str(xmin_UTM)+ " round xmin: " +str(round_xmin)+ " y_min: " +str(ymin_UTM)+" round y_min: " + str(round_ymin)
# now we need to figure out where the xllocs and ylocs are
xUTMlocs = np.zeros(2*n_target_tics)
yUTMlocs = np.zeros(2*n_target_tics)
xlocs = np.zeros(2*n_target_tics)
ylocs = np.zeros(2*n_target_tics)
new_x_labels = []
new_y_labels = []
for i in range(0,2*n_target_tics):
xUTMlocs[i] = round_xmin+(i)*dy_spacing_rounded
yUTMlocs[i] = round_ymin+(i)*dy_spacing_rounded
#xlocs[i] = (xUTMlocs[i]-XMin)
xlocs[i] = xUTMlocs[i]
# need to account for the rows starting at the upper boundary
ylocs[i] = YMax-(yUTMlocs[i]-YMin)
new_x_labels.append( str(xUTMlocs[i]).split(".")[0] )
new_y_labels.append( str(yUTMlocs[i]).split(".")[0] )
print xUTMlocs
print xlocs
print yUTMlocs
print ylocs
print new_x_labels
print new_y_labels
return xlocs,ylocs,new_x_labels,new_y_labels
#==============================================================================
#==============================================================================
def LogStretchDensityPlot(FileName, thiscmap='gray',colorbarlabel='Elevation in meters',clim_val = (0,0)):
import matplotlib.pyplot as plt
import matplotlib.lines as mpllines
label_size = 20
#title_size = 30
axis_size = 28
# Set up fonts for plots
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Liberation Sans']
rcParams['font.size'] = label_size
# get the data
raster = ReadRasterArrayBlocks(FileName)
# get the log of the raster
raster = np.log10(raster)
# now get the extent
extent_raster = GetRasterExtent(FileName)
x_min = extent_raster[0]
x_max = extent_raster[1]
y_min = extent_raster[2]
y_max = extent_raster[3]
# make a figure, sized for a ppt slide
fig = plt.figure(1, facecolor='white',figsize=(10,7.5))
# make room for the colorbar
fig.subplots_adjust(bottom=0.2)
fig.subplots_adjust(top=0.9)
#fig.subplots_adjust(left=0.2)
#fig.subplots_adjust(right=0.8)
ax1 = fig.add_subplot(1,1,1)
im = ax1.imshow(raster, thiscmap, extent = extent_raster)
print "The is the extent raster data element"
print extent_raster
print "now I am in the mapping routine"
print "x_min: " + str(x_min)
print "x_max: " + str(x_max)
print "y_min: " + str(y_min)
print "y_max: " + str(y_max)
# now get the tick marks
n_target_tics = 5
xlocs,ylocs,new_x_labels,new_y_labels = GetTicksForUTM(FileName,x_max,x_min,y_max,y_min,n_target_tics)
plt.xticks(xlocs, new_x_labels, rotation=60) #[1:-1] skips ticks where we have no data
plt.yticks(ylocs, new_y_labels)
print "The x locs are: "
print xlocs
print "The x labels are: "
print new_x_labels
# some formatting to make some of the ticks point outward
for line in ax1.get_xticklines():
line.set_marker(mpllines.TICKDOWN)
#line.set_markeredgewidth(3)
for line in ax1.get_yticklines():
line.set_marker(mpllines.TICKLEFT)
#line.set_markeredgewidth(3)
plt.xlim(x_min,x_max)
plt.ylim(y_max,y_min)
plt.xlabel('Easting (m)',fontsize = axis_size)
plt.ylabel('Northing (m)', fontsize = axis_size)
ax1.set_xlabel("Easting (m)")
ax1.set_ylabel("Northing (m)")
# set the colour limits
print "Setting colour limits to "+str(clim_val[0])+" and "+str(clim_val[1])
if (clim_val == (0,0)):
print "I don't think I should be here"
im.set_clim(0, np.max(raster))
else:
print "Now setting colour limits to "+str(clim_val[0])+" and "+str(clim_val[1])
im.set_clim(clim_val[0],clim_val[1])
cbar = fig.colorbar(im, orientation='horizontal')
cbar.set_label(colorbarlabel)
#plt.tight_layout()
plt.show()
#==============================================================================
#==============================================================================
def BasicDensityPlot(FileName, thiscmap='gray',colorbarlabel='Elevation in meters',clim_val = (0,0)):
import matplotlib.pyplot as plt
import matplotlib.lines as mpllines
label_size = 20
#title_size = 30
axis_size = 28
# Set up fonts for plots
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Liberation Sans']
rcParams['font.size'] = label_size
# get the data
raster = ReadRasterArrayBlocks(FileName)
# now get the extent
extent_raster = GetRasterExtent(FileName)
x_min = extent_raster[0]
x_max = extent_raster[1]
y_min = extent_raster[2]
y_max = extent_raster[3]
# make a figure, sized for a ppt slide
fig = plt.figure(1, facecolor='white',figsize=(10,7.5))
# make room for the colorbar
#fig.subplots_adjust(bottom=0.1)
#fig.subplots_adjust(top=0.9)
#fig.subplots_adjust(left=0.2)
#fig.subplots_adjust(right=0.8)
ax1 = fig.add_subplot(1,1,1)
im = ax1.imshow(raster, thiscmap, extent = extent_raster)
print "The is the extent raster data element"
print extent_raster
print "now I am in the mapping routine"
print "x_min: " + str(x_min)
print "x_max: " + str(x_max)
print "y_min: " + str(y_min)
print "y_max: " + str(y_max)
# now get the tick marks
n_target_tics = 5
xlocs,ylocs,new_x_labels,new_y_labels = GetTicksForUTM(FileName,x_max,x_min,y_max,y_min,n_target_tics)
plt.xticks(xlocs, new_x_labels, rotation=60) #[1:-1] skips ticks where we have no data
plt.yticks(ylocs, new_y_labels)
print "The x locs are: "
print xlocs
print "The x labels are: "
print new_x_labels
# some formatting to make some of the ticks point outward
for line in ax1.get_xticklines():
line.set_marker(mpllines.TICKDOWN)
#line.set_markeredgewidth(3)
for line in ax1.get_yticklines():
line.set_marker(mpllines.TICKLEFT)
#line.set_markeredgewidth(3)
plt.xlim(x_min,x_max)
plt.ylim(y_max,y_min)
plt.xlabel('Easting (m)',fontsize = axis_size)
plt.ylabel('Northing (m)', fontsize = axis_size)
ax1.set_xlabel("Easting (m)")
ax1.set_ylabel("Northing (m)")
# set the colour limits
print "Setting colour limits to "+str(clim_val[0])+" and "+str(clim_val[1])
if (clim_val == (0,0)):
print "I don't think I should be here"
im.set_clim(0, np.max(raster))
else:
print "Now setting colour limits to "+str(clim_val[0])+" and "+str(clim_val[1])
im.set_clim(clim_val[0],clim_val[1])
cbar = fig.colorbar(im, orientation='vertical')
cbar.set_label(colorbarlabel)
#plt.tight_layout()
plt.show()
#==============================================================================
#==============================================================================
# My attempt to drape two plots
def DrapedPlot(FileName,DrapeFilename):
import matplotlib.pyplot as plt
import matplotlib.lines as mpllines
label_size = 20
#title_size = 30
axis_size = 28
# Set up fonts for plots
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Liberation Sans']
rcParams['font.size'] = label_size
# get the data
raster = ReadRasterArrayBlocks(FileName)
DrapeRaster = ReadRasterArrayBlocks(DrapeFilename)
# now get the extent
extent_raster = GetRasterExtent(FileName)
x_min = extent_raster[0]
x_max = extent_raster[1]
y_min = extent_raster[2]
y_max = extent_raster[3]
# make a figure, sized for a ppt slide
fig = plt.figure(1, facecolor='white',figsize=(10,7.5))
# make room for the colorbar
#fig.subplots_adjust(bottom=0.1)
#fig.subplots_adjust(top=0.9)
ax1 = fig.add_subplot(1,1,1)
im = ax1.imshow(raster, cmap='gray', extent = extent_raster)
# now get the tick marks
n_target_tics = 5
xlocs,ylocs,new_x_labels,new_y_labels = GetTicksForUTM(FileName,x_max,x_min,y_max,y_min,n_target_tics)
plt.xticks(xlocs, new_x_labels, rotation=60) #[1:-1] skips ticks where we have no data
plt.yticks(ylocs, new_y_labels)
# some formatting to make some of the ticks point outward
for line in ax1.get_xticklines():
line.set_marker(mpllines.TICKDOWN)
#line.set_markeredgewidth(3)
for line in ax1.get_yticklines():
line.set_marker(mpllines.TICKLEFT)
#line.set_markeredgewidth(3)
plt.xlim(x_min,x_max)
plt.ylim(y_max,y_min)
plt.xlabel('Easting (m)',fontsize = axis_size)
plt.ylabel('Northing (m)', fontsize = axis_size)
ax1.set_xlabel("Easting (m)")
ax1.set_ylabel("Northing (m)")
#plt.hold(True)
im2 = ax1.imshow(DrapeRaster, cmap=plt.cm.jet, alpha=.6, extent = extent_raster)
#im.set_clim(0, np.max(raster))
#cbar = fig.colorbar(im, orientation='vertical')
#cbar.set_label("Elevation in meters")
#plt.tight_layout()
plt.show()
#==============================================================================
#==============================================================================
# Make a simple hillshade plot
def Hillshade(raster_file, azimuth, angle_altitude):
array = ReadRasterArrayBlocks(raster_file,raster_band=1)
x, y = np.gradient(array)
slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))
aspect = np.arctan2(-x, y)
azimuthrad = np.azimuth*np.pi / 180.
altituderad = np.angle_altitude*np.pi / 180.
shaded = np.sin(altituderad) * np.sin(slope)\
+ np.cos(altituderad) * np.cos(slope)\
* np.cos(azimuthrad - aspect)
return 255*(shaded + 1)/2
#==============================================================================
def round_to_n(x, n):
if n < 1:
raise ValueError("number of significant digits must be >= 1")
# Use %e format to get the n most significant digits, as a string.
format = "%." + str(n-1) + "e"
as_string = format % x
return float(as_string)
def read_headers(input_file):
with open(input_file+'.hdr','r') as f:
return [float(h) if not h.isalpha() else h for h in [l.split()[1] for l in f.readlines()]] #isdigit() does not catch floats
def read_bin(filename):
import sys
import numpy as np
with open(filename + '.flt', "rb") as f:
raster_data = np.fromstring(f.read(), 'f')
if sys.byteorder == 'big':
raster_data = raster_data.byteswap() #ensures data is little endian
return raster_data
def read_flt(input_file):
if input_file.endswith('.flt') or input_file.endswith('.hdr'):
input_file = input_file[:-4]
else:
print 'Incorrect filename'
return 0,0 #exits module gracefully
headers = read_headers(input_file)
#read the data as a 1D array and reshape it to the dimensions in the header
raster_array = read_bin(input_file).reshape(headers[1], headers[0])
raster_array = raster_array.reshape(headers[1], headers[0]) #rows, columns
return raster_array, headers
def read_ascii_raster(ascii_raster_file):
import numpy as np
with open(ascii_raster_file) as f:
header_data = [float(f.next().split()[1]) for x in xrange(6)] #read the first 6 lines
raster_data = np.genfromtxt(ascii_raster_file, delimiter=' ', skip_header=6)
raster_data = raster_data.reshape(header_data[1], header_data[0]) #rows, columns
return raster_data, header_data
# this gets the extent of the asc for use with plotting
# It returns a list with 4 elements, x_min, x_max, y_min,y_max
def get_raster_extent_asc(header):
x_min = header[2]
y_min = header[3]
spacing = header[4]
n_cols = header[0]
n_rows = header[1]
x_max = x_min+n_cols*spacing
y_max = y_min+n_rows*spacing
extent = [x_min,x_max,y_min,y_max]
return extent
# This function makes a simple density plot of a raster
def simple_density_plot_asc(rfname):
import numpy as np, matplotlib.pyplot as plt
from matplotlib import rcParams
import matplotlib.colors as colors
import matplotlib.cm as cmx
label_size = 20
#title_size = 30
axis_size = 28
# Set up fonts for plots
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Liberation Sans']
rcParams['font.size'] = label_size
# get the data
raster,header = read_ascii_raster(rfname)
# now get the extent
extent_raster = get_raster_extent_asc(header)
#print extent_raster
# make a figure, sized for a ppt slide
fig = plt.figure(1, facecolor='white',figsize=(10,7.5))
ax1 = fig.add_subplot(1,1,1)
im = ax1.imshow(raster, cmap='gray', extent = extent_raster)
ax1.set_xlabel("Easting (m)")
ax1.set_ylabel("Northing (m)")
im.set_clim(0, np.max(raster))
cbar = fig.colorbar(im, orientation='horizontal')
cbar.set_label("Elevation in meters")
plt.show()
def format_ticks_for_UTM_imshow(hillshade_header,x_max,x_min,y_max,y_min,n_target_tics):
import numpy as np
xmax_UTM = hillshade_header[2]+x_max*hillshade_header[4]
xmin_UTM = hillshade_header[2]+x_min*hillshade_header[4]
# need to be careful with the ymax_UTM since the rows go from the top
# but the header index is to bottom corner
print "yll: "+str(hillshade_header[3])+" and nrows: " +str(hillshade_header[1]) + " dx: "+str(hillshade_header[4])
ymax_from_bottom = hillshade_header[1]-y_min
ymin_from_bottom = hillshade_header[1]-y_max
ymax_UTM = hillshade_header[3]+ymax_from_bottom*hillshade_header[4]
ymin_UTM = hillshade_header[3]+ymin_from_bottom*hillshade_header[4]
print "now UTM, xmax: " +str(xmax_UTM)+" x_min: " +str(xmin_UTM)+" y_maxb: " +str(ymax_UTM)+" y_minb: " +str(ymin_UTM)
dy_fig = ymax_UTM-ymin_UTM
dx_fig = xmax_UTM-xmin_UTM
dx_spacing = dx_fig/n_target_tics
dy_spacing = dy_fig/n_target_tics
if (dx_spacing>dy_spacing):
dy_spacing = dx_spacing
str_dy = str(dy_spacing)
str_dy = str_dy.split('.')[0]
n_digits = str_dy.__len__()
nd = int(n_digits)
first_digit = float(str_dy[0])
print "str_dy: " +str_dy+ " n_digits: " +str(nd)+" first_digit: " + str(first_digit)
dy_spacing_rounded = first_digit*pow(10,(nd-1))
print "n_digits: "+str(n_digits)+" dy_spacing: " +str(dy_spacing) + " and rounded: "+str(dy_spacing_rounded)
str_xmin = str(xmin_UTM)
str_ymin = str(ymin_UTM)
str_xmin = str_xmin.split('.')[0]
str_ymin = str_ymin.split('.')[0]
xmin_UTM = float(str_xmin)
ymin_UTM = float(str_ymin)
n_digx = str_xmin.__len__()
n_digy = str_ymin.__len__()
front_x = str_xmin[:(n_digx-nd+1)]
front_y = str_ymin[:(n_digy-nd+1)]
print "xmin: " + str_xmin + " ymin: " + str_ymin + " n_digx: " + str(n_digx)+ " n_digy: " + str(n_digy)
print "frontx: " +front_x+" and fronty: "+ front_y
round_xmin = float(front_x)*pow(10,nd-1)
round_ymin = float(front_y)*pow(10,nd-1)
print "x_min: " +str(xmin_UTM)+ " round xmin: " +str(round_xmin)+ " y_min: " +str(ymin_UTM)+" round y_min: " + str(round_ymin)
# now we need to figure out where the xllocs and ylocs are
xUTMlocs = np.zeros(2*n_target_tics)
yUTMlocs = np.zeros(2*n_target_tics)
xlocs = np.zeros(2*n_target_tics)
ylocs = np.zeros(2*n_target_tics)
new_x_labels = []
new_y_labels = []
for i in range(0,2*n_target_tics):
xUTMlocs[i] = round_xmin+(i)*dy_spacing_rounded
yUTMlocs[i] = round_ymin+(i)*dy_spacing_rounded
xlocs[i] = (xUTMlocs[i]-hillshade_header[2])/hillshade_header[4]
# need to account for the rows starting at the upper boundary
ylocs[i] = hillshade_header[1]-((yUTMlocs[i]-hillshade_header[3])/hillshade_header[4])
new_x_labels.append( str(xUTMlocs[i]).split(".")[0] )
new_y_labels.append( str(yUTMlocs[i]).split(".")[0] )
return xlocs,ylocs,new_x_labels,new_y_labels
def vectorize(hillshade_file, m_value_file):
import matplotlib.pyplot as pp
import numpy as np
import matplotlib.colors as colors
import matplotlib.cm as cmx
from matplotlib import rcParams
#get data
hillshade, hillshade_header = read_flt(hillshade_file)
m_values, m_values_header = read_flt(m_value_file)
#handle plotting hillshades which are larger than the m_value raster
#cannot cope with m_value raster larger than the hillshade
corrected_x = 0
corrected_y = 0
if (hillshade_header[0] != m_values_header[0]) or (hillshade_header[1] != m_values_header[1]):
corrected_x = (m_values_header[2] - hillshade_header[2]) / hillshade_header[4]
corrected_y = (((m_values_header[3] / m_values_header[4] )+ m_values_header[1]) - ((hillshade_header[3] / hillshade_header[4]) + hillshade_header[1])) * -1
#ignore nodata values
hillshade = np.ma.masked_where(hillshade == -9999, hillshade)
m_values = np.ma.masked_where(m_values == -9999, m_values)
#fonts
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Liberation Sans']
rcParams['font.size'] = 12
fig, ax = pp.subplots()
ax.imshow(hillshade, vmin=0, vmax=255, cmap=cmx.gray)
xlocs, xlabels = pp.xticks()
ylocs, ylabels = pp.yticks()
new_x_labels = np.linspace(hillshade_header[2],hillshade_header[2]+(hillshade_header[1]*hillshade_header[4]), len(xlocs))
new_y_labels = np.linspace(hillshade_header[3],hillshade_header[3]+(hillshade_header[0]*hillshade_header[4]), len(ylocs))
new_x_labels = [str(x).split('.')[0] for x in new_x_labels] #get rid of decimal places in axis ticks
new_y_labels = [str(y).split('.')[0] for y in new_y_labels][::-1] #invert y axis
pp.xticks(xlocs[1:-1], new_x_labels[1:-1], rotation=30) #[1:-1] skips ticks where we have no data
pp.yticks(ylocs[1:-1], new_y_labels[1:-1])
pp.xlabel('Easting (m)')
pp.ylabel('Northing (m)')
# SET UP COLOURMAPS
jet = pp.get_cmap('jet')
m_MIN = np.min(m_values)
m_MAX = np.max(m_values)
cNorm_m_values = colors.Normalize(vmin=m_MIN, vmax=m_MAX)
scalarMap_m_values = cmx.ScalarMappable(norm=cNorm_m_values, cmap=jet)
for i in xrange(len(m_values)):
for j in xrange(len(m_values[0])):
if m_values[i][j] > 0:
colorVal = scalarMap_m_values.to_rgba(m_values[i][j])
pp.scatter(j + corrected_x, i + corrected_y, marker=".", color=colorVal,edgecolors='none')
# Configure final plot
sm = pp.cm.ScalarMappable(cmap=jet,norm=pp.normalize(vmin=m_MIN, vmax=m_MAX))
sm._A = []
cbar = pp.colorbar(sm)
cbar.set_label('M Values')
pp.show()
