Commit 0b542bb8 authored by Florian Ziemen's avatar Florian Ziemen
Browse files

clouds.py on laptop conflicting with mistral

parent b0773191
# 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, 60),
"timestamp_location" : [0.025, 0.94],
"colorbar_font_color" : (1, 1, 1)},
"world_no-legend" : {"latlon" : (0,0),
"scale": 200,
"ts_range":(-50, 60),
"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}
region = files[0]
files = files[1:]
else:
if sphere:
for r in regions.keys():
if world in r:
regions.pop(r, True)
print (regions)
print (region)
if "nwp" in files[1]:
scale = 4
else:
scale = 3600
print ("scaling precipitation by ", scale)
#import numpy as np
from paraview.simple import *
LoadPlugin('/usr/local/lib/paraview-5.8/plugins/CDIReader/CDIReader.so', remote=False, ns=globals())
LoadPlugin('/usr/local/lib/paraview-5.8/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 'Legacy VTK Reader'
planar_coastlines_intermedres_riversvtk = LegacyVTKReader(FileNames=['/Users/flo/Paraview/Planar-Coastlines-0...360-for-ParaView/Planar_coastlines_intermedres_norivers.vtk'])
# show data in view
planar_coastlines_intermedres_riversvtkDisplay = Show(planar_coastlines_intermedres_riversvtk, renderView1, 'GeometryRepresentation')
# trace defaults for the display properties.
planar_coastlines_intermedres_riversvtkDisplay.Representation = 'Surface'
planar_coastlines_intermedres_riversvtkDisplay.ColorArrayName = [None, '']
planar_coastlines_intermedres_riversvtkDisplay.OSPRayScaleFunction = 'PiecewiseFunction'
planar_coastlines_intermedres_riversvtkDisplay.SelectOrientationVectors = 'None'
planar_coastlines_intermedres_riversvtkDisplay.ScaleFactor = 36.0
planar_coastlines_intermedres_riversvtkDisplay.SelectScaleArray = 'None'
planar_coastlines_intermedres_riversvtkDisplay.GlyphType = 'Arrow'
planar_coastlines_intermedres_riversvtkDisplay.GlyphTableIndexArray = 'None'
planar_coastlines_intermedres_riversvtkDisplay.GaussianRadius = 1.8
planar_coastlines_intermedres_riversvtkDisplay.SetScaleArray = [None, '']
planar_coastlines_intermedres_riversvtkDisplay.ScaleTransferFunction = 'PiecewiseFunction'
planar_coastlines_intermedres_riversvtkDisplay.OpacityArray = [None, '']
planar_coastlines_intermedres_riversvtkDisplay.OpacityTransferFunction = 'PiecewiseFunction'
planar_coastlines_intermedres_riversvtkDisplay.DataAxesGrid = 'GridAxesRepresentation'
planar_coastlines_intermedres_riversvtkDisplay.PolarAxes = 'PolarAxesRepresentation'
# update the view to ensure updated data information
renderView1.Update()
# create a new 'Programmable Filter'
programmableFilter1 = ProgrammableFilter(Input=planar_coastlines_intermedres_riversvtk)
programmableFilter1.Script = ''
programmableFilter1.RequestInformationScript = ''
programmableFilter1.RequestUpdateExtentScript = ''
programmableFilter1.PythonPath = ''
# Properties modified on programmableFilter1
programmableFilter1.Script = """
from math import pi, sin, cos, sqrt
pdi = self.GetPolyDataInput()
pdo = self.GetPolyDataOutput()
newPoints = vtk.vtkPoints()
numPoints = pdi.GetNumberOfPoints()
print("PolydataMollweide: found %d points to transform."%numPoints)
def real_mollweide (lon, lat):
lon = lon * pi / 180.
lat = lat * pi / 180.
lon_center = 0 * pi / 180.
scale_by = 120.
if lat == pi/2.:
theta = pi/2
elif lat == -pi/2.:
theta = -pi/2
else:
theta = get_theta (lon, lat, lat)
for i in range(5):
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)
for i in range(0, numPoints):
coord = pdi.GetPoint(i)
x, y, z = coord[:3]
x,y = real_mollweide(x,y)
z = z
newPoints.InsertPoint(i, x, y, z)
pdo.SetPoints(newPoints)"""
programmableFilter1.RequestInformationScript = ''
programmableFilter1.RequestUpdateExtentScript = ''
programmableFilter1.PythonPath = ''
# show data in view
programmableFilter1Display = Show(programmableFilter1, renderView1, 'GeometryRepresentation')
# trace defaults for the display properties.
programmableFilter1Display.Representation = 'Surface'
programmableFilter1Display.ColorArrayName = [None, '']
programmableFilter1Display.OSPRayScaleFunction = 'PiecewiseFunction'
programmableFilter1Display.SelectOrientationVectors = 'None'
programmableFilter1Display.ScaleFactor = 66.94908752441407
programmableFilter1Display.SelectScaleArray = 'None'
programmableFilter1Display.GlyphType = 'Arrow'
programmableFilter1Display.GlyphTableIndexArray = 'None'
programmableFilter1Display.GaussianRadius = 3.3474543762207034
programmableFilter1Display.SetScaleArray = [None, '']
programmableFilter1Display.ScaleTransferFunction = 'PiecewiseFunction'
programmableFilter1Display.OpacityArray = [None, '']
programmableFilter1Display.OpacityTransferFunction = 'PiecewiseFunction'
programmableFilter1Display.DataAxesGrid = 'GridAxesRepresentation'
programmableFilter1Display.PolarAxes = 'PolarAxesRepresentation'
# change solid color
programmableFilter1Display.AmbientColor = [0.0, 0.0, 0.0]
programmableFilter1Display.DiffuseColor = [0.0, 0.0, 0.0]
# Properties modified on programmableFilter1Display
programmableFilter1Display.Opacity = 0.4
# hide data in view
Hide(planar_coastlines_intermedres_riversvtk, renderView1)
# update the view to ensure updated data information
# 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
# 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 = '%.3f'
log_clouds=True
if (log_clouds):
# Rescale transfer function
prPWF.RescaleTransferFunction(0.001, 5.0)
prLUT.RescaleTransferFunction(0.001, 5.0)
# convert to log space
# prPWF.MapControlPointsToLogSpace()
# Properties modified on prPWF
# prPWF.UseLogScale = 1
# convert to log space
prLUT.MapControlPointsToLogSpace()
# Properties modified on prLUT
prLUT.UseLogScale = 1
# set color bar visibility
prLUTColorBar.Visibility = 1
# show color legend
icon_readerDisplay.SetScalarBarVisibility(renderView1, True)
icon_reader2Display.SetScalarBarVisibility(renderView1, True)
print (region)
def disable_colorbars_if_requested():
if region and "no-legend" in region:
print("de-activating colorbars")
icon_readerDisplay.SetScalarBarVisibility(renderView1, False)
icon_reader2Display.SetScalarBarVisibility(renderView1, False)
disable_colorbars_if_requested()
# ----------------------------------------------------------------
# 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'
basepath = '/tmp'
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))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment