diff --git a/src/e3sm_quickview/pipeline.py b/src/e3sm_quickview/pipeline.py index f7cb479..d475144 100644 --- a/src/e3sm_quickview/pipeline.py +++ b/src/e3sm_quickview/pipeline.py @@ -209,8 +209,11 @@ def Update(self, data_file, conn_file, force_reload=False): ) # Step 1: Extract and transform atmospheric data - atmos_extract = EAMTransformAndExtract( # noqa: F821 - registrationName="AtmosExtract", Input=self.data + atmos_center = EAMCenterMeridian( + registrationName="AtmosCenter", Input=self.data + ) + atmos_extract = EAMExtract( # noqa: F821 + registrationName="AtmosExtract", Input=atmos_center ) atmos_extract.LongitudeRange = [-180.0, 180.0] atmos_extract.LatitudeRange = [-90.0, 90.0] @@ -224,7 +227,11 @@ def Update(self, data_file, conn_file, force_reload=False): atmos_proj.Projection = self.projection atmos_proj.Translate = 0 atmos_proj.UpdatePipeline() - self.moveextents = atmos_proj.GetDataInformation().GetBounds() + + atmos_polydata = ExtractSurface( + registrationName="AtmosPolyData", Input=atmos_proj + ) + self.moveextents = atmos_polydata.GetDataInformation().GetBounds() # Step 3: Load and process continent outlines if self.globe is None: @@ -270,7 +277,7 @@ def Update(self, data_file, conn_file, force_reload=False): grid_proj.UpdatePipeline() # Step 8: Cache all projected views for rendering - self.views["atmosphere_data"] = OutputPort(atmos_proj, 0) + self.views["atmosphere_data"] = OutputPort(atmos_polydata, 0) self.views["continents"] = OutputPort(cont_proj, 0) self.views["grid_lines"] = OutputPort(grid_proj, 0) diff --git a/src/e3sm_quickview/plugins/eam_projection.py b/src/e3sm_quickview/plugins/eam_projection.py index 8e866b0..e60cc16 100644 --- a/src/e3sm_quickview/plugins/eam_projection.py +++ b/src/e3sm_quickview/plugins/eam_projection.py @@ -1,14 +1,20 @@ from paraview.simple import * from paraview.util.vtkAlgorithm import * from vtkmodules.numpy_interface import dataset_adapter as dsa -from vtkmodules.vtkCommonCore import vtkPoints +from vtkmodules.vtkCommonCore import ( + vtkPoints, + vtkIdList, +) from vtkmodules.vtkCommonDataModel import ( vtkPolyData, vtkCellArray, vtkPlane, ) from vtkmodules.vtkCommonTransforms import vtkTransform -from vtkmodules.vtkFiltersCore import vtkAppendFilter +from vtkmodules.vtkFiltersCore import ( + vtkAppendFilter, + vtkGenerateIds, +) from vtkmodules.vtkFiltersGeneral import ( vtkTransformFilter, vtkTableBasedClipDataSet, @@ -52,7 +58,6 @@ def ProcessPoint(point, radius): z = rho * math.cos(math.radians(phi)) return [x, y, z] - @smproxy.filter() @smproperty.input(name="Input") @smdomain.datatype( @@ -408,74 +413,189 @@ def RequestData(self, request, inInfo, outInfo): @smproxy.filter() @smproperty.input(name="Input") @smdomain.datatype( - dataTypes=["vtkPolyData", "vtkUnstructuredGrid"], composite_data_supported=False + dataTypes=["vtkPolyData"], composite_data_supported=False ) @smproperty.xml( """ - - + + + + """ ) -class EAMCenterMeridian(VTKPythonAlgorithmBase): +class EAMExtract(VTKPythonAlgorithmBase): def __init__(self): super().__init__( nInputPorts=1, nOutputPorts=1, outputType="vtkUnstructuredGrid" ) - self.project = 0 - self.cmeridian = 0.0 + self.longrange = [-180.0, 180.0] + self.latrange = [-90.0, 90.0] - def SetCentralMeridian(self, meridian): - if self.cmeridian != meridian: - self.cmeridian = meridian + def SetLongitudeRange(self, min, max): + if self.longrange[0] != min or self.longrange[1] != max: + self.longrange = [min, max] + self.Modified() + + def SetLatitudeRange(self, min, max): + if self.latrange[0] != min or self.latrange[1] != max: + self.latrange = [min, max] self.Modified() def RequestData(self, request, inInfo, outInfo): inData = self.GetInputData(inInfo, 0, 0) outData = self.GetOutputData(outInfo, 0) - - if self.cmeridian == 0: + if self.longrange == [-180.0, 180] and self.latrange == [-90, 90]: outData.ShallowCopy(inData) return 1 - split = (self.cmeridian - 180) if self.cmeridian > 0 else (self.cmeridian + 180) + box = vtkPVBox() + box.SetReferenceBounds( + self.longrange[0], + self.longrange[1], + self.latrange[0], + self.latrange[1], + -1.0, + 1.0, + ) + box.SetUseReferenceBounds(True) + extract = vtkPVClipDataSet() + extract.SetClipFunction(box) + extract.InsideOutOn() + extract.ExactBoxClipOn() + extract.SetInputData(inData) + extract.Update() - planeL = vtkPlane() - planeL.SetOrigin([split, 0.0, 0.0]) - planeL.SetNormal([-1, 0, 0]) - clipL = vtkTableBasedClipDataSet() - clipL.SetClipFunction(planeL) - clipL.SetInputData(inData) - clipL.Update() + outData.ShallowCopy(extract.GetOutput()) + return 1 - planeR = vtkPlane() - planeR.SetOrigin([split, 0.0, 0.0]) - planeR.SetNormal([1, 0, 0]) - clipR = vtkTableBasedClipDataSet() - clipR.SetClipFunction(planeR) - clipR.SetInputData(inData) - clipR.Update() - transFunc = vtkTransform() - transFunc.Translate(360, 0, 0) - transform = vtkTransformFilter() - transform.SetInputData(clipL.GetOutput()) - transform.SetTransform(transFunc) - transform.Update() - append = vtkAppendFilter() - append.AddInputData(clipR.GetOutput()) - append.AddInputData(transform.GetOutput()) - append.Update() +@smproxy.filter() +@smproperty.input(name="Input") +@smproperty.xml( + """ + + + + Sets the central meridian. + Commonly used central meridians (longitudes) (- represents West, + represents East, + 0 is Greenwitch prime meridian): + - 0 (Prime Meridian): Standard "Western" view. + - -90/-100: Centered on North America. + - 100/110: Centered on Asia. + - -150/-160: Centered on the Pacific Ocean. + - 20: Often used to center Europe and Africa. + + + """ +) +@smdomain.datatype( + dataTypes=["vtkPolyData", "vtkUnstructuredGrid"], composite_data_supported=False +) +class EAMCenterMeridian(VTKPythonAlgorithmBase): + '''Cuts an unstructured grid and re-arranges the pieces such that + the specified meridian is in the middle. Note that the mesh is + specified with bounds [0, 360], but the meridian is specified in the more + common bounds [-180, 180]. + ''' + def __init__(self): + super().__init__( + nInputPorts=1, nOutputPorts=1, outputType="vtkUnstructuredGrid" + ) + # common values: + self._center_meridian = 0 + self._cached_output = None + + def SetMeridian(self, center_meridian_): + ''' + Specifies the central meridian (longitude in the middle of the map) + ''' + self._center_meridian = center_meridian_ + self.Modified() - transFunc = vtkTransform() - transFunc.Translate(-(180 + split), 0, 0) - transform = vtkTransformFilter() - transform.SetInputData(append.GetOutput()) - transform.SetTransform(transFunc) - transform.Update() + def GetMeridian(self): + ''' + Returns the central meridian + ''' + return self._center_meridian - outData.ShallowCopy(transform.GetOutput()) + def RequestData(self, request, inInfo, outInfo): + inData = self.GetInputData(inInfo, 0, 0) + inPoints = inData.GetPoints() + inCellArray = inData.GetCells() + + outData = self.GetOutputData(outInfo, 0) + if self._cached_output and self._cached_output.GetMTime() > inPoints.GetMTime() and \ + self._cached_output.GetMTime() > inCellArray.GetMTime(): + # only scalars have been added or removed + cached_cell_data = self._cached_output.GetCellData() + + in_cell_data = inData.GetCellData() + + outData.ShallowCopy(self._cached_output) + out_cell_data = outData.GetCellData() + + out_cell_data.Initialize() + for i in range(in_cell_data.GetNumberOfArrays()): + in_array = in_cell_data.GetArray(i) + cached_array = cached_cell_data.GetArray(in_array.GetName()) + if cached_array and cached_array.GetMTime() >= in_array.GetMTime(): + # this scalar has been seen before + # simply add a reference in the outData + out_cell_data.AddArray(in_array) + else: + # this scalar is new + # we have to fill in the additional cells resulted from the clip + out_array = in_array.NewInstance() + array0 = cached_cell_data.GetArray(0) + out_array.SetNumberOfComponents(array0.GetNumberOfComponents()) + out_array.SetNumberOfTuples(array0.GetNumberOfTuples()) + out_array.SetName(in_array.GetName()) + out_cell_data.AddArray(out_array) + outData.cell_data[out_array.GetName()] = inData.cell_data[i][self._cached_output.cell_data['PedigreeIds']] + else: + generate_ids = vtkGenerateIds() + generate_ids.SetInputData(inData) + generate_ids.PointIdsOff() + generate_ids.SetCellIdsArrayName("PedigreeIds") + + cut_meridian = self._center_meridian + 180 + plane = vtkPlane() + plane.SetOrigin([cut_meridian, 0.0, 0.0]) + plane.SetNormal([-1, 0, 0]) + # vtkClipPolyData hangs + clipL = vtkTableBasedClipDataSet() + clipL.SetClipFunction(plane) + clipL.SetInputConnection(generate_ids.GetOutputPort()) + clipL.Update() + + plane.SetNormal([1, 0, 0]) + clipR = vtkTableBasedClipDataSet() + clipR.SetClipFunction(plane) + clipR.SetInputConnection(generate_ids.GetOutputPort()) + clipR.Update() + + transFunc = vtkTransform() + transFunc.Translate(-360, 0, 0) + transform = vtkTransformFilter() + transform.SetInputData(clipR.GetOutput()) + transform.SetTransform(transFunc) + transform.Update() + + append = vtkAppendFilter() + append.AddInputData(clipL.GetOutput()) + append.AddInputData(transform.GetOutput()) + append.Update() + outData.ShallowCopy(append.GetOutput()) + self._cached_output = outData.NewInstance() + self._cached_output.ShallowCopy(outData) return 1 diff --git a/src/e3sm_quickview/plugins/eam_reader.py b/src/e3sm_quickview/plugins/eam_reader.py index 37bf7b1..2c60339 100644 --- a/src/e3sm_quickview/plugins/eam_reader.py +++ b/src/e3sm_quickview/plugins/eam_reader.py @@ -221,6 +221,20 @@ def _markmodified(*args, **kwars): """ ) +@smproperty.xml( + """ + + + + If True, the points of the dataset will be float, otherwise they will be float or double depending + on the type of corner_lat and corner_lon variables in the connectivity file. + + + """ +) class EAMSliceSource(VTKPythonAlgorithmBase): def __init__(self): VTKPythonAlgorithmBase.__init__( @@ -230,6 +244,7 @@ def __init__(self): self._DataFileName = None self._ConnFileName = None + self._ForceFloatPoints = True self._dirty = False # Variables for dimension sliders @@ -448,7 +463,12 @@ def _build_geometry(self, meshdata): lat = meshdata[latdim][:].data.flatten() lon = meshdata[londim][:].data.flatten() - coords = np.empty((len(lat), 3), dtype=np.float64) + + if self._ForceFloatPoints: + points_type = np.float32 + else: + points_type = np.float64 if lat.dtype == np.float64 else np.float32 + coords = np.empty((len(lat), 3), dtype=points_type) coords[:, 0] = lon coords[:, 1] = lat coords[:, 2] = 0.0 @@ -579,6 +599,11 @@ def SetConnFileName(self, fname): self._populate_variable_metadata() self.Modified() + def SetForceFloatPoints(self, forceFloatPoints): + if self._ForceFloatPoints != forceFloatPoints: + self._ForceFloatPoints = forceFloatPoints + self.Modified() + def SetSlicing(self, slice_str): # Parse JSON string containing dimension slices and update self._slices