Skip to content

Instantly share code, notes, and snippets.

@briochemc
Last active June 15, 2019 08:14
Show Gist options
  • Save briochemc/59a84810f70911aab6a23293aebd92c7 to your computer and use it in GitHub Desktop.
Save briochemc/59a84810f70911aab6a23293aebd92c7 to your computer and use it in GitHub Desktop.
using Makie
using AbstractPlotting
using FileIO, Colors
# Load World Ocean Atlas 2D fields for PO₄, Si(OH)₄, and NO₄
using WorldOceanAtlasTools
lat, lon, PO₄_2D = WorldOceanAtlasTools.WOA13_surface_map("p", 0, "1°")
lat, lon, NO₄_2D = WorldOceanAtlasTools.WOA13_surface_map("n", 0, "1°")
lat, lon, SiOH₄_2D = WorldOceanAtlasTools.WOA13_surface_map("i", 0, "1°")
# Add the image layer to each plot
# (TODO update with my heatmaps)
# Julia_colors_RGB
darker_purple = RGB(0.584, 0.345, 0.698)
darker_green = RGB(0.22 , 0.596, 0.149)
darker_red = RGB(0.796, 0.235, 0.2 )
# convert z to RGB array
function array_to_RGB(A, col)
minimum(A)
a = minimum(A[.~isnan.(A)])
b = maximum(A[.~isnan.(A)])
[v_to_RGB(v, a, b, col) for v in A]
end
function v_to_RGB(v, a, b, col)
if isnan(v)
return RGB(0.0,0.0,0.0)
else
α = (v - a) / (b - a)
return α * RGB(1.0,1.0,1.0) + (1-α) * col
end
end
map1 = array_to_RGB(PO₄_2D, darker_red)
map2 = array_to_RGB(NO₄_2D, darker_purple)
map3 = array_to_RGB(SiOH₄_2D, darker_green)
# Reverse them because they are the wrong way it seems?
map1 = reverse(map1, dims=1)
map2 = reverse(map2, dims=1)
map3 = reverse(map3, dims=1)
# Optional add equator white line
#map1[90:91,:] .= RGB(1.0,1.0,1.0)
#map2[90:91,:] .= RGB(1.0,1.0,1.0)
#map3[90:91,:] .= RGB(1.0,1.0,1.0)
# Create the 3 speric meshes placed in a (3D) triangle
m1 = GLNormalUVMesh(Sphere(Point3f0(0), 1f0), 100)
m2 = GLNormalUVMesh(Sphere(Point3f0(0), 1f0), 100)
m3 = GLNormalUVMesh(Sphere(Point3f0(0), 1f0), 100)
# add `, show_axis = false` later to the `mesh` calls to remove axes.
s = mesh(m1, color = map1, shading = false, show_axis = false)
sph1 = s[end]
sph2 = mesh!(m2, color = map2, shading = false, show_axis = false)[end]
sph3 = mesh!(m3, color = map3, shading = false, show_axis = false)[end]
# Rotate the earths individually using quaternion rotations
# First, rotate globes for viewer to face the oceans (IND, PAC, ATL)
zaxis = Vec3f0(0, 0, 1) # z-axis to rotate to correct ocean
zrot1 = qrotation(zaxis, 1.9π) # ATL
zrot2 = qrotation(zaxis, 1.2π) # IND
zrot3 = qrotation(zaxis, 0.40π) # PAC
# Second, rotate globes to "face" the camera's eye position
xyaxis = Vec3f0(-1/sqrt(2), 1/sqrt(2), 0)
xyrot =qrotation(xyaxis, -35.26/180*π)
# Rotate them all by conposing the rotations
rotate!(sph1, xyrot * zrot1)
rotate!(sph2, xyrot * zrot2)
rotate!(sph3, xyrot * zrot3)
# Translate the spheres to symmetric positions
d = 2.4 # distance of the globe centers from origin
translate!(sph1, Point3f0(d, 0, 0))
translate!(sph2, Point3f0(0, d, 0))
translate!(sph3, Point3f0(0, 0, d))
# TODO Make sure the camera is centered
# This is not working yet.
# I mean, it works in the plot I am interactively looking at
# but if I call the scene `s` again, it has not changed...
# And so when I save `s`, it does not have the updated cam
#s
#update_cam!(s, s.camera, Vec3f0(4), Vec3f0(0))
#s.camera.eyeposition[] = Vec3f0(4)
#rotate_cam!(s, Vec3f0(0,0,0.1)
# TODO Adding cycling arrows
# (Maybe leave to external tool like GIMP or Inkscape
save("AIBECS_logo.png", s, resolution=(1000,1000))
# Load World Ocean Atlas 2D fields for PO₄, Si(OH)₄, and NO₄
using WorldOceanAtlasTools
product_year = 2018
tracer = "p"
period = 0
resolution = "1°"
ds = WorldOceanAtlasTools.WOA_Dataset(product_year, tracer, period, resolution, verbose=false) ;
field = "an"
field3D, lat, lon, depth = WorldOceanAtlasTools.get_gridded_3D_field(ds, tracer, field) ;
map_2D = field3D[:,:,1] ;
# Julia_colors_RGB
using Colors
darker_purple = RGB(0.584, 0.345, 0.698)
darker_green = RGB(0.22 , 0.596, 0.149)
darker_red = RGB(0.796, 0.235, 0.2 )
# colrmap
C(g::ColorGradient, levels) = [(g[z].r, g[z].g, g[z].b) for z=range(0.0,1.0,length=length(levels))]
g = cgrad([:white, darker_red])
levels = 0:0.2:1.6
# Use Cartopy
ENV["MPLBACKEND"]="qt5agg"
"qt5agg"
using PyPlot, PyCall, Plots
ccrs = pyimport("cartopy.crs")
clf()
# Define the projection
ax = subplot(projection=ccrs.NearsidePerspective(central_latitude=10.0, central_longitude=170.0, satellite_height=100000000.0))
# Add coastline
# TODO fill it with black instead of coastlines?
ax.coastlines()
# Make the data cyclic
lon_cyc = [lon; 360+lon[1]] # making it cyclic for Cartopy
map_cyc = hcat(map_2D, map_2D[:,1])
map_cyc = map(x -> ismissing(x) ? NaN : min(x, levels[end]), map_cyc)
# contour
p = contourf(lon_cyc, lat, map_cyc, levels=levels, transform=ccrs.PlateCarree(), zorder=-1, colors=C(g,levels))
#ax.gridlines(xlocs=collect(-180:3:180), ylocs=collect(-90:3:90), linewidth = 2, color="black", alpha = 1.0)
Pyplot.savefig("AIBECS_logo_red.svg")
% A simple cycle
% Author : Jerome Tremblay
\documentclass[tikz,border=3mm]{standalone}
\usepackage{tikz}
\newcommand{\arcarrow}[8]{%
% inner radius, middle radius, outer radius, start angle,
% end angle, tip protusion angle, options, text
\pgfmathsetmacro{\rin}{#1}
\pgfmathsetmacro{\rmid}{#2}
\pgfmathsetmacro{\rout}{#3}
\pgfmathsetmacro{\astart}{#4}
\pgfmathsetmacro{\aend}{#5}
\pgfmathsetmacro{\atip}{#6}
\fill[#7] (\astart:\rin) arc (\astart:\aend:\rin)
-- (\aend+\atip:\rmid) -- (\aend:\rout) arc (\aend:\astart:\rout)
-- (\astart+\atip:\rmid) -- cycle;
\draw[#8] (\astart:\rin) arc (\astart:\aend:\rin)
-- (\aend+\atip:\rmid) -- (\aend:\rout) arc (\aend:\astart:\rout)
-- (\astart+\atip:\rmid) -- cycle;
}
\begin{document}
\begin{tikzpicture}
\def \ballsize {.4\textwidth}
\def \blanksize {.4\textwidth}
%\def \blanksize {.5\textwidth} % for straight arrows
\def \radius {5}
\def \radiusarrows {4}
\def \rin {3.7}
\def \rout {4.3}
\def \tip {-3}
\def \margin {35} % margin in angles, depends on the radius
\def \marginbis {80} % margin in angles, depends on the radius
\def \marginter {28} % margin in angles, depends on the radius
% Arc Arrows
%\fill[even odd rule,black!20!white] circle (\rin) circle (\rout);
%\arcarrow{\rin}{\radiusarrows}{\rout}{ 90}{ -30+\marginter+5}{\tip}{black!0!white}
%\arcarrow{\rin}{\radiusarrows}{\rout}{-30}{-150+\marginter+5}{\tip}{black!0!white}
%\arcarrow{\rin}{\radiusarrows}{\rout}{210}{ 90+\marginter+5}{\tip}{black!0!white}
\arcarrow{\rin}{\radiusarrows}{\rout}{ 90-\margin}{ -30+\margin-\tip}{\tip}{white}{black!100!white}
\arcarrow{\rin}{\radiusarrows}{\rout}{-30-\margin}{-150+\margin-\tip}{\tip}{white}{black!100!white}
\arcarrow{\rin}{\radiusarrows}{\rout}{210-\margin}{ 90+\margin-\tip}{\tip}{white}{black!100!white}
%\arcarrow{\rin}{\radiusarrows}{\rout}{ 90}{ -30+\marginbis+5}{\tip}{black!0!white}
%\arcarrow{\rin}{\radiusarrows}{\rout}{-30}{-150+\marginbis+5}{\tip}{black!0!white}
%\arcarrow{\rin}{\radiusarrows}{\rout}{210}{ 90+\marginbis+5}{\tip}{black!0!white}
% Globes
\node[circle, fill=white, inner sep=0pt, minimum size=\blanksize] (GB_blank) at ( {90}:\radius) {};
\node[circle, fill=white, inner sep=0pt, minimum size=\blanksize] (RB_blank) at ( {210}:\radius) {};
\node[circle, fill=white, inner sep=0pt, minimum size=\blanksize] (PB_blank) at ( {-30}:\radius) {};
\node (GB) at ( {90}:\radius) {\includegraphics[width=\ballsize]{green_ball.eps}} ;
\node (RB) at ({210}:\radius) {\includegraphics[width=\ballsize]{red_ball.eps}} ;
\node (PB) at ({330}:\radius) {\includegraphics[width=\ballsize]{purple_ball.eps}};
%% Straight arrows
%\draw[->, ultra thick, black] (GB_blank) -- (PB_blank) ;
%\draw[->, ultra thick, black] (PB_blank) -- (RB_blank) ;
%\draw[->, ultra thick, black] (RB_blank) -- (GB_blank) ;
\end{tikzpicture}
\end{document}
# Load World Ocean Atlas 2D fields for PO₄, Si(OH)₄, and NO₄
using WorldOceanAtlasTools
product_year = 2018
#tracer = "p"
tracer = "i"
#tracer = "n"
period = 0
resolution = "1°"
ds = WorldOceanAtlasTools.WOA_Dataset(product_year, tracer, period, resolution, verbose=false) ;
field = "an"
field3D, lat, lon, depth = WorldOceanAtlasTools.get_gridded_3D_field(ds, tracer, field) ;
map_2D = field3D[:,:,1] ;
# Use Cartopy
ENV["MPLBACKEND"]="qt5agg"
using PyPlot, PyCall
ccrs = pyimport("cartopy.crs")
clf()
# Julia_colors_RGB
using Colors, Plots
#c = RGB(0.584, 0.345, 0.698) # darker purple
c = RGB(0.22 , 0.596, 0.149) # darker green
#c = RGB(0.796, 0.235, 0.2 ) # darker red
white = RGB(1.0, 1.0, 1.0)
# colrmap
C(g::ColorGradient, levels) = [(g[z].r, g[z].g, g[z].b) for z=range(0.0,1.0,length=length(levels))]
g = cgrad([white, c])
#levels = 0:0.2:2
levels = range(-10.0, 50.0, length=10)
#levels = 0:0.2:2
# Define the projection
central_lat, central_lon = 43.7384, 7.4246 # Monaco
central_lat, central_lon = 33.6846, -117.8265 + 360 # Irvine
central_lat, central_lon = -33.8688, 151.2093 # Sydney
ax = subplot(projection=ccrs.NearsidePerspective(central_latitude=central_lat, central_longitude=central_lon, satellite_height=100000000.0))
# Make background transparent?
#ax.outline_patch.set_visible(false)
ax.background_patch.set_visible(false)
# Add black land
cfeature = pyimport("cartopy.feature")
#ax.add_feature(cfeature.COASTLINE, edgecolor="#000000") # black coast lines
ax.add_feature(cfeature.LAND, facecolor="#000000") # gray land
# Make the data cyclic
lon_cyc = [lon; 360+lon[1]] # making it cyclic for Cartopy
map_cyc = hcat(map_2D, map_2D[:,1])
map_cyc = map(x -> ismissing(x) ? NaN : min(x, levels[end]), map_cyc)
# contour
p = PyPlot.contourf(lon_cyc, lat, map_cyc, levels=levels, transform=ccrs.PlateCarree(), zorder=-1, colors=C(g,levels))
#ax.gridlines(xlocs=collect(-180:3:180), ylocs=collect(-90:3:90), linewidth = 2, color="black", alpha = 1.0)
#
PyPlot.savefig("green_ball.svg")
PyPlot.savefig("green_ball.eps")
# Load World Ocean Atlas 2D fields for PO₄, Si(OH)₄, and NO₄
using WorldOceanAtlasTools
product_year = 2018
#tracer = "p"
#tracer = "i"
tracer = "n"
period = 0
resolution = "1°"
ds = WorldOceanAtlasTools.WOA_Dataset(product_year, tracer, period, resolution, verbose=false) ;
field = "an"
field3D, lat, lon, depth = WorldOceanAtlasTools.get_gridded_3D_field(ds, tracer, field) ;
map_2D = field3D[:,:,1] ;
# Use Cartopy
ENV["MPLBACKEND"]="qt5agg"
using PyPlot, PyCall
ccrs = pyimport("cartopy.crs")
clf()
# Julia_colors_RGB
using Colors, Plots
c = RGB(0.584, 0.345, 0.698) # darker purple
#c = RGB(0.22 , 0.596, 0.149) # darker green
#c = RGB(0.796, 0.235, 0.2 ) # darker red
white = RGB(1.0, 1.0, 1.0)
# colrmap
C(g::ColorGradient, levels) = [(g[z].r, g[z].g, g[z].b) for z=range(0.0,1.0,length=length(levels))]
g = cgrad([white, c])
#levels = range(0.0, 1.5, length=10)
#levels = range(0.0, 50.0, length=10)
levels = range(-1.5, 8.0, length=10)
# Define the projection
central_lat, central_lon = 43.7384, 7.4246 # Monaco
#central_lat, central_lon = 33.6846, -117.8265 + 360 # Irvine
#central_lat, central_lon = -33.8688, 151.2093 # Sydney
ax = subplot(projection=ccrs.NearsidePerspective(central_latitude=central_lat, central_longitude=central_lon, satellite_height=100000000.0))
# Make background transparent?
#ax.outline_patch.set_visible(false)
ax.background_patch.set_visible(false)
# Add black land
cfeature = pyimport("cartopy.feature")
#ax.add_feature(cfeature.COASTLINE, edgecolor="#000000") # black coast lines
ax.add_feature(cfeature.LAND, facecolor="#000000") # gray land
# Make the data cyclic
lon_cyc = [lon; 360+lon[1]] # making it cyclic for Cartopy
map_cyc = hcat(map_2D, map_2D[:,1])
map_cyc = map(x -> ismissing(x) ? NaN : min(x, levels[end]), map_cyc)
# contour
p = PyPlot.contourf(lon_cyc, lat, map_cyc, levels=levels, transform=ccrs.PlateCarree(), zorder=-1, colors=C(g,levels))
#ax.gridlines(xlocs=collect(-180:3:180), ylocs=collect(-90:3:90), linewidth = 2, color="black", alpha = 1.0)
#
PyPlot.savefig("purple_ball.svg")
PyPlot.savefig("purple_ball.eps")
PyPlot.savefig("purple_ball.pdf")
# Load World Ocean Atlas 2D fields for PO₄, Si(OH)₄, and NO₄
using WorldOceanAtlasTools
product_year = 2018
tracer = "p"
#tracer = "i"
#tracer = "n"
period = 0
resolution = "1°"
ds = WorldOceanAtlasTools.WOA_Dataset(product_year, tracer, period, resolution, verbose=false) ;
field = "an"
field3D, lat, lon, depth = WorldOceanAtlasTools.get_gridded_3D_field(ds, tracer, field) ;
map_2D = field3D[:,:,1] ;
# Use Cartopy
ENV["MPLBACKEND"]="qt5agg"
using PyPlot, PyCall
ccrs = pyimport("cartopy.crs")
clf()
# Julia_colors_RGB
using Colors, Plots
#c = RGB(0.584, 0.345, 0.698) # darker purple
#c = RGB(0.22 , 0.596, 0.149) # darker green
c = RGB(0.796, 0.235, 0.2 ) # darker red
white = RGB(1.0, 1.0, 1.0)
# colrmap
C(g::ColorGradient, levels) = [(g[z].r, g[z].g, g[z].b) for z=range(0.0,1.0,length=length(levels))]
g = cgrad([white, c])
levels = range(-0.2, 1.5, length=10)
#levels = range(0.0, 50.0, length=10)
#levels = 0:0.2:2
# Define the projection
#central_lat, central_lon = 43.7384, 7.4246 # Monaco
central_lat, central_lon = 33.6846, -117.8265 + 360 # Irvine
#central_lat, central_lon = -33.8688, 151.2093 # Sydney
ax = subplot(projection=ccrs.NearsidePerspective(central_latitude=central_lat, central_longitude=central_lon, satellite_height=100000000.0))
# Make background transparent?
#ax.outline_patch.set_visible(false)
ax.background_patch.set_visible(false)
# Add black land
cfeature = pyimport("cartopy.feature")
#ax.add_feature(cfeature.COASTLINE, edgecolor="#000000") # black coast lines
ax.add_feature(cfeature.LAND, facecolor="#000000") # gray land
# Make the data cyclic
lon_cyc = [lon; 360+lon[1]] # making it cyclic for Cartopy
map_cyc = hcat(map_2D, map_2D[:,1])
map_cyc = map(x -> ismissing(x) ? NaN : min(x, levels[end]), map_cyc)
# contour
p = PyPlot.contourf(lon_cyc, lat, map_cyc, levels=levels, transform=ccrs.PlateCarree(), zorder=-1, colors=C(g,levels))
#ax.gridlines(xlocs=collect(-180:3:180), ylocs=collect(-90:3:90), linewidth = 2, color="black", alpha = 1.0)
#
PyPlot.savefig("red_ball.svg")
PyPlot.savefig("red_ball.eps")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment