[thredds] Thredds WMS support for large source grids

• To: Jon Blower <j.d.blower@xxxxxxxxxxxxx>
• Subject: [thredds] Thredds WMS support for large source grids
• From: Ben Caradoc-Davies <Ben.Caradoc-Davies@xxxxxxxx>
• Date: Mon, 5 Sep 2011 17:12:45 +0800

Jon,

You might remember writing (in PixelMap):

// TODO: watch out for overflows (would only happen with a very large grid!)
int sourceGridIndex = j * this.sourceGridISize + i;

http://mailman.unidata.ucar.edu/mailing_lists/archives/thredds/2011/msg00312.html

It also includes:

- access changes in HorizontalCoordSys to support unit testing

http://siss2.anu.edu.au/thredds/godiva2/godiva2.html?server=http://siss2.anu.edu.au/thredds/wms/ga/test/PRISM_UTM55_wgs84.nc

The test file has a single ubyte variable on a 461276x1505407 grid.

Kind regards,

--
Ben Caradoc-Davies <Ben.Caradoc-Davies@xxxxxxxx>
Software Engineering Team Leader
CSIRO Earth Science and Resource Engineering
Australian Resources Research Centre

diff --git a/src/META-INF/MANIFEST.MF b/src/META-INF/MANIFEST.MF
deleted file mode 100644
index d7c96d6..0000000
--- a/src/META-INF/MANIFEST.MF
+++ /dev/null
@@ -1,9 +0,0 @@
-Manifest-Version: 1.0
-Ant-Version: Apache Ant 1.7.1
-Created-By: 17.1-b03 (Sun Microsystems Inc.)
-Built-By: Unidata - edavis
-Built-On: 2011-02-04 20:58:25
-Implementation-Title: Source of ncWMS for TDS
-Implementation-Version: 1.0.tds.4.2.20110204.2058
-Implementation-Vendor: University of Reading
-
diff --git a/src/uk/ac/rdg/resc/ncwms/cdm/CdmUtils.java 
b/src/uk/ac/rdg/resc/ncwms/cdm/CdmUtils.java
index ba42d03..90dfa3a 100644
--- a/src/uk/ac/rdg/resc/ncwms/cdm/CdmUtils.java
+++ b/src/uk/ac/rdg/resc/ncwms/cdm/CdmUtils.java
@@ -84,14 +84,16 @@ public final class CdmUtils
 
     static
     {
+        Range zero_range = null;
         try
         {
-            ZERO_RANGE = new Range(0, 0);
+            zero_range = new Range(0, 0);
         }
         catch(InvalidRangeException ire)
         {
             // Won't happen so we do nothing
         }
+        ZERO_RANGE = zero_range;
     }
 
     /** Enforce non-instantiability */
@@ -218,7 +220,7 @@ public final class CdmUtils
     public static DataReadingStrategy 
getOptimumDataReadingStrategy(NetcdfDataset nc)
     {
         String fileType = nc.getFileTypeId();
-        return fileType.equals("netCDF") || fileType.equals("HDF4")
+        return fileType.equals("netCDF") || fileType.equals("HDF4") || 
fileType.equals("HDF5")
             ? DataReadingStrategy.SCANLINE
             : DataReadingStrategy.BOUNDING_BOX;
     }
