Climate Rendering

This page contains some tips for plotting climate data in VisIt.

Combining Ocean Data and Satellite Images

Sometimes you will want to combine plots of ocean data with real color imagery from satellites. For instance, you might show ocean temperature while using satellite images to color the continents. VisIt currently does not have an easy way to do this but there is a workaround. You can create 3D "image" files that use the Truecolor plot in VisIt. In this case, you could create the image of the data that you want to map onto the globe and then convert it to a VTK file that looks like a globe before plotting the data in VisIt using the Truecolor plot.

  1. Get the Pseudocolor of the ocean temperature. You can save this out using VisIt and then scale it to the same size as your satellite data image. Be sure your data uses the same map projection as the satellite and height images you're using so all of the continents line up.
  2. Bring the Pseudocolor image, satellite image, and a mask image of the continents into an image editing program like gimp and combine them as desired. Or, you could write a program to do this.
  3. Run the makeglobe program on the resulting JPEG image to create a globe that you can use in VisIt
  4. Run VisIt and open the globe VTK file and create a Truecolor plot
Combined climate globe.png
Combined satellite and temperature data

makeglobe

Makeglobe program is a simple VTK program that turns a JPEG image into a VTK file that can be used with VisIt's Truecolor plot. Makeglobe takes a JPEG image for the colors and a JPEG image called height.jpg for the height values that are used to displace mountain ranges. The resolutions of both images must match.

  • The makeglobe program currently does not take ANY lat/lon coordinates into account. It merely creates a sphere with the same number of cells as your JPEG image and maps the data from the JPEG image to the sphere.
Satellite rgb.jpg
Satellite RGB image
Satellite height.jpg
Height image


#include <vtkSphereSource.h>
#include <vtkPolyDataWriter.h>
#include <vtkPolyData.h>

#include <vtkUnsignedCharArray.h>
#include <vtkJPEGReader.h>
#include <vtkImageData.h>
#include <vtkCellData.h>
#include <vtkPointData.h>

#include <math.h>

//#define X_RES 675
//#define Y_RES 338

// Earth medium
#define X_RES 1350
#define Y_RES 675

// Mars image
//#define X_RES 460
//#define Y_RES 230

// Jupiter image
//#define X_RES 1799
//#define Y_RES 600


const float radius = 6378.f; // ~ Earth
const float pole_factor = 0.996;

//const float radius = 2110.f; // Mars

//const float radius = 71492.f; // Jupiter
//const float pole_factor = 0.9;

int
main(int argc, char *argv[])
{
    int npts = X_RES * (Y_RES-1) + 2;
    int ncells = X_RES * Y_RES;
    vtkPoints *points = vtkPoints::New();
    points->SetNumberOfPoints(npts);
    float *fptr = (float *)points->GetVoidPointer(0);

    vtkPolyData *pd = vtkPolyData::New();
    pd->SetPoints(points);
    pd->Allocate(ncells);

    for(int ip = 0; ip < Y_RES+1; ++ip)
    {
        if(ip == 0)
        {
            *fptr++ = 0.;
            *fptr++ = -radius * pole_factor;
            *fptr++ = 0.;
        }
        else if(ip == Y_RES)
        {
            *fptr++ = 0.;
            *fptr++ = radius * pole_factor;
            *fptr++ = 0.;
        }
        else
        {
            float poleAngle = float(ip) / float(Y_RES) * M_PI;
            float sign = (poleAngle > M_PI_2) ? -1. : 1;
            float y = radius * cos(poleAngle) * -1;
            float yrad = radius * sin(poleAngle);

            for(int ri = 0; ri < X_RES; ++ri)
            {
                float angle = float(ri) / float(X_RES) * -2. * M_PI;
                float x = yrad * cos(angle);
                float z = yrad * sin(angle);

                *fptr++ = x;
                *fptr++ = y * pole_factor;
                *fptr++ = z;
            }
        }
    }
    vtkIdType verts[4];

    int i;
    for(i = 0; i < X_RES; ++i)
    {
        if(i < X_RES-1)
        {
            verts[0] = 0;
            verts[1] = i+2;
            verts[2] = i+1;
        }
        else
        {
            verts[0] = 0;
            verts[1] = 1;
            verts[2] = i+1;
        }
        pd->InsertNextCell(VTK_TRIANGLE, 3, verts);
    }

    int row = 1;
    int nextrow = X_RES + 1;
    for(int j = 0; j < Y_RES - 2; ++j)
    {
        for(i = 0; i < X_RES; ++i)
        {
            if(i < X_RES-1)
            {
                verts[0] = row + i;
                verts[1] = row + i + 1;
                verts[2] = nextrow + i + 1;
                verts[3] = nextrow + i;
            }
            else
            {
                verts[0] = row + i;
                verts[1] = row;
                verts[2] = nextrow;
                verts[3] = nextrow + i;
            }
            pd->InsertNextCell(VTK_QUAD, 4, verts);
        }

        row += X_RES;
        nextrow += X_RES;
    }

    int last = X_RES * (Y_RES-1) + 1;
    for(i = 0; i < X_RES; ++i)
    {
        if(i < X_RES-1)
        {
            verts[0] = row + i;
            verts[1] = row + i + 1;
            verts[2] = last;
        }
        else
        {
            verts[0] = row + i;
            verts[1] = row;
            verts[2] = last;
        }
        pd->InsertNextCell(VTK_TRIANGLE, 3, verts);
    }

    vtkJPEGReader *imgr = vtkJPEGReader::New();
    imgr->SetFileName(argv[1]);
    imgr->Update();
    vtkImageData *img = imgr->GetOutput();

    if(img->GetPointData()->GetNumberOfArrays() > 0)
    {
        vtkDataArray *data = img->GetPointData()->GetArray(0);
        vtkUnsignedCharArray *ucdata = vtkUnsignedCharArray::New();
        ucdata->SetNumberOfComponents(4);
        ucdata->SetNumberOfTuples(data->GetNumberOfTuples());
        ucdata->SetName(data->GetName());
        unsigned char *ucptr = (unsigned char *)ucdata->GetVoidPointer(0);
        for(int i = 0; i < data->GetNumberOfTuples(); ++i)
        {
            float old[4];
            data->GetTuple(i, old);
            *ucptr++ = (unsigned char)old[0];
            *ucptr++ = (unsigned char)old[1];
            *ucptr++ = (unsigned char)old[2];
            *ucptr++ = 255;
        }

        pd->GetCellData()->AddArray(ucdata);
    }
#if 1
    // Deform the surface using a second image.
    vtkJPEGReader *imgr2 = vtkJPEGReader::New();
    imgr2->SetFileName("height.jpg");
    imgr2->Update();
    vtkImageData *img2 = imgr2->GetOutput();
    fptr = (float *)points->GetVoidPointer(0);
    fptr += 3;
    if(img2->GetPointData()->GetNumberOfArrays() > 0)
    {
        vtkDataArray *hdata = img2->GetPointData()->GetArray(0);
        for(i = 0; i < npts-2; ++i, fptr += 3)
        {
            float dx = fptr[0];
            float dy = fptr[1];
            float dz = fptr[2];
            float inv_r = 1. / sqrt(dx*dx + dy*dy + dz*dz);
            float vx = dx * inv_r;
            float vy = dy * inv_r;
            float vz = dz * inv_r;
            float h[4];
            hdata->GetTuple(i, h);
#define EXAGERRATED_SCALE 120.f
            float altitude = (h[0] / 255.) * EXAGERRATED_SCALE;
            fptr[0] += vx * altitude;
            fptr[1] += vy * altitude;
            fptr[2] += vz * altitude;
        }        
    }
#endif

    vtkPolyDataWriter *dsw = vtkPolyDataWriter::New();
    dsw->SetInput(pd);
    dsw->SetFileName(argv[2]);
    dsw->SetFileTypeToBinary();
    dsw->Write();

    return 0;
}

Makefile

Here is an old Makefile that was used to build the makeglobe program. If you try to use it, you will need to adapt it to your system and you may need to make adjustments to it to account for newer VTK library versions. In addition, this Makefile is geared towards building on Linux using g++.

CXX=g++
CXXFLAGS=-g -O2
VISITDIR=/usr/gapps/visit
VTKDIR=$(VISITDIR)/vtk/2003.10.28/linux_rhel3_gcc_3.2.3
VTKLIBDIR=$(VTKDIR)/lib
MESADIR=$(VISITDIR)/mesa/current/linux_rhel3_gcc_3.2.3
MESALIBDIR=$(MESADIR)/lib

CPPFLAGS=-I. -I$(VTKDIR) -I$(VTKDIR)/Common -I$(VTKDIR)/Filtering -I$(VTKDIR)/Graphics \
-I$(VTKDIR)/Hybrid -I$(VTKDIR)/IO -I$(VTKDIR)/Imaging -I$(VTKDIR)/Rendering \
-I/usr/local/include

LDFLAGS=-L$(VTKDIR)/lib -L$(MESADIR)/lib -Wl,-rpath,$(VTKDIR)/lib -Wl,-rpath,$(MESADIR)/lib

