How do I contour projected data?

In order for Gsharp to contour your 2D data, it must be presented as a grid of data and not scattered irregularly. Whether you will need to interpolate or not depends on the shape of your input data and on the projection you are performing.

The lettered steps are detailed below:

A - Project grid definition

With the mollweide and mercator projections lines of latitude remain horizontal and hence projected grids remain as grids. There is no need to interpolate our grid, we just need to change the datasets that define the location of the grid - xgrid and ygrid.

  #Read in your 24 by 36 grid
  import_ascii("gridfile.dat",,,,,,"grid",,24,36);
  longs=range(-180//180,24); dumlats=repeatx(0,24);
  lats=range(-90//90,36); dumlongs=repeatx(0,36);

  #do_mollweide can be found in libproj.gsl
  do_mollweide(longs,dumlats,"ugrid","dumlats");
  do_mollweide(dumlongs,lats,"dumlongs","vgrid");

  create Viewport page_1.view;
  create Domain page_1.view.domain;
  create Graph page_1.view.domain.graph
   ( XuNcolorDataGrid = "grid",
     XuNgraphType = "2DContour",
     XuNxData = "ugrid",
     XuNyData = "vgrid"
   );

B - Convert grid to x,y,z

The projection routines take scatter data not grids. A grid can be converted to x,y,z with the following commands:

  #Read in your 24 by 36 grid
  import_ascii("gridfile.dat",,,,,,"grid",,24,36);
  longs=range(-180//180,24); dumlats=repeatx(0,24);
  lats=range(-90//90,36); dumlongs=repeatx(0,36);

  #Convert to x,y,z
  x = list(repeaty(longs,size(lats)));
  y = list(repeatx(lats,size(longs)));
  z = list(grid);

C - Project x,y into u,v

Project x,y as you would any other 1D data, such as a worldmap or grid. For example:

  do_orthographic(12,55,x,y,"u","v");

D - Interpolate u,v,z into grid

There are two methods for interpolating scattered data into a grid, bilinear() and bivariate() - you choose.

  ugrid = range(u,40);
  vgrid = range(v,40);
  gridp = bilinear(x,y,z,ugrid,vgrid);
  #gridp = bivariate(x,y,z,ugrid,vgrid);

N.B. It is possible that you collected your data so that if formed a grid when it had been projected. If so replace this step with commands to reshape z into a grid. Possibly something like this:

  ugrid = unique(u);  vgrid = unique(v);
  gridp = reshape(x,40,40,1);

E - Calculate the projection boundary

It is possible that the desired plot is not rectangular, where as gridp is. You need to calculated the extremes of the projection and then plot gridp within this region.

Methods for calculating this extreme vary depend on the projection and the data. For example:

  #Orthographic and Stereographic plots require a circular
  #boundary of radius 1 and centre (0,0).
  i = range(0//2*pi,360);
  rx = sin(i);  ry = cos(i);

All done

You should now have a grid, gridp, with grid definitions ugrid and vgrid and possibly a region define in rx and ry. These can be plotted like so:

  create Viewport page_1.view;
  create Domain page_1.view.domain;
  create Graph page_1.view.domain.graph
   ( XuNcolorDataGrid = "gridp",
     XuNgraphType = "2DContour",
     #XuNregionXData = "rx",
     #XuNregionYData = "ry"
     XuNxData = "ugrid",
     XuNyData = "vgrid"
   );

Now you can add your worldmap and/or grid.