diff --git a/src/uk/ac/rdg/resc/ncwms/coords/HorizontalCoordSys.java 
b/src/uk/ac/rdg/resc/ncwms/coords/HorizontalCoordSys.java
index a3b3a6b..111b890 100644
--- a/src/uk/ac/rdg/resc/ncwms/coords/HorizontalCoordSys.java
+++ b/src/uk/ac/rdg/resc/ncwms/coords/HorizontalCoordSys.java
@@ -41,8 +41,8 @@ import ucar.unidata.geoloc.ProjectionImpl;
  */
 public abstract class HorizontalCoordSys
 {
-    private int xAxisSize;
-    private int yAxisSize;
+    protected int xAxisSize;
+    protected int yAxisSize;
 
     /** Protected constructor to limit direct instantiation to subclasses */
     protected HorizontalCoordSys() {}
diff --git a/src/uk/ac/rdg/resc/ncwms/coords/LatLonCoordSys.java 
b/src/uk/ac/rdg/resc/ncwms/coords/LatLonCoordSys.java
index 710a885..4ce9ceb 100644
--- a/src/uk/ac/rdg/resc/ncwms/coords/LatLonCoordSys.java
+++ b/src/uk/ac/rdg/resc/ncwms/coords/LatLonCoordSys.java
@@ -41,6 +41,8 @@ public final class LatLonCoordSys extends OneDCoordSys
     LatLonCoordSys(OneDCoordAxis lonAxis, OneDCoordAxis latAxis)
     {
         super(lonAxis, latAxis, null);
+        xAxisSize = lonAxis.getSize();
+        yAxisSize = latAxis.getSize();
     }
 
     /**
diff --git a/src/uk/ac/rdg/resc/ncwms/coords/PixelMap.java 
b/src/uk/ac/rdg/resc/ncwms/coords/PixelMap.java
index 2198e08..aabf317 100644
--- a/src/uk/ac/rdg/resc/ncwms/coords/PixelMap.java
+++ b/src/uk/ac/rdg/resc/ncwms/coords/PixelMap.java
@@ -74,19 +74,30 @@ public final class PixelMap implements 
Iterable<PixelMap.PixelMapEntry>
     private static final Logger logger = 
LoggerFactory.getLogger(PixelMap.class);
 
     /**
-     * <p>Maps points in the source grid to points in the target grid.  Each 
entry
-     * in this array represents a mapping from one source grid point (high four
-     * bytes) to one target grid point (low four bytes).</p>
-     * <p>We use an array of longs instead of an array of objects (e.g. 
Pair<Integer>
-     * or similar) because each object carries a memory overhead.</p>
+     * Source grid points corresponding to target grid points in
+     * {@link PixelMap#pixelMapTargetEntries}. Only the first
+     * {@link PixelMap#numEntries} elements are in use.
+     */
+    private long[] pixelMapSourceEntries;
+
+    /**
+     * Target grid points corresponding to source grid points in
+     * {@link PixelMap#pixelMapSourceEntries}. Only the first
+     * {@link PixelMap#numEntries} elements are in use.
+     */
+    private int[] pixelMapTargetEntries;
+
+    /**
+     * The number of entries in this pixel map, that is, the number of elements
+     * of {@link PixelMap#pixelMapSourceEntries} and
+     * {@link PixelMap#pixelMapTargetEntries} that are in use.
      */
-    private long[] pixelMapEntries;
-    /** The number of entries in this pixel map */
     private int numEntries = 0;
+
     /**
-     * The array of pixel map entries will grow in size as required by this
-     * number of longs.  This means that the array of pixel map entries will
-     * never be more than {@code chunkSize - 1} greater than it has to be.
+     * The arrays of pixel map entries will grow in size as required by this
+     * number of elements.  This means that the arrays of pixels map entries 
will
+     * never be more than {@code chunkSize - 1} greater than they have to be.
      */
     private final int chunkSize;
 
@@ -126,7 +137,8 @@ public final class PixelMap implements 
Iterable<PixelMap.PixelMapEntry>
         this.chunkSize = pointList.size() < 1000
             ? pointList.size()
             : pointList.size() / 10;
-        this.pixelMapEntries = new long[this.chunkSize];
+        this.pixelMapSourceEntries = new long[this.chunkSize];
+        this.pixelMapTargetEntries = new int[this.chunkSize];
 
         if (pointList instanceof HorizontalGrid)
         {
@@ -233,23 +245,40 @@ public final class PixelMap implements 
Iterable<PixelMap.PixelMapEntry>
         if (j < this.minJIndex) this.minJIndex = j;
         if (j > this.maxJIndex) this.maxJIndex = j;
 
-        // Calculate a single integer representing this grid point in the 
source grid
-        // TODO: watch out for overflows (would only happen with a very large 
grid!)
-        int sourceGridIndex = j * this.sourceGridISize + i;
+        // Calculate a single long integer representing this grid point in the 
source grid
+        long sourceGridIndex = ((long) j) * ((long) sourceGridISize) + ((long) 
i);
 
-        // See if we need to grow the array of pixel map entries
-        if (this.numEntries >= this.pixelMapEntries.length)
-        {
-            long[] newArray = new long[this.pixelMapEntries.length + 
this.chunkSize];
-            System.arraycopy(this.pixelMapEntries, 0, newArray, 0, 
this.pixelMapEntries.length);
-            this.pixelMapEntries = newArray;
+        if (pixelMapSourceEntries.length != pixelMapTargetEntries.length) {
+            throw new RuntimeException(
+                    "Internal error: array size mismatch in PixelMap "
+                            + "(this should never happen)");
+        }
+
+        // See if we need to grow the arrays of pixel map entries
+        if (numEntries >= pixelMapSourceEntries.length) {
+            // grow source array
+            {
+                long[] newPixelMapSourceEntries = new 
long[pixelMapSourceEntries.length
+                        + chunkSize];
+                System.arraycopy(this.pixelMapSourceEntries, 0,
+                        newPixelMapSourceEntries, 0,
+                        pixelMapSourceEntries.length);
+                pixelMapSourceEntries = newPixelMapSourceEntries;
+            }
+            // grow target array
+            {
+                int[] newPixelMapTargetEntries = new 
int[pixelMapTargetEntries.length
+                        + chunkSize];
+                System.arraycopy(pixelMapTargetEntries, 0,
+                        newPixelMapTargetEntries, 0,
+                        pixelMapTargetEntries.length);
+                pixelMapTargetEntries = newPixelMapTargetEntries;
+            }
         }
 
-        // Make an entry in the pixel map.  The source grid index becomes the
-        // high four bytes of the entry (long value) and the targetGridIndex
-        // becomes the low four bytes.
-        this.pixelMapEntries[this.numEntries] = (long)sourceGridIndex << 32 | 
targetGridIndex;
-        this.numEntries++;
+        pixelMapSourceEntries[numEntries] = sourceGridIndex;
+        pixelMapTargetEntries[numEntries] = targetGridIndex;
+        numEntries++;
     }
 
     /**
@@ -313,39 +342,21 @@ public final class PixelMap implements 
Iterable<PixelMap.PixelMapEntry>
     public int getNumUniqueIJPairs()
     {
         int count = 0;
-        for (PixelMapEntry pme : this) count++;
+        for (@SuppressWarnings("unused") PixelMapEntry pme : this) {
+            count++;
+        }
         return count;
     }
 
     /**
-     * Gets the sum of the lengths of each row of data points,
-     * {@literal i.e.} sum(imax - imin + 1).  This is the number of data 
points that will
-     * be extracted by the {@link DataReadingStrategy#SCANLINE SCANLINE} data
-     * reading strategy.
-     * @return the sum of the lengths of each row of data points
-     * @todo could reinstate this by moving the Scanline-generating code from
-     * DataReadingStrategy to this class and counting the lengths of the 
scanlines
-     */
-    /*public int getSumRowLengths()
-    {
-        int sumRowLengths = 0;
-        for (Row row : this.pixelMap.values())
-        {
-            sumRowLengths += (row.getMaxIIndex() - row.getMinIIndex() + 1);
-        }
-        return sumRowLengths;
-    }*/
-
-    /**
      * Gets the size of the i-j bounding box that encompasses all data.  This 
is
      * the number of data points that will be extracted using the
      * {@link DataReadingStrategy#BOUNDING_BOX BOUNDING_BOX} data reading 
strategy.
      * @return the size of the i-j bounding box that encompasses all data.
      */
-    public int getBoundingBoxSize()
-    {
-        return (this.maxIIndex - this.minIIndex + 1) *
-               (this.maxJIndex - this.minJIndex + 1);
+    public long getBoundingBoxSize() {
+        return ((long) (maxIIndex - minIIndex + 1))
+                * ((long) (maxJIndex - minJIndex + 1));
     }
 
     /**
@@ -356,32 +367,28 @@ public final class PixelMap implements 
Iterable<PixelMap.PixelMapEntry>
     {
         return new Iterator<PixelMapEntry>()
         {
-            /** Index in the array of entries */
+            /** Index in the arrays of entries */
             private int index = 0;
 
             @Override
             public boolean hasNext() {
-                return index < PixelMap.this.numEntries;
+                return index < numEntries;
             }
 
             @Override
             public PixelMapEntry next() {
-                long packed = PixelMap.this.pixelMapEntries[this.index];
-                // Unpack the source grid and target grid indices
-                int[] unpacked = unpack(packed);
-                this.index++;
-                final int sourceGridIndex = unpacked[0];
+                final long sourceGridIndex = pixelMapSourceEntries[index];
                 final List<Integer> targetGridIndices = new 
ArrayList<Integer>();
-                targetGridIndices.add(unpacked[1]);
+                targetGridIndices.add(pixelMapTargetEntries[index]);
+                index++;
 
                 // Now find all the other entries that use the same source grid
                 // index
                 boolean done = false;
-                while (!done && this.hasNext()) {
-                    packed = PixelMap.this.pixelMapEntries[this.index];
-                    unpacked = unpack(packed);
-                    if (unpacked[0] == sourceGridIndex) {
-                        targetGridIndices.add(unpacked[1]);
+                while (!done && hasNext()) {
+                    long otherSourceGridIndex = pixelMapSourceEntries[index];
+                    if (otherSourceGridIndex == sourceGridIndex) {
+                        targetGridIndices.add(pixelMapTargetEntries[index]);
                         this.index++;
                     } else {
                         done = true;
@@ -392,12 +399,12 @@ public final class PixelMap implements 
Iterable<PixelMap.PixelMapEntry>
 
                     @Override
                     public int getSourceGridIIndex() {
-                        return sourceGridIndex % PixelMap.this.sourceGridISize;
+                        return (int) (sourceGridIndex % sourceGridISize);
                     }
 
                     @Override
                     public int getSourceGridJIndex() {
-                        return sourceGridIndex / PixelMap.this.sourceGridISize;
+                        return (int) (sourceGridIndex / sourceGridISize);
                     }
 
                     @Override
@@ -406,6 +413,7 @@ public final class PixelMap implements 
Iterable<PixelMap.PixelMapEntry>
                     }
 
                 };
+
             }
 
             @Override
@@ -414,16 +422,7 @@ public final class PixelMap implements 
Iterable<PixelMap.PixelMapEntry>
             }
 
         };
-    }
 
-    /** Unpacks a long integer into two 4-byte integers (first value in array
-     * are the high 4 bytes of the long). */
-    private static int[] unpack(long packed)
-    {
-        return new int[] {
-            (int)(packed >> 32),
-            (int)(packed & 0xffffffff)
-        };
     }
 
 }
diff --git a/test/uk/ac/rdg/resc/ncwms/coords/PixelMapTest.java 
b/test/uk/ac/rdg/resc/ncwms/coords/PixelMapTest.java
new file mode 100644
index 0000000..b763f89
--- /dev/null
+++ b/test/uk/ac/rdg/resc/ncwms/coords/PixelMapTest.java
@@ -0,0 +1,277 @@
+package uk.ac.rdg.resc.ncwms.coords;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import junit.framework.Assert;
+
+import org.junit.Test;
+
+import ucar.nc2.constants.AxisType;
+import uk.ac.rdg.resc.ncwms.coords.PixelMap.PixelMapEntry;
+
+/**
+ * Tests for {@link PixelMap} index handling.
+ * 
+ * @author Ben Caradoc-Davies (CSIRO Earth Science and Resource Engineering)
+ */
+public class PixelMapTest {
+
+    /**
+     * Test for a small source grid.
+     * 
+     * @throws Exception
+     */
+    @Test
+    public void small() throws Exception {
+        test(SourceGrid.Small);
+    }
+
+    /**
+     * Test for a medium source grid.
+     * 
+     * @throws Exception
+     */
+    @Test
+    public void medium() throws Exception {
+        test(SourceGrid.Medium);
+    }
+
+    /**
+     * Test for a large source grid.
+     * 
+     * @throws Exception
+     */
+    @Test
+    public void large() throws Exception {
+        test(SourceGrid.Large);
+    }
+
+    /**
+     * Test that expected indices are found for a selection of arbitrary 
points.
+     */
+    @Test
+    public void indices() throws Exception {
+        testIndices(SourceGrid.Small, 0.51, 0.48, 2352, 7225);
+        testIndices(SourceGrid.Small, 0.01, 0.02, 46, 301);
+        testIndices(SourceGrid.Small, 0.97, 0.99, 4474, 14902);
+        testIndices(SourceGrid.Small, 0.03, 0.96, 138, 14451);
+        testIndices(SourceGrid.Small, 0.98, 0.05, 4520, 753);
+    }
+
+    /**
+     * Run tests for a source grid of a given size.
+     * 
+     * @param sourceGrid
+     * @throws Exception
+     */
+    public void test(SourceGrid sourceGrid) throws Exception {
+        HorizontalCoordSys cs = sourceGrid.buildCoordSys();
+        PointList targetPoints = sourceGrid.buildTargetPointList();
+        for (HorizontalPosition targetPoint : targetPoints.asList()) {
+            testPoint(cs, targetPoint);
+        }
+        testGrid(cs, targetPoints);
+    }
+
+    /**
+     * Test that a single point has consistent indices.
+     * 
+     * @param cs
+     * @param targetPoint
+     * @throws Exception
+     */
+    public void testPoint(HorizontalCoordSys cs, HorizontalPosition 
targetPoint)
+            throws Exception {
+        PixelMap pixelMap = new PixelMap(cs,
+                PointList.fromPoint((LonLatPosition) targetPoint));
+        Assert.assertEquals(pixelMap.getMinIIndex(), pixelMap.getMaxIIndex());
+        Assert.assertEquals(pixelMap.getMinJIndex(), pixelMap.getMaxJIndex());
+        int count = 0;
+        for (PixelMapEntry entry : pixelMap) {
+            Assert.assertEquals(pixelMap.getMinIIndex(),
+                    entry.getSourceGridIIndex());
+            Assert.assertEquals(pixelMap.getMinJIndex(),
+                    entry.getSourceGridJIndex());
+            count++;
+        }
+        Assert.assertEquals(1, count);
+    }
+
+    /**
+     * Test that expected indices are found for an arbitrary selection of
+     * points.
+     * 
+     * @param sourceGrid
+     * @param targetPointLonFraction
+     * @param targetPointLatFraction
+     * @param expectedSourceGridIIndex
+     * @param expectedSourceGridJIndex
+     * @throws Exception
+     */
+    public void testIndices(SourceGrid sourceGrid,
+            double targetPointLonFraction, double targetPointLatFraction,
+            int expectedSourceGridIIndex, int expectedSourceGridJIndex)
+            throws Exception {
+        PixelMap pixelMap = new PixelMap(sourceGrid.buildCoordSys(),
+                PointList.fromPoint((LonLatPosition) sourceGrid.buildTestPoint(
+                        targetPointLonFraction, targetPointLatFraction)));
+        Assert.assertEquals(pixelMap.getMinIIndex(), pixelMap.getMaxIIndex());
+        Assert.assertEquals(pixelMap.getMinJIndex(), pixelMap.getMaxJIndex());
+        int count = 0;
+        for (PixelMapEntry entry : pixelMap) {
+            Assert.assertEquals(pixelMap.getMinIIndex(),
+                    entry.getSourceGridIIndex());
+            Assert.assertEquals(pixelMap.getMinJIndex(),
+                    entry.getSourceGridJIndex());
+            Assert.assertEquals(expectedSourceGridIIndex,
+                    entry.getSourceGridIIndex());
+            Assert.assertEquals(expectedSourceGridJIndex,
+                    entry.getSourceGridJIndex());
+            count++;
+        }
+        Assert.assertEquals(1, count);
+    }
+
+    /**
+     * Test that grid point indices are in the expected row-major order.
+     * 
+     * @param cs
+     * @param targetPoints
+     * @throws Exception
+     */
+    public void testGrid(HorizontalCoordSys cs, PointList targetPoints)
+            throws Exception {
+        PixelMap pixelMap = new PixelMap(cs, targetPoints);
+        int count = 0;
+        int lastSourceGridIIndex = 0;
+        int lastSourceGridJIndex = 0;
+        for (PixelMapEntry entry : pixelMap) {
+            Assert.assertTrue(lastSourceGridIIndex <= entry
+                    .getSourceGridIIndex());
+            Assert.assertTrue(lastSourceGridJIndex < entry
+                    .getSourceGridJIndex()
+                    || lastSourceGridIIndex < entry.getSourceGridIIndex());
+            count++;
+            lastSourceGridIIndex = entry.getSourceGridIIndex();
+            lastSourceGridJIndex = entry.getSourceGridJIndex();
+        }
+        Assert.assertEquals(targetPoints.size(), count);
+    }
+
+    /**
+     * A longitude/latitude source grid.
+     */
+    public enum SourceGrid {
+
+        Small(141.84, 152.32, 4613, -44.20, -9.98, 15054), //
+        Medium(141.84, 152.32, 92255, -44.20, -9.98, 301081), //
+        Large(141.84, 152.32, 461301, -44.20, -9.98, 1505403);
+
+        private double lonMin;
+
+        private double lonMax;
+
+        private int lonSize;
+
+        private double latMin;
+
+        private double latMax;
+
+        private int latSize;
+
+        /**
+         * Constructor a source grid for given longitude and latitude range and
+         * grid size.
+         * 
+         * @param lonMin
+         * @param lonMax
+         * @param lonSize
+         * @param latMin
+         * @param latMax
+         * @param latSize
+         */
+        private SourceGrid(double lonMin, double lonMax, int lonSize,
+                double latMin, double latMax, int latSize) {
+            this.lonMin = lonMin;
+            this.lonMax = lonMax;
+            this.lonSize = lonSize;
+            this.latMin = latMin;
+            this.latMax = latMax;
+            this.latSize = latSize;
+        }
+
+        /**
+         * Return the coordinate system representing this source grid.
+         */
+        public HorizontalCoordSys buildCoordSys() {
+            return new LatLonCoordSys(buildCoordAxisAxis(lonMin, lonMax,
+                    lonSize, AxisType.Lon), buildCoordAxisAxis(latMin, latMax,
+                    latSize, AxisType.Lat));
+        }
+
+        /**
+         * Return a coordinate axis.
+         * 
+         * @param min
+         * @param max
+         * @param size
+         * @param axisType
+         * @return
+         */
+        private static OneDCoordAxis buildCoordAxisAxis(double min, double max,
+                int size, AxisType axisType) {
+            return new Regular1DCoordAxis(min, (max - min) / (size - 1), size,
+                    axisType);
+        }
+
+        /**
+         * Return a test point inside the limits of the grid calculated from
+         * fractions spans of longitude and latitude.
+         * 
+         * @param lonFraction
+         *            longitude fraction from 0.0 to 1.0
+         * @param latFraction
+         *            latitude fraction from 0.0 to 1.0)
+         * @return
+         */
+        public LonLatPosition buildTestPoint(double lonFraction,
+                double latFraction) {
+            return new LonLatPositionImpl(lonMin + lonFraction
+                    * (lonMax - lonMin), latMin + latFraction
+                    * (latMax - latMin));
+        }
+
+        /**
+         * Return an arbitrary target grid in row-major order. The target grid
+         * is contained in the source grid but does not necessarily share
+         * points. It is assumed to be much coarser than the source grid.
+         * 
+         * @return
+         */
+        public PointList buildTargetPointList() {
+            // these target grid parameters are arbitrary; should be coarser
+            // than source grid
+            int lonSize = 101;
+            int latSize = 103;
+            double lonFractionMin = 0.01;
+            double lonFractionMax = 0.98;
+            double latFractionMin = 0.03;
+            double latFractionMax = 0.95;
+            List<HorizontalPosition> positions = new 
ArrayList<HorizontalPosition>();
+            for (int i = 0; i < lonSize; i++) {
+                for (int j = 0; j < latSize; j++) {
+                    positions.add(buildTestPoint(
+                            lonFractionMin + i
+                                    * (lonFractionMax - lonFractionMin)
+                                    / (lonSize - 1), latFractionMin + j
+                                    * (latFractionMax - latFractionMin)
+                                    / (latSize - 1)));
+                }
+            }
+            return PointList.fromList(positions, CrsHelper.CRS_84);
+        }
+
+    }
+
+}

• Follow-Ups:
  • Re: [thredds] Thredds WMS support for large source grids
    • From: Jon Blower
  • Re: [thredds] Thredds WMS support for large source grids
    • From: Jon Blower