Commit 6c6fcaa3 authored by Florian Ziemen's avatar Florian Ziemen
Browse files

basic plots

parent 628cc3ee
#!/bin/bash
#SBATCH --job-name=clouds # Specify job name
#SBATCH --partition=compute2 # Specify partition name
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=36
#SBATCH --time=02:00:00 # Set a limit on the total run time
#SBATCH --mail-type=FAIL # Notify user by email in case of job failure
#SBATCH --account=bb1153 # Charge resources on this project account
#SBATCH --output=clouds.eo%j # File name for standard output
#SBATCH --error=clouds.eo%j # File name for standard error output
#SBATCH --constraint=256G
set -vx
module load paraview/5.7-insitu
which pvbatch
md5sum $(which pvbatch)
scriptdir=$(pwd)n
if [ "$1" == MOLL ] ; then moll="MOLL world" ; fi
if [ "$1" == SPHERE ] ; then moll=SPHERE ; fi
for x in $@ ; do
if [ "$x" == MOLL -o "$x" == SPHERE ] ; then
continue
fi
pr=${x/ts/pr}
pvbatch clouds.py $moll $x $pr
done
# srun -l --propagate=STACK,CORE --cpu_bind=cores \
# --distribution=block:cyclic
# state file generated using paraview version 5.7.0-RC3-442-g1b4ad14
# ----------------------------------------------------------------
# setup views used in the visualization
# ----------------------------------------------------------------
# trace generated using paraview version 5.7.0-RC3-442-g1b4ad14
#
# To ensure correct image size when batch processing, please search
# for and uncomment the line `# renderView*.ViewSize = [*,*]`
#### import the simple module from the paraview
import glob
import subprocess as sp
import sys
import os
from math import pi, sin, cos, tan, sqrt
def sphere_proj (lat, lon) :
r = 10000
theta = (90 - lat) * pi / 180.
phi = lon * pi / 180.
z = r * cos (theta)
x = r * sin(theta) * cos (phi)
y = r * sin(theta) * sin (phi)
return (x, y, z)
def mollweide_proj (lat, lon):
lat = lat * pi / 180.
lon = lon * pi / 180.
lon_center = 0 * pi / 180.
scale_by = 120. # Used for the size of the final map. This matches Niklas' Projection of the ICON data.
theta = lat
if abs(lat) != (pi/2):
for i in range(1):# only do one call for Niklas' projection
theta = get_theta (lon, lat, theta)
x = scale_by * 2 * sqrt(2) / pi * (lon - lon_center) * cos (theta)
y = scale_by * sqrt(2) * sin (theta)
return (x, y)
def get_theta(lon, lat, theta):
theta_new = theta - ( 2 * theta + sin (2* theta) - pi * sin ( lat)) / ( 2+ 2* cos (2*theta))
return (theta_new)
files = sys.argv[1:]
mollweide = False
if files[0] == "MOLL" :
mollweide = True
files = files[1:]
sphere = False
if files[0] == "SPHERE" :
sphere = True
files = files[1:]
regions = {"barbados": {"latlon" : (14,-56),
"scale": 30,
"ts_range":(20, 45),
"timestamp_location" : [0.025, 0.94],
"colorbar_font_color" : (1, 1, 1)},
"europe": {"latlon" : (50,0),
"scale": 70,
"ts_range":(-25, 20),
"timestamp_location" : [0.025, 0.94],
"colorbar_font_color" : (1, 1, 1)},
"madagascar": {"latlon" : (-20, 60),
"scale": 70,
"ts_range":(0, 45),
"timestamp_location" : [0.025, 0.94],
"colorbar_font_color" : (1, 1, 1)},
"SE_asia": {"latlon" : (11, 130),
"scale": 70,
"ts_range":(-25, 35),
"timestamp_location" : [0.025, 0.94],
"colorbar_font_color" : (0, 0, 0)},
"world" : {"latlon" : (0,10),
"scale": 210,
"ts_range":(-50, 45),
"timestamp_location" : [0.025, 0.94],
"colorbar_font_color" : (1, 1, 1)}}
if (mollweide):
for key in regions.keys():
lat, lon = regions[key]['latlon']
x , y = mollweide_proj (lat ,lon )
regions[key]["position"] = (x,y,1000)
regions[key]["focus"] = (x,y,0)
if (sphere):
for key in regions.keys():
lat, lon = regions[key]['latlon']
x , y , z = sphere_proj (lat ,lon )
regions[key]["position"] = (x,y,z)
regions[key]["focus"] = (0,0,0)
def sort_location(obj, loc):
"""Make it a bit easier to assign a location to the legend"""
if type(loc) is type ("abc"):
print("String location "+ loc)
obj.WindowLocation = loc
else:
print("Array location ", loc)
obj.WindowLocation = "AnyLocation"
obj.Position = loc
region = False
if regions.get(files[0], region):
r = regions.get(files[0], region)
regions = { files[0] : r}
files = files[1:]
else:
if sphere:
regions.pop("world", True)
print (regions)
if "nwp" in files[1]:
scale = 4
else:
scale = 3600
print ("scaling precipitation by ", scale)
import numpy as np
from paraview.simple import *
# LoadPlugin('/sw/rhel6-x64/paraview/paraview-5.7.0-gcc71-osmesa/lib64/paraview-5.7/plugins/CDIReader/CDIReader.so', remote=False, ns=globals())
# LoadPlugin('/sw/rhel6-x64/paraview/paraview-5.7.0-gcc71-osmesa/lib64/paraview-5.7/plugins/NetCDFTimeAnnotationPlugin/NetCDFTimeAnnotationPlugin.so', remote=False, ns=globals())
# LoadPlugin('/sw/rhel6-x64/paraview/paraview-5.7.0-gcc71-osmesa/lib64/paraview-5.7/plugins/EmbossingRepresentations/EmbossingRepresentations.so', remote=False, ns=globals())
# LoadPlugin('/sw/rhel6-x64/paraview/paraview-5.7.0-gcc71-osmesa/lib64/paraview-5.7/plugins/StreamLinesRepresentation/StreamLinesRepresentation.so', remote=False, ns=globals())
# LoadPlugin('/sw/rhel6-x64/paraview/paraview-5.7.0-gcc71-osmesa/lib64/paraview-5.7/plugins/SurfaceLIC/SurfaceLIC.so', remote=False, ns=globals())
#### disable automatic camera reset on 'Show'
paraview.simple._DisableFirstRenderCameraReset()
# get the material library
materialLibrary1 = GetMaterialLibrary()
# Create a new 'Render View'
renderView1 = CreateView('RenderView')
renderView1.ViewSize = [1134, 906]
renderView1.InteractionMode = '2D'
renderView1.AxesGrid = 'GridAxes3DActor'
renderView1.OrientationAxesVisibility = 0
renderView1.CenterOfRotation = [-0.0300445556640625, 0.0, 0.0]
renderView1.StereoType = 'Crystal Eyes'
#renderView1.CameraPosition = [7.991613089729558, 112.69542620767417, 1466.0518274342055]
#renderView1.CameraFocalPoint = [7.991613089729558, 112.69542620767417, 0.0]
renderView1.CameraFocalDisk = 1.0
renderView1.Background = [0.32, 0.34, 0.43]
renderView1.BackEnd = 'OSPRay raycaster'
renderView1.OSPRayMaterialLibrary = materialLibrary1
SetActiveView(None)
# ----------------------------------------------------------------
# setup view layouts
# ----------------------------------------------------------------
# create new layout object 'Layout #1'
layout1 = CreateLayout(name='Layout #1')
layout1.AssignView(0, renderView1)
# ----------------------------------------------------------------
# restore active view
SetActiveView(renderView1)
# ----------------------------------------------------------------
# ----------------------------------------------------------------
# setup the data processing pipelines
# ----------------------------------------------------------------
# create a new 'CDIReader'
print(files[0])
icon_reader = CDIReader(FileNames=[files[0]])
icon_reader.Dimensions = '(time, ncells)'
icon_reader.CellArrayStatus = ['ts']
if (mollweide):
icon_reader.SetProjection = 'Mollweide Projection'
elif sphere:
icon_reader.SetProjection = 'Spherical Projection'
icon_reader.LayerThickness = 50
print(icon_reader.TimestepValues)
print(files[1])
icon_reader2 = CDIReader(FileNames=[files[1]])
icon_reader2.Dimensions = '(time, ncells)'
icon_reader2.CellArrayStatus = ['cllvi', 'clivi']
if (mollweide):
icon_reader2.SetProjection = 'Mollweide Projection'
elif sphere :
icon_reader2.SetProjection = 'Spherical Projection'
icon_reader2.LayerThickness = 50
calculator1 = Calculator(Input=icon_reader2)
calculator1.AttributeType = 'Cell Data'
calculator1.ResultArrayName = 'pr'
calculator1.Function = 'cllvi+5*clivi'
# create a new 'NetCDF Time Annotation'
netCDFTimeAnnotation1 = NetCDFTimeAnnotation(Input=icon_reader)
netCDFTimeAnnotation1.Expression = '"%02i-%02i-%02i %02i:%02i UTC" % (Date[0], Date[1], Date[2], Date[3], Date[4])' #
# ----------------------------------------------------------------
# setup the visualization in view 'renderView1'
# ----------------------------------------------------------------
# show data from icon_reader
icon_readerDisplay = Show(icon_reader, renderView1)
# get color transfer function/color map for 'ts'
tsLUT = GetColorTransferFunction('ts')
tsLUT.AutomaticRescaleRangeMode = 'Never'
tsLUT.ApplyPreset('Inferno (matplotlib)', True)
tsLUT.NanColor = [0.0, 1.0, 0.0]
tsLUT.ScalarRangeInitialized = 1.0
if not (sphere or mollweide) :
tsr = regions["world"]["ts_range"]
tsLUT.RescaleTransferFunction(tsr[0],tsr[1])
# get opacity transfer function/opacity map for 'ts'
tsPWF = GetOpacityTransferFunction('ts')
tsPWF.Points = [250.0, 0.0, 0.5, 0.0, 320.0, 1.0, 0.5, 0.0]
tsPWF.ScalarRangeInitialized = 1
# trace defaults for the display properties.
icon_readerDisplay.Representation = 'Surface'
icon_readerDisplay.ColorArrayName = ['CELLS', 'ts']
icon_readerDisplay.LookupTable = tsLUT
icon_readerDisplay.OSPRayScaleFunction = 'PiecewiseFunction'
icon_readerDisplay.SelectOrientationVectors = 'None'
icon_readerDisplay.ScaleFactor = 59.99998779296875
icon_readerDisplay.SelectScaleArray = 'None'
icon_readerDisplay.GlyphType = 'Arrow'
icon_readerDisplay.GlyphTableIndexArray = 'None'
icon_readerDisplay.GaussianRadius = 2.9999993896484374
icon_readerDisplay.SetScaleArray = [None, '']
icon_readerDisplay.ScaleTransferFunction = 'PiecewiseFunction'
icon_readerDisplay.OpacityArray = [None, '']
icon_readerDisplay.OpacityTransferFunction = 'PiecewiseFunction'
icon_readerDisplay.DataAxesGrid = 'GridAxesRepresentation'
icon_readerDisplay.PolarAxes = 'PolarAxesRepresentation'
icon_readerDisplay.ScalarOpacityFunction = tsPWF
icon_readerDisplay.ScalarOpacityUnitDistance = 2.432556203035958
icon_readerDisplay.BumpMapInputDataArray = [None, '']
icon_readerDisplay.ExtrusionInputDataArray = ['CELLS', 'clivi']
icon_readerDisplay.InputVectors = [None, '']
icon_readerDisplay.SelectInputVectors = [None, '']
icon_readerDisplay.WriteLog = ''
# show data from netCDFTimeAnnotation1
netCDFTimeAnnotation1Display = Show(netCDFTimeAnnotation1, renderView1)
netCDFTimeAnnotation1Display.Color = [1., 1., 1.]
netCDFTimeAnnotation1Display.Bold = 1
netCDFTimeAnnotation1Display.Modified()
# show data from icon_reader2
icon_reader2Display = Show(calculator1, renderView1)
# get color transfer function/color map for 'pr'
prLUT = GetColorTransferFunction('pr')
prLUT.AutomaticRescaleRangeMode = 'Never'
prLUT.EnableOpacityMapping = 1
prLUT.RGBPoints = [0.0, 1.0, 1.0, 1.0, 2., 1.0, 1.0, 1.0] # value, r, g, b, value , r, g, b
prLUT.ColorSpace = 'Lab'
prLUT.ScalarRangeInitialized = 1.0
# get opacity transfer function/opacity map for 'pr'
prPWF = GetOpacityTransferFunction('pr')
prPWF.Points = [0.0, 0.0, 0.5, 0.0, 2., 1.0, 0.5, 0.0] # same as above with r = opacity, g = 0.5, b=0 (whyever)
prPWF.ScalarRangeInitialized = 1
# trace defaults for the display properties.
icon_reader2Display.Representation = 'Surface'
icon_reader2Display.ColorArrayName = ['CELLS', 'pr']
icon_reader2Display.LookupTable = prLUT
icon_reader2Display.OSPRayScaleFunction = 'PiecewiseFunction'
icon_reader2Display.SelectOrientationVectors = 'None'
icon_reader2Display.ScaleFactor = 67.87527770996094
icon_reader2Display.SelectScaleArray = 'None'
icon_reader2Display.GlyphType = 'Arrow'
icon_reader2Display.GlyphTableIndexArray = 'None'
icon_reader2Display.GaussianRadius = 3.3937638854980468
icon_reader2Display.SetScaleArray = [None, '']
icon_reader2Display.ScaleTransferFunction = 'PiecewiseFunction'
icon_reader2Display.OpacityArray = [None, '']
icon_reader2Display.OpacityTransferFunction = 'PiecewiseFunction'
icon_reader2Display.DataAxesGrid = 'GridAxesRepresentation'
icon_reader2Display.PolarAxes = 'PolarAxesRepresentation'
icon_reader2Display.ScalarOpacityFunction = prPWF
icon_reader2Display.ScalarOpacityUnitDistance = 2.751897456869179
icon_reader2Display.BumpMapInputDataArray = [None, '']
icon_reader2Display.ExtrusionInputDataArray = ['CELLS', 'clivi']
icon_reader2Display.InputVectors = [None, '']
icon_reader2Display.SelectInputVectors = [None, '']
icon_reader2Display.WriteLog = ''
# setup the color legend parameters for each legend in this view
# get color legend/bar for tsLUT in view renderView1
tsLUTColorBar = GetScalarBar(tsLUT, renderView1)
tsLUTColorBar.WindowLocation = 'UpperRightCorner'
tsLUTColorBar.Title = 'Surface temperature in deg C '
tsLUTColorBar.ComponentTitle = ''
tsLUTColorBar.RangeLabelFormat = '%.0f'
# set color bar visibility
tsLUTColorBar.Visibility = 1
# get color legend/bar for prLUT in view renderView1
prLUTColorBar = GetScalarBar(prLUT, renderView1)
prLUTColorBar.Title = u'clouds (liquid + 5*ice) in kg / m²'
prLUTColorBar.ComponentTitle = ''
prLUTColorBar.RangeLabelFormat = '%.0f'
# set color bar visibility
prLUTColorBar.Visibility = 1
# show color legend
icon_readerDisplay.SetScalarBarVisibility(renderView1, True)
# show color legend
icon_reader2Display.SetScalarBarVisibility(renderView1, True)
# ----------------------------------------------------------------
# setup color maps and opacity mapes used in the visualization
# note: the Get..() functions create a new object, if needed
# ----------------------------------------------------------------
# ----------------------------------------------------------------
# finally, restore active source
SetActiveSource(icon_reader)
# ----------------------------------------------------------------
reader=GetActiveSource()
tsteps=reader.TimestepValues
print ("Timesteps: ", ", ".join((str (x) for x in tsteps)))
view = GetActiveView()
#view.UseGradientBackground = 1
#view.Background2 = [0.4, 0.4, 0.4] # top
view.Background = [0.2, 0.2, 0.2] # bottom
var="ts-pr_"
basepath = '/work/bk1040/k202134/DYAMOND_WINTER/Images'
renderView1.UseLight = 0 # NO 3D, so we don't need fancy light, and can have full colors
my_id = os.path.basename(os.getcwd())+"_"+var
if len(files) == 1 :
my_id += "_"+os.path.basename(files[0][:-3])
else:
my_id += "_"+os.path.basename(files[0][:-3])+"--"+os.path.basename(files[-1][:-3])
for n,t in enumerate (tsteps):
print ("rendering for time %f"%t)
view.ViewTime = t
prLUTColorBar.LabelFontSize = 16
tsLUTColorBar.LabelFontSize = 16
prLUTColorBar.TitleFontSize = 16
tsLUTColorBar.TitleFontSize = 16
if not (mollweide or sphere):
renderView1.CameraPosition = [0., 0.0, 1500.0] # die Höhe ist für ’ne Paralllel-Projektion eigentlich egal, es sei denn, es sind Schichten / Objekte im Weg.
renderView1.CameraFocalPoint = [0.0, 0.0, 0.0]
renderView1.CameraParallelScale = 149.9 # dann sind die Ränder leicht beschnitten, und die fehlenden Zellen tauchen nicht auf.
netCDFTimeAnnotation1Display.WindowLocation = 'AnyLocation' # top-left
netCDFTimeAnnotation1Display.Position = [0.01, 0.46]
netCDFTimeAnnotation1Display.FontSize = 4
icon_readerDisplay.SetScalarBarVisibility(renderView1, False)
icon_reader2Display.SetScalarBarVisibility(renderView1, False)
SaveScreenshot(basepath + '/%s_4p_%06d.png'%(my_id, n), ImageResolution=(4096,2048))
icon_readerDisplay.SetScalarBarVisibility(renderView1, True)
icon_reader2Display.SetScalarBarVisibility(renderView1, True)
tsLUTColorBar.WindowLocation = 'AnyLocation' # center left for crystal sphere
prLUTColorBar.WindowLocation = 'AnyLocation' # center left for crystal sphere
SaveScreenshot(basepath + '/%s_4p_cb_%06d.png'%(my_id, n), ImageResolution=(4096,2048))
else:
netCDFTimeAnnotation1Display.FontSize = 4
sort_location(tsLUTColorBar, [0.92, 0.05])
sort_location(prLUTColorBar, [0.92, 0.63])
prLUTColorBar.LabelFontSize = 8
tsLUTColorBar.LabelFontSize = 8
prLUTColorBar.TitleFontSize = 8
tsLUTColorBar.TitleFontSize = 8
prLUTColorBar.TitleBold = 1
prLUTColorBar.LabelBold = 1
tsLUTColorBar.TitleBold = 1
tsLUTColorBar.LabelBold = 1
for region in sorted(regions.keys()):
conf = regions[region]
prLUTColorBar.TitleColor = conf["colorbar_font_color"]
prLUTColorBar.LabelColor = conf["colorbar_font_color"]
tsLUTColorBar.TitleColor = conf["colorbar_font_color"]
tsLUTColorBar.LabelColor = conf["colorbar_font_color"]
sort_location ( netCDFTimeAnnotation1Display, conf["timestamp_location"])
renderView1.CameraPosition = conf["position"]
renderView1.CameraFocalPoint = conf["focus"]
renderView1.CameraParallelScale = conf["scale"]
if sphere:
renderView1.CameraViewUp =(0,0,10000)
tsr = conf["ts_range"]
tsLUT.RescaleTransferFunction(tsr[0],tsr[1])
SaveScreenshot(basepath + '/%s_%s_4k_%06d.png'%(my_id, region, n), ImageResolution=(3840,2160))
SaveScreenshot(basepath + '/%s_%s_hd_%06d.png'%(my_id, region, n), ImageResolution=(1920,1080))
#!/bin/bash
#SBATCH --job-name=ts-pr # Specify job name
#SBATCH --partition=compute2 # Specify partition name
#SBATCH --partition=gpu # Specify partition name
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=36
#SBATCH --time=02:00:00 # Set a limit on the total run time
......@@ -8,7 +8,7 @@
#SBATCH --account=bb1153 # Charge resources on this project account
#SBATCH --output=pvjob.eo%j # File name for standard output
#SBATCH --error=pvjob.eo%j # File name for standard error output
#SBATCH --constraint=256G
#SBATCH --constraint=1024G
set -vx
......@@ -18,9 +18,11 @@ which pvbatch
md5sum $(which pvbatch)
scriptdir=$(pwd)n
if [ $1 == MOLL ] ; then moll=MOLL ; fi
if [ "$1" == MOLL ] ; then moll="MOLL world" ; fi
if [ "$1" == SPHERE ] ; then moll=SPHERE ; fi
for x in $@ ; do
if [ $x == MOLL ] ; then
if [ "$x" == MOLL -o "$x" == SPHERE ] ; then
continue
fi
pr=${x/ts/pr}
......
......@@ -14,37 +14,92 @@ import glob
import subprocess as sp
import sys
import os
from math import pi, sin, cos, tan, sqrt
def sphere_proj (lat, lon) :
r = 10000
theta = (90 - lat) * pi / 180.
phi = lon * pi / 180.
z = r * cos (theta)
x = r * sin(theta) * cos (phi)
y = r * sin(theta) * sin (phi)
return (x, y, z)
def mollweide_proj (lat, lon):
lat = lat * pi / 180.
lon = lon * pi / 180.
lon_center = 0 * pi / 180.
scale_by = 120. # Used for the size of the final map. This matches Niklas' Projection of the ICON data.
theta = lat
if abs(lat) != (pi/2):
for i in range(1):# only do one call for Niklas' projection
theta = get_theta (lon, lat, theta)
x = scale_by * 2 * sqrt(2) / pi * (lon - lon_center) * cos (theta)
y = scale_by * sqrt(2) * sin (theta)
return (x, y)
def get_theta(lon, lat, theta):
theta_new = theta - ( 2 * theta + sin (2* theta) - pi * sin ( lat)) / ( 2+ 2* cos (2*theta))
return (theta_new)
files = sys.argv[1:]
mollweide = False
if files[0] == "MOLL" :
mollweide = True
files = files[1:]
regions = {"barbados": {"position" : (-94,31,1000),
"focus" : (-94,31,0),
"scale": 16,
sphere = False
if files[0] == "SPHERE" :
sphere = True
files = files[1:]
regions = {"barbados": {"latlon" : (14,-56),
"scale": 30,
"ts_range":(20, 45),
"timestamp_location" : "UpperLeftCorner"},
"europe": {"position" : (8,113,1000),
"focus" : (8,113,0),
"scale": 32,
"ts_range":(-20, 20),
"timestamp_location" : [0.01, 0.86]},
"madagascar": {"position" : (80,-45,1000),
"focus" : (80,-45,0),
"scale": 64,
"ts_range":(0, 45),
"timestamp_location" : "UpperLeftCorner"},
"SE_asia": {"position" : (240,40,1000),
"focus" : (240,40,0),
"scale": 64,
"ts_range":(-20, 35),
"timestamp_location" : "LowerLeftCorner"},
"world" : {"position" : (0,0,1000),
"focus" : (0,0,0),
"scale": 195,
"timestamp_location" : [0.025, 0.94],
"colorbar_font_color" : (1, 1, 1)},
"europe": {"latlon" : (50,0),
"scale": 70,
"ts_range":(-25, 20),
"timestamp_location" : [0.025, 0.94],
"colorbar_font_color" : (1, 1, 1)},
"madagascar": {"latlon" : (-20, 60),
"scale": 70,
"ts_range":(0, 45),
"timestamp_location" : [0.025, 0.94],
"colorbar_font_color" : (1, 1, 1)},
"SE_asia": {"latlon" : (11, 130),
"scale": 70,
"ts_range":(-25, 35),
"timestamp_location" : [0.025, 0.94],
"colorbar_font_color" : (0, 0, 0)},
"world" : {"latlon" : (0,10),
"scale": 210,
"ts_range":(-50, 45),
"timestamp_location" : [0.01, 0.86]}}
"timestamp_location" : [0.025, 0.94],
"colorbar_font_color" : (1, 1, 1)}}
if (mollweide):
for key in regions.keys():
lat, lon = regions[key]['latlon']
x , y = mollweide_proj (lat ,lon )
regions[key]["position"] = (x,y,1000)
regions[key]["focus"] = (x,y,0)
if (sphere):
for key in regions.keys():
lat, lon = regions[key]['latlon']
x , y , z = sphere_proj (lat ,lon )
regions[key]["position"] = (x,y,z)
regions[key]["focus"] = (0,0,0)
def sort_location(obj, loc):
"""Make it a bit easier to assign a location to the legend"""
......@@ -59,10 +114,21 @@ def sort_location(obj, loc):
region = False
if regions.get(files[0], region):
region = files[0]
r = regions.get(files[0], region)
regions = { files[0] : r}