IOLIBS=-lvtkIO -lvtkpng -lvtktiff -lvtkjpeg -lvtkDICOMParser -lvtkzlib -lvtkexpat
LIBS=-lvtkCommon -lvtkGraphics -lvtkRendering -lvtkImaging -lvtkFiltering $(IOLIBS) \
-lvtkfreetype -lvtkftgl -lGL -lMesaGL -lOSMesa

PROGRAM=makeglobe

SRC=$(PROGRAM).C
OBJ=$(SRC:.C=.o)

all: $(PROGRAM)

$(PROGRAM): $(OBJ)
	$(CXX) -o $(PROGRAM) $(OBJ) $(CXXFLAGS) $(LDFLAGS) $(LIBS)

clean:
	$(RM) $(OBJ) $(PROGRAM)

.C.o:
	$(CXX) $(CXXFLAGS) $(CPPFLAGS) -c $<

Elevating satellite data by height data

Sometimes a globe is not the visualization of choice for climate data yet you may still want to elevate satellite data by height data. This can be achieved in VisIt using CMFE expressions and an Elevate operator.

  1. Run VisIt and open the Satellite_rgb.jpg image file
  2. Create the expression height = conn_cmfe(<Satellite_height.jpg:intensity>, ImageMesh) * 0.15
  3. Create a Truecolor plot of color
  4. Apply an Elevate operator and make it elevate by the new height expression
  5. Draw the plots
  • You can export this elevated geometry as a VTK file if you want to simplify your setup for next time. This can be especially helpful if you have applied coordinate transforms that let you register the image to another coordinate system.
Elevated map.png
Satellite image elevated by height

Transforming coordinates

It is possible to transform lat/lon coordinates to give different map projections using VisIt's expressions and the Displace operator. Here is an example where we use the Displace operator and expressions to transform to a coordinate system that places the north pole at the center of the visualization.

These expressions create a displacement vector that we can add to the original mesh's coordinates. The displacement is such that we'll deform the mesh from a regular grid (2D) into a circular shape around the north pole. The expressions assume that the mesh's extents are: lon[-180,180] and lat[-90,90].

Here are some vector expressions that displace the coordinate system. Be sure to make separate expressions when you enter them into VisIt's Expressions window:

Name Definition Purpose
c coords(mesh193x96) This expression captures the mesh's coordinates in 'c'
thetaradius {2. * 3.14159 * (c[0] + 180) / 360., 1. - ((c[1] + 90) / 180.)} This expression normalizes the mesh coordinates in c to [0,1] and then transforms them into a theta,radius pair with ranges [0,2*pi] and [0,1]
offset {thetaradius[1] * cos(thetaradius[0]), thetaradius[1] * sin(thetaradius[0])} - c Transform the thetaradius values into polar coordinates and subtract the original mesh coordinates
  1. Create a Pseudocolor plot of the data
  2. Add a Displace operator and use the offset expression as the displacement variable
  3. Draw
Latlon map.png Polar map.png
Original plot with lat/lon Transformed polar plot


Polar Plots

If using the code above to create 3D globes from a 2D image, you can also transform those to create polar plots. These plots will generate equal spacing in Latitude (causes distortion near the equator).

Here are some vector expressions that displace the coordinate system in 3D. Be sure to omit the comments when you enter them into VisIt's Expressions window:

Name Definition Purpose
phi -rad2deg(acos(coord(mesh)[1]/sqrt(coord(mesh)[0]^2+coord(mesh)[1]^2+coord(mesh)[2]^2)))+90 get angles (essentially Latitude and Longitude)
theta if(gt(coord(mesh)[2],0),180-rad2deg(acos(coord(mesh)[0]/sqrt(coord(mesh)[0]^2+coord(mesh)[2]^2))),-1*(180-rad2deg(acos(coord(mesh)[0]/sqrt(coord(mesh)[0]^2+coord(mesh)[2]^2))))) get angles (essentially Latitude and Longitude)
r_prime 1-abs(phi)/90 get new coordinates
x_prime -r_prime*cos(deg2rad(theta)) get new coordinates
z_prime r_prime*sin(deg2rad(theta)) get new coordinates
shift {x_prime-coord(mesh)[0],-(coord(mesh)[1]*1),z_prime-coord(mesh)[2]} determine how much to offset or shift the mesh
  1. Create a Pseudocolor plot of the data
  2. Add a Displace operator and use the shift expression as the displacement variable
  3. Draw
Polar2 map.jpg
Transformed polar plot