In [1]:
import pandas as pd
df = pd.read_csv('../raw_data/metadata_clean.csv')
lon, lat = df['lon'], df['lat']
In [14]:
from osgeo import gdal
dataset1 = gdal.Open('../maps/grdn40w121_1/w001001.adf')
dataset2 = gdal.Open('../maps/grdn39w121_1/w001001.adf')
arr = np.append(dataset1.ReadAsArray(), dataset2.ReadAsArray(), axis=0,)
geotransform = dataset1.GetGeoTransform()
origin = geotransform[0], geotransform[3]
dlon, dlat = geotransform[1], geotransform[5]
# Sun direction
azimuth=315.0
# Sun angle
altitude=45.0
# Elevation exageration
z=1.0
# Resolution
scale=1.0
xres = dlon*111.045e3 * np.cos(np.deg2rad(39.)) # resolution in m correct for latitude
yres = -dlat*111.045e3 # resolution in m
print xres, yres
deg2rad = 3.141592653589793 / 180.0
# Exclude 2 pixels around the edges which are usually NODATA.
# Also set up structure for a 3x3 window to process the slope
# throughout the grid
window = []
for row in range(3):
for col in range(3):
window.append(arr[row:(row + arr.shape[0] - 2), col:(col + arr.shape[1] - 2)])
# Process each cell
x = ((z*window[0] + z*window[3] + z*window[3] + z*window[6]) - \
(z*window[2] + z*window[5] + z*window[5] + z*window[8])) / (8.0*xres*scale);
y = ((z*window[6] + z*window[7] +z* window[7] + z*window[8])
- (z*window[0] + z*window[1] + z*window[1] + z*window[2])) / (8.0*yres*scale);
# b = window[2] + 2*window[5] + window[8] - window[0] - 2*window[3] - window[6]/(8.0*xres)
# c = window[0] + 2*window[1] + window[2] - window[6] - 2*window[7] - window[8]/(8.0*yres)
# Calculate slope
slope = 90.0 - np.rad2deg(np.arctan(np.sqrt(x*x + y*y)))
#slope = 90 - np.rad2deg(np.arctan(np.sqrt(b**2 + c**2)))
# Calculate aspect
aspect = np.arctan2(x, y)
# Calculate the shaded relief
shaded = np.sin(altitude * deg2rad) * np.sin(slope * deg2rad) \
+ np.cos(altitude * deg2rad) * np.cos(slope * deg2rad) \
* np.cos((azimuth - 90.0) * deg2rad - aspect);
shaded = shaded * 255
In [47]:
def PlotShaded(arr, vmin, vmax, cbar=False):
plt.figure(figsize=(8,8))
extent = [origin[0], origin[0]+dlon*shaded.shape[1], origin[1]+dlat*shaded.shape[0], origin[1]]
#im = plt.imshow(arr,extent=extent,aspect=1, cmap='gray' , vmin=0)
im = plt.imshow(arr,extent=extent,aspect=1, cmap='gray' , vmin=vmin, vmax=vmax,)
if cbar:
plt.colorbar(label='Slope [deg]')
lon, lat = df['lon'], df['lat']
plt.scatter(lon, lat, s=5, c='r', marker='+')
plt.xlim(-120.8, -120 )
plt.ylim(38.5, 39.5)
plt.gca().ticklabel_format(useOffset=False)
plt.ylabel('Latitude')
plt.xlabel('Longitude')
n_sites = df['site_id'].max()
for i_site in range(1, n_sites+1):
idx = df['site_id']==i_site
lats, lons = df['lat'][idx], df['lon'][idx]
mean_lon, mean_lat = np.mean(lons), np.mean(lats)
plt.text(mean_lon+.007, mean_lat+.007, r'{\bf %i}'%i_site, color='r', fontsize=14)
PlotShaded(shaded, vmin=0, vmax=None, cbar=False)
PlotShaded(90-slope, vmin=0, vmax=30, cbar=True)
In [15]:
plt.figure(figsize=(12,10))
n_sites = df['site_id'].max()
for i_site in range(1, n_sites+1):
plt.subplot(5, 3, i_site)
idx = df['site_id']==i_site
lats, lons = df['lat'][idx], df['lon'][idx]
# Read out the metadata.
geotransform = dataset1.GetGeoTransform()
origin = geotransform[0], geotransform[3]
dlon, dlat = geotransform[1], geotransform[5]
# Plot the topographical data
im = plt.imshow(shaded, extent=extent, aspect=1, cmap='gray',)
#plt.colorbar()
lon, lat = lons[idx], lats[idx]
plt.scatter(lon, lat, s=5, c='r', marker='+')
dlon_total = lon.max()-lon.min()
dlat_total = lat.max()-lat.min()
plt.xlim(lon.min()-.4*dlon_total, lon.max()+.4*dlon_total)
plt.ylim(lat.min()-.4*dlat_total, lat.max()+.4*dlat_total)
plt.gca().ticklabel_format(useOffset=False)
n_sites = df['site_id'].max()
idx = np.where(df['site_id'].values==i_site)[0]
node_lats, node_lons = df['lat'][idx].values, df['lon'][idx].values
plt.text(.05, .05, 'Site %i'%(i_site), transform=plt.gca().transAxes, color='r')
for i_node in range(0, len(idx)):
plt.text(node_lons[i_node]+.00005, node_lats[i_node]+.00005, '%i'%(i_node+1), color='r')
In [69]:
from scipy.interpolate import RegularGridInterpolator as RGI
def InterpolateGrid(lon, lat, arr):
# Use bilinear interpolation
lats = np.linspace(origin[1]+dlat*arr.shape[0] , origin[1], arr.shape[0])
lons = np.linspace(origin[0], origin[0]+dlon*arr.shape[1], arr.shape[1])
rgi = RGI( (lats, lons), arr[::-1,] )
return rgi((lat,lon))
# Check that the DEM data is pretty close to the actual elevation
dem_elevation = InterpolateGrid(df['lon'].values, df['lat'].values, arr)
dem_slopes = InterpolateGrid(df['lon'].values, df['lat'].values, 90-slope)
dem_aspect = InterpolateGrid(df['lon'].values, df['lat'].values, np.rad2deg(aspect))
dem_aspect[dem_aspect<0] += 360
# Plot the elevation error compared with the sensor metadata values
plt.plot(dem_elevation - df['elevation'].values, linestyle='', marker='o', markersize=2)
print 'Sum of Absolute Elevation Errors', np.sum(np.abs(dem_elevation - df['elevation'].values))
plt.ylabel('DEM-GPS Elevation [m]')
plt.xlabel('Sensor Number')
running_count = 0
for i_site in range(1, n_sites+1):
idx = np.where(df['site_id']==i_site)[0]
node_count = len(idx)
plt.vlines(running_count, -150,250,alpha=.4)
running_count += node_count
plt.text(running_count-.5*(node_count), -145, '%i'%i_site, horizontalalignment='center', fontsize=8)
In [74]:
df['slope'] = dem_slopes
df['aspect'] = dem_aspect
df.to_csv('../raw_data/metadata_clean_with_slope_aspect.csv')
In [ ]:
In [ ]:
In [ ]:
In [ ]: