Skip to content

Commit

Permalink
adding various aggregation methods for raster merging + refactor
Browse files Browse the repository at this point in the history
  • Loading branch information
@moovida
moovida committed Nov 24, 2023
1 parent 0a80af7 commit b2c92b8
Show file tree
Hide file tree
Showing 7 changed files with 284 additions and 116 deletions.
6 changes: 4 additions & 2 deletions gears/src/main/java/org/hortonmachine/gears/io/stac/HMStacCollection.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
import org.geotools.stac.client.SearchQuery;
import org.hortonmachine.gears.libs.modules.HMRaster;
import org.hortonmachine.gears.libs.modules.HMRaster.HMRasterWritableBuilder;
import org.hortonmachine.gears.libs.modules.HMRaster.MergeMode;
import org.hortonmachine.gears.libs.monitor.DummyProgressMonitor;
import org.hortonmachine.gears.libs.monitor.IHMProgressMonitor;
import org.hortonmachine.gears.utils.CrsUtilities;
Expand Down Expand Up @@ -159,11 +160,12 @@ public List<HMStacItem> searchItems() throws Exception {
* @param bandName the name o the band to extract.
* @param items the list of items containing the various assets to read from.
* @param allowTransform if true, allows datasets of different projections to be transformed and merged together.
* @param mergeMode the merge mode to use in case of multiple value per cell.
* @return the final raster.
* @throws Exception
*/
public static HMRaster readRasterBandOnRegion( RegionMap latLongRegionMap, String bandName, List<HMStacItem> items,
boolean allowTransform, IHMProgressMonitor pm ) throws Exception {
boolean allowTransform, MergeMode mergeMode, IHMProgressMonitor pm ) throws Exception {

if (!allowTransform) {
List<String> epsgs = items.stream().map(( i ) -> i.getEpsg().toString()).distinct().collect(Collectors.toList());
Expand Down Expand Up @@ -217,7 +219,7 @@ public static HMRaster readRasterBandOnRegion( RegionMap latLongRegionMap, Strin
}

GridCoverage2D readRaster = asset.readRaster(readRegion);
outRaster.mapRasterSum(null, HMRaster.fromGridCoverage(readRaster));
outRaster.mapRaster(null, HMRaster.fromGridCoverage(readRaster), mergeMode);
pm.worked(1);
}
pm.done();
Expand Down
2 changes: 1 addition & 1 deletion gears/src/main/java/org/hortonmachine/gears/io/subsampling/OmsRasterSubsampler.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,7 @@ public void process() throws Exception {

HMRaster outHMRaster = new HMRaster.HMRasterWritableBuilder().setCrs(raster.getCrs()).setName("subsampled")
.setRegion(newRegionMap).setNoValue(raster.getNovalue()).build();
outHMRaster.mapRasterSubst(pm, raster);
outHMRaster.mapRaster(pm, raster, HMRaster.MergeMode.INSERT_ON_NOVALUE);

outRaster = outHMRaster.buildCoverage();
raster.close();
Expand Down
194 changes: 115 additions & 79 deletions gears/src/main/java/org/hortonmachine/gears/libs/modules/HMRaster.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,10 +31,12 @@
import org.geotools.coverage.grid.GridGeometry2D;
import org.geotools.geometry.jts.JTS;
import org.geotools.referencing.CRS;
import org.hortonmachine.gears.libs.exceptions.ModelsRuntimeException;
import org.hortonmachine.gears.libs.monitor.DummyProgressMonitor;
import org.hortonmachine.gears.libs.monitor.IHMProgressMonitor;
import org.hortonmachine.gears.utils.RegionMap;
import org.hortonmachine.gears.utils.coverage.CoverageUtilities;
import org.hortonmachine.gears.utils.math.NumericsUtilities;
import org.locationtech.jts.geom.Coordinate;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
Expand Down Expand Up @@ -62,10 +64,46 @@ public class HMRaster implements AutoCloseable {
private GridCoverage2D originalCoverage;

/**
* Raster map that can be initialized and used to count occurrences per cell.
* This count can then be used to average the summ of multiple occurrences.
* Support Rasters for the aggregation methods.
*/
private HMRaster sumRaster = null;
private HMRaster countRaster = null;
private CategoriesInCell[][] categoriesRaster = null;


/**
* Enumeration representing the merge modes for combining raster values.
*
* <p>These are used in the {@link #mapRaster(IHMProgressMonitor, HMRaster, MergeMode)} method.
*/
public static enum MergeMode {
/**
* Sum the values of the mapped rasters.
*/
SUM,

/**
* Average the values of the mapped rasters.
*/
AVG,

/**
* Substitute the values everytime. Last values wins.
*/
SUBSTITUTE,

/**
* Insert the values only if the cell contains novalue. First value wins.
*/
INSERT_ON_NOVALUE,

/**
* Collects all the values in the cell and keeps track of the post present one.
*
* <p>This requires the call of {@link #applyMostPopular()} to perform the proper substitution.
*/
MOST_POPULAR_VALUE
}

public static interface RasterCellProcessor {
void processCell( int col, int row, double value, int cols, int rows ) throws Exception;
Expand Down Expand Up @@ -399,46 +437,20 @@ public void close() throws Exception {
}
}

/**
* Writes the values of the coverage into the current raster, summing multiple occurrences.
*
* <p> In this case a bit matrix that tracks how many values per pixels are recorded is created
* and accessible through {@link #getCountRaster()}.
* Sine the values are summed, the number is needed for averaging.
*
* @param pm optional Process monitor.
* @param otherRaster the raster to map over the current raster.
* @throws IOException
*/
public void mapRasterSum( IHMProgressMonitor pm, HMRaster otherRaster ) throws Exception {
mapRaster(pm, otherRaster, 0);
}

/**
* Writes the values of the coverage into the current raster, only where the current does not have values.
*
* @param pm optional Process monitor.
* @param otherRaster the raster to map over the current raster.
* @throws IOException
*/
public void mapRasterSubst( IHMProgressMonitor pm, HMRaster otherRaster ) throws Exception {
mapRaster(pm, otherRaster, 1);
}

/**
* Writes the values of the coverage into the current raster, summing multiple occurrences.
*
* @param pm optional Process monitor.
* @param otherRaster the raster to map over the current raster.
* @param valuesCountRaster a bit matrix that tracks how many values per pixels are recorded (they are summed, so the number is needed for averaging).
* @param mergeMode optional merge mode parameter. null/0 = sum of both, 1 = place other only in novalues
* @param mergeMode optional merge mode parameter. If null MergeMode.SUM is used.
* @throws IOException
*/
public void mapRaster( IHMProgressMonitor pm, HMRaster otherRaster, Integer mergeMode ) throws Exception {
public void mapRaster( IHMProgressMonitor pm, HMRaster otherRaster, MergeMode mergeMode ) throws Exception {
if (pm == null)
pm = new DummyProgressMonitor();

int _mergeMode = 0;
MergeMode _mergeMode = MergeMode.SUM;
if (mergeMode != null) {
_mergeMode = mergeMode;
}
Expand Down Expand Up @@ -469,11 +481,16 @@ public void mapRaster( IHMProgressMonitor pm, HMRaster otherRaster, Integer merg
if (fromRow < 0)
fromRow = 0;

if (_mergeMode == 0 && countRaster == null) {
countRaster = new HMRasterWritableBuilder().setName("valuescount").setDoInteger(true).setRegion(regionMap).setCrs(crs)
.setInitialValue(0).build();
if (_mergeMode == MergeMode.AVG && sumRaster == null) {
// set the current sum and count to the current raster to start with
sumRaster = new HMRasterWritableBuilder().setTemplate(this).setCopyValues(true).setName("sum").build();
countRaster = new HMRasterWritableBuilder().setName("valuescount").setDoShort(true).setRegion(regionMap).setCrs(crs)
.setInitialValue((short) 1).build();
} else if (_mergeMode == MergeMode.MOST_POPULAR_VALUE && categoriesRaster == null) {
categoriesRaster = new CategoriesInCell[rows][cols];
}


pm.beginTask("Patch raster...", toRow - fromRow); //$NON-NLS-1$
// fill the points of the current raster picking form the
// other raster via nearest neighbor interpolation
Expand All @@ -483,23 +500,49 @@ public void mapRaster( IHMProgressMonitor pm, HMRaster otherRaster, Integer merg
Coordinate coordinate = getWorld(c, r);
double otherRasterValue = otherRaster.getValue(coordinate);
if (!otherRaster.isNovalue(otherRasterValue)) {
if (countRaster != null) {
int count = countRaster.getIntValue(c, r);
count++;
countRaster.setValue(c, r, count);
}
double thisRasterValue = getValue(c, r);

boolean thisIsNovalue = isNovalue(thisRasterValue);
if (thisIsNovalue) {
// if the current is novalue, then use the other
if (!thisIsNovalue && mergeMode == MergeMode.INSERT_ON_NOVALUE){
// value exists already, ignore the new one
continue;
}

if (_mergeMode == MergeMode.SUBSTITUTE || mergeMode == MergeMode.INSERT_ON_NOVALUE) {
// you want to always substitute, use the other
thisRasterValue = otherRasterValue;
} else if (_mergeMode == 0) {
// thisvalue is NOT novalue + mergemode is sum
} else if (_mergeMode == MergeMode.SUM) {
// just sum with any previous value
if (thisIsNovalue) {
thisRasterValue = 0;
}
thisRasterValue = thisRasterValue + otherRasterValue;
} else if (_mergeMode == MergeMode.AVG) {
double sum = sumRaster.getValue(c, r);
short count = countRaster.getShortValue(c, r);
if(sumRaster.isNovalue(sum)){
sum = 0;
count = 0;
}
sum += otherRasterValue;
sumRaster.setValue(c, r, sum);
count++;
countRaster.setValue(c, r, count);
thisRasterValue = sum / count;
} else if (_mergeMode == MergeMode.MOST_POPULAR_VALUE) {
CategoriesInCell categoriesInCell = categoriesRaster[r][c];
if(categoriesInCell == null){
categoriesInCell = new CategoriesInCell();
if(!thisIsNovalue){
categoriesInCell.addValue((int) thisRasterValue);
}
categoriesRaster[r][c] = categoriesInCell;
}
categoriesInCell.addValue((int) otherRasterValue);
thisRasterValue = categoriesInCell.mostPresentValue;
} else {
thisRasterValue = otherRasterValue;
// throw new ModelsRuntimeException("This should never happen.", this);
// thisRasterValue = otherRasterValue;
throw new ModelsRuntimeException("This should never happen.", this);
}
setValue(c, r, thisRasterValue);
}
Expand All @@ -510,40 +553,6 @@ public void mapRaster( IHMProgressMonitor pm, HMRaster otherRaster, Integer merg
pm.done();
}

/**
* @return the raster that contains the occurrences of data per cell.
*/
public HMRaster getCountRaster() {
return countRaster;
}

/**
* Apply averaging based on the {@link #countRaster} map. <b>Note that this modifies the original raster.</b>
*
* <p>This only makes sense if multiple mapping of the raster occurred (using {@link #mapRaster(IHMProgressMonitor, HMRaster, boolean)}.
*
* @param pm and optional progress monitor.
* @throws Exception
*/
public void applyCountAverage( IHMProgressMonitor pm ) throws Exception {
if (pm == null)
pm = new DummyProgressMonitor();

pm.beginTask("Averaging raster...", rows);
for( int r = 0; r < rows; r++ ) {
for( int c = 0; c < cols; c++ ) {
double value = getValue(c, r);
if (!isNovalue(value)) {
int count = countRaster.getIntValue(c, r);
double newValue = value / count;
setValue(c, r, newValue);
}
}
pm.worked(1);
}
pm.done();
}

/**
* Get the values surrounding the current col/row.
*
Expand Down Expand Up @@ -606,6 +615,33 @@ public void printData() {
}
}

private static class CategoriesInCell {
int mostPresentValue = 0;
private int valuesCount = 0;

private int arraySize = 3;
private int[] allValues = new int[arraySize];

void addValue(int value) {
allValues[valuesCount] = value;
valuesCount++;
if (valuesCount == arraySize) {
// we need to grow the array
int[] newArray = new int[arraySize + 2];
System.arraycopy(allValues, 0, newArray, 0, arraySize);
arraySize = arraySize + 2;
allValues = newArray;
}
if(valuesCount == 1){
mostPresentValue = value;
}else{
// find the most present value in the array
mostPresentValue = NumericsUtilities.getMostPopular(allValues, valuesCount);
}
}

}

public static class HMRasterWritableBuilder {
private String name = "newraster";

Expand Down
4 changes: 2 additions & 2 deletions gears/src/main/java/org/hortonmachine/gears/modules/r/mosaic/OmsMosaic.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@
import org.hortonmachine.gears.libs.modules.HMConstants;
import org.hortonmachine.gears.libs.modules.HMModel;
import org.hortonmachine.gears.libs.modules.HMRaster;
import org.hortonmachine.gears.libs.modules.HMRaster.MergeMode;
import org.hortonmachine.gears.utils.CrsUtilities;
import org.hortonmachine.gears.utils.RegionMap;
import org.hortonmachine.gears.utils.coverage.CoverageUtilities;
Expand Down Expand Up @@ -191,12 +192,11 @@ public void process() throws Exception {
HMRaster raster = data.getRaster();

pm.message("Patch map " + index++);
outHMRaster.mapRasterSum(pm, raster);
outHMRaster.mapRaster(pm, raster, MergeMode.AVG);

raster.close();

}
outHMRaster.applyCountAverage(pm);
outRaster = outHMRaster.buildCoverage();
outHMRaster.close();

Expand Down
27 changes: 27 additions & 0 deletions gears/src/main/java/org/hortonmachine/gears/utils/math/NumericsUtilities.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -426,4 +426,31 @@ public static double[] range2Bins( double min, double max, double step, boolean
return bins;
}

/**
* Finds the most popular element in an integer array.
*
* <p>If two values are equally popular, the first one is returned.</p>
*
* @param the array of integers to check.
* @param length the length of the array to check. Might be shorter than the array's size
* @return the most popular value.
*/
public static int getMostPopular(int[] a, int length) {
int maxCount = 1;
int popular = a[0];
for (int i = 0; i < (length - 1); i++) {
int value = a[i];
int valueCount = 1;
for (int j = 0; j < length; j++) {
if (i != j && value == a[j])
valueCount++;
}
if (valueCount > maxCount) {
popular = value;
maxCount = valueCount;
}
}
return popular;
}

}
Loading

0 comments on commit b2c92b8

Please sign in to comment.