Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have a 3d raster dataset of atmospheric variables (lat x lon x altitude). I would like to sample this raster dataset along a 3d path.

Any idea how I can achieve this, preferable with qgis or postgis or rasdaman?

Edit: Actually, my data is 4d (time x lat x lon x altitude), but having a 3d solution would already be great.

share|improve this question
Could you please explain the sense in which a 3D "polygon" can be interpreted as a path? Do you mean that its boundary is intended as a path, or are you perhaps asking about sampling the raster throughout the polygon itself (including its interior)? And when you use "polygon," does that presuppose (as in the conventional mathematical sense) that it is planar? – whuber Dec 5 '12 at 21:07
@whuber of course you're right, this doesn't make sense. Forget about the polygon, I simply want to sample along a 3d path. – andreash Dec 5 '12 at 21:23
You probably need to bake your own trilinear interpolation routine with NumPy arrays. I have a bilinear interpolation solution for GDAL (2D) if you need inspiration. – Mike T Dec 5 '12 at 22:02
See also SciPy's multivariate interpolation functions, although I've experience some oddities. Worth checking though. – Mike T Dec 5 '12 at 22:46
thanks @MikeToews, I'm actually quite familiar with SciPy and interpolation ;) However, I didn't want to have to deal with projection issues etc. and thus would favour some GIS-based solution where I don't have to worry about such stuff. – andreash Dec 6 '12 at 6:54

Some linear interpolation stuff like this?

    # -----------------------------------------------------------
    # Pseudo code adopted from R
    # -----------------------------------------------------------
    # Linear interpolation of a 3D line segment in a fixed time interval
    # startT Start in Time
    # startX Start in X
    # startY Start in Y
    # startZ Start in Y
    # endT   End in Time
    # endX   End in X
    # endY   End in Y
    # endZ   End in Z
    # steps  if == 0 steps are calculater from by the time interval
    # -----------------------------------------------------------
    dataFrame function sampleLineTime4D(
       int startT, int endT, 
       double startX, double endX, 
       double startY, double endY, 
       double startZ, double endZ, 
       int steps=0) {

      # gradient
      double dt = endT - startT;
      double dx = endX - startX;
      double dy = endY - startY;
      double dz = endZ - startZ;
      if (steps == 0) steps = dt;

      # relative stepsize
      double st = dt / steps;
      double sx = dx / steps;
      double sy = dy / steps;
      double sz = dz / steps;

      # definition of an empty vector set for T X Y and Z
      Vector<double> vx = new Vector();
      Vector<double> vx = new Vector();
      Vector<double> vy = new Vector();
      Vector<double> vz = new Vector();

      # Iteration over te steps
      for (int s=0; s<steps; s++) {
      } # Ende der schleife 

      return (dataFrame(t=rt, x=rx, y=ry, z=rz));
share|improve this answer

You could certainly load this raster data into Rasdaman (using Postgresql as the underlying database), then use either the Petascope interface to provide a WCS/WCPS service and query the data that way, or use the internal RASQL query language to do the same querying. See the Jacobs University pages for examples of using Rasdaman for multi-dimensional data.

The details of how to do this are a little convoluted to detail here, but the Rasdaman development community is very active, and there is plenty of documentation.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.