In [19]:
class PrettyFloat(float):
def __repr__(self):
return '%.2f' % self
def pretty_floats(obj):
if isinstance(obj, float):
return PrettyFloat(obj)
elif isinstance(obj, dict):
return dict((k, pretty_floats(v)) for k, v in obj.items())
elif isinstance(obj, (list, tuple)):
return map(pretty_floats, obj)
return obj
import json, glob, os
base_url = '/static/html/ExRidges/'
default_name = 'not_available.html'
base_dir = '/home/mchin/workspace/sarah/Revised_Portal_Dec2016/'
active_ridges = base_dir + 'Active_ridges/ActiveRidges.geojson'
active_ridges_analyzed = base_dir + 'Shapefiles/ActiveRidges_analyzed.geojson'
excluded_tier_extinct_ridges = base_dir + 'Shapefiles/ExcludedTier_Extinct_Ridges.geojson'
primary_tier_extinct_ridges = base_dir + 'Shapefiles/PrimaryTier_Extinct_Ridges.geojson'
secondary_tier_extinct_ridges = base_dir + 'Shapefiles/SecondaryTier_Extinct_Ridges.geojson'
RED = [1.0, 0.0, 0.0, 1.0]
PINK = [255.0/255, 192.0/255, 203.0/255, 1.0]
YELLOW = [1.0,1.0,0.0,1.0]
BLUE = [0.0,0.0,1.0,1.0]
DARK_ORANGE = [1.0, 140.0/255, 0]
PURPLE = [160.0/255, 32.0/255.0, 240.0/255.0]
new_data = {"type": "FeatureCollection"}
new_data["features"] = []
cnt = 0
for r in [active_ridges,
active_ridges_analyzed,
excluded_tier_extinct_ridges,
primary_tier_extinct_ridges,
secondary_tier_extinct_ridges]:
with open(r) as f:
data = json.load(f)
for f in data["features"]:
cnt += 1
feature = {"type": "Feature"}
feature["geometry"] = f["geometry"]
feature["properties"] = {}
feature["properties"]["width"] = 3.0
if r == active_ridges:
feature["properties"]["color"] = DARK_ORANGE
feature["properties"]["width"] = 1.0
feature["properties"]["Tier"] = 'None'
feature["properties"]['Ocean'] = 'None'
feature["properties"]['RidgeTypeFlag'] = 'Active'
elif r == active_ridges_analyzed:
feature["properties"]["color"] = PURPLE
try:
feature["properties"]['name'] = f["properties"]['NAME'][0:f["properties"]['NAME'].index('<')]
except:
feature["properties"]['name'] = f["properties"]['NAME']
feature["properties"]['RidgeTypeFlag'] = 'AnalyzedActive'
feature["properties"]['AnalyzedActiveRidge'] = True
elif r in [excluded_tier_extinct_ridges,
primary_tier_extinct_ridges,
secondary_tier_extinct_ridges]:
feature["properties"] = f["properties"]
feature["properties"]["width"] = 3.0
feature["properties"]['RidgeTypeFlag'] = 'Extinct'
if r == excluded_tier_extinct_ridges:
feature["properties"]["color"] = YELLOW
elif r == primary_tier_extinct_ridges:
feature["properties"]["color"] = RED
elif r == secondary_tier_extinct_ridges:
feature["properties"]["color"] = PINK
else:
feature["properties"]["color"] = [0.0,0.0,0.0,1.0]
feature["properties"]["name"] = f["properties"]["SegmentFul"]
del feature["properties"]["SegmentFul"]
#feature["properties"]["ID"] = f["properties"]["SEG_ID"]
feature["properties"]["ID"] = str(cnt)
#del feature["properties"]["SEG_ID"]
full_file_path = ''
if 'htmlPage' in f["properties"]:
full_file_path = base_dir+'ExRidges_HTML_pages/{}'.format(f["properties"]["htmlPage"])
elif 'HTMLPage' in f["properties"]:
full_file_path = base_dir+'ExRidges_HTML_pages/{}'.format(f["properties"]["HTMLPage"])
if os.path.isfile(full_file_path):
feature["properties"]["URL"] = base_url+"ExRidges_HTML_pages/" + os.path.basename(full_file_path)
else:
feature["properties"]["URL"] = None
new_data["features"].append(feature)
for t in new_data["features"]:
if not t['properties']:
print t
with open(base_dir + 'ExRidges_new.geojson', "w") as fo:
fo.write(json.dumps(pretty_floats(new_data)))
print 'done'
In [ ]:
ROOT_PATH = "/mnt/workspace/sarah/Portal_ExRidges"
import sys
sys.path.append(ROOT_PATH)
import split_grid, shading
import struct, math, gzip, os
import numpy
from osgeo import gdal
from gdalconst import *
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.colors import hsv_to_rgb
%matplotlib inline
from pylab import *
count = 10
show_img_flag = False;
In [ ]:
input_file = ROOT_PATH + '/GEBCO_2014_1D.nc'
output_dir = ROOT_PATH + '/GEBCO_2014_1D_pieces/'
split_grid.split_grid(input_file,output_dir,count)
In [ ]:
lon_inc = 360/count
lat_inc = 180/count
#x_idx = [1]
#y_idx = [4]
x_idx = range(count)
y_idx = range(count)
root_dir = ROOT_PATH
output_dir= root_dir + '/results/'
def toRGB(data):
hsv = data.split('-')
return hsv_to_rgb([float(hsv[0])/360.0, float(hsv[1]), float(hsv[2])])
def gmtColormap(filePath):
steps = []
colors = []
with open(filePath,'r') as f:
lines = f.readlines()
for l in lines:
ls = l.split()
if ls[0] == '#' or ls[0] == "B" or ls[0] == "F" or ls[0] == "N":
continue
else:
steps.append(float(ls[0]))
steps.append(float(ls[2]))
colors.append(toRGB(ls[1]))
colors.append(toRGB(ls[3]))
color_list = []
for i in range(len(steps)):
color_list.append((float(steps[i]-steps[0])/(steps[-1]-steps[0]), colors[i]))
#print color_list
return LinearSegmentedColormap.from_list('cmap', color_list, N=1024)
def gmtColormap_1(filePath):
steps = []
colors = []
with open(filePath,'r') as f:
lines = f.readlines()
for l in lines:
ls = l.split()
if ls[0] == '#' or ls[0] == "B" or ls[0] == "F" or ls[0] == "N":
continue
else:
steps.append(float(ls[0]))
steps.append(float(ls[4]))
colors.append(hsv_to_rgb([float(ls[1])/360.0,float(ls[2]),float(ls[3])]))
colors.append(hsv_to_rgb([float(ls[5])/360.0,float(ls[6]),float(ls[7])]))
color_list = []
for i in range(len(steps)):
color_list.append((float(steps[i]-steps[0])/(steps[-1]-steps[0]), colors[i]))
#print color_list
return LinearSegmentedColormap.from_list('cmap', color_list, N=1024)
#cmap = gmtColormap(ROOT_PATH+"/Gebco_BathymetryTopo.cpt")
cmap = gmtColormap_1(ROOT_PATH+"/topo.cpt")
if not os.path.exists(output_dir):
os.makedirs(output_dir)
for i in x_idx:
for j in y_idx:
dataset = gdal.Open(root_dir + '/GEBCO_2014_1D_pieces/{0}_{1}.tif'.format(i,j), GA_ReadOnly )
band = dataset.GetRasterBand(1)
r = band.ReadAsArray( 0, 0, band.XSize, band.YSize, band.XSize, band.YSize)
#r = band.ReadAsArray( 0, 0, band.XSize, band.YSize, 3600, 1800)
fig = plt.figure(figsize=(36, 18),dpi=240, frameon=False)
ax = plt.Axes(fig, [0., 0., 1., 1.])
ax.set_axis_off()
fig.add_axes(ax)
plt.ioff()
'''
You must be careful with the Nan data in the grid.
'''
where_are_NaNs = isnan(r)
r[where_are_NaNs] = 0
rgb = shading.shade(r,cmap=cmap,intensity=shading.intensity(r),vmin=-7000, vmax=7000)
plt.imshow(rgb)
#plt.imshow(r,cmap=cmap,vmin=-10927.5, vmax=8726)
fig.savefig(output_dir+'{0}_{1}.tif'.format(i,j),pad_inches=0,dpi=240)
'''
If the image saved by matplotlib has the annoying frame margin, use the following imagemagic command to get rid of it.
You might need to use "transparent=True," when saving the image with matplotlib.
'''
#cmd = "convert "+output_dir+'{0}_{1}_topo15.tif '.format(i,j) +" -trim "+ output_dir+'{0}_{1}_topo15_trimmed.tif'.format(i,j)
#os.system(cmd)
ulx=-180+i*lon_inc
uly=90-j*lat_inc
lrx=-180+(i+1)*lon_inc
lry=90-(j+1)*lat_inc
'''
Georeference the image pieces.
'''
cmd = "gdal_translate -a_ullr {0} {1} {2} {3} ".format(ulx,uly,lrx,lry) + output_dir+'{0}_{1}.tif '.\
format(i,j) + output_dir+'{0}_{1}_ref.tif '.format(i,j)
print cmd
print os.system(cmd)
if not show_img_flag:
plt.close()
os.system("rm "+output_dir+'{0}_{1}.tif '.format(i,j))
print 'done'
In [36]:
def gmtColormap_1(filePath):
steps = []
colors = []
with open(filePath,'r') as f:
lines = f.readlines()
for l in lines:
ls = l.split()
if ls[0] == '#' or ls[0] == "B" or ls[0] == "F" or ls[0] == "N":
continue
else:
steps.append(float(ls[0]))
steps.append(float(ls[4]))
colors.append(hsv_to_rgb([float(ls[1])/360.0,float(ls[2]),float(ls[3])]))
colors.append(hsv_to_rgb([float(ls[5])/360.0,float(ls[6]),float(ls[7])]))
color_list = []
print len(steps), len(colors)
for i in range(len(steps)):
color_list.append((float(steps[i]-steps[0])/(steps[-1]-steps[0]), colors[i]))
#print color_list
return LinearSegmentedColormap.from_list('cmap', color_list, N=1024)
# test gmtColormap() function
cmap = gmtColormap_1(ROOT_PATH + "/topo.cpt")
dataset = gdal.Open(root_dir + '/GEBCO_2014_1D_pieces/5_5.tif', GA_ReadOnly )
band = dataset.GetRasterBand(1)
#r = band.ReadAsArray( 0, 0, band.XSize, band.YSize, band.XSize, band.YSize)
r = band.ReadAsArray( 0, 0, band.XSize, band.YSize, 3600, 1800)
fig = plt.figure(figsize=(12, 6),dpi=240, frameon=False)
ax = plt.Axes(fig, [0., 0., 1., 1.])
ax.set_axis_off()
fig.add_axes(ax)
plt.ioff()
'''
You must be careful with the Nan data in the grid.
'''
where_are_NaNs = isnan(r)
r[where_are_NaNs] = 0
rgb = shading.shade(r,cmap=cmap,intensity=shading.intensity(r),vmin=-7000, vmax=7000)
plt.imshow(rgb)
plt.show()
In [33]:
fig = plt.figure(figsize=(8, 3),dpi=300,frameon=False)
ax1 = fig.add_axes([0.05, 0.80, 0.9, 0.15])
cmap = gmtColormap_1(ROOT_PATH + "/topo.cpt")
norm = matplotlib.colors.Normalize(vmin=-7000, vmax=7000, clip=False)
cb1 = matplotlib.colorbar.ColorbarBase(ax1, cmap=cmap, norm=norm, orientation='horizontal',
extend='both',spacing='uniform',
ticks=[-6000,-4000, -2000, 0, 2000, 4000, 6000])
cb1.set_label('meters')
plt.show()
In [ ]:
In [ ]: