Skip to content

Commit 62df2b9

Browse files
TonioFTonioF
authored
implemented coastal cloud distinction (#62)
1 parent b240a46 commit 62df2b9

File tree

2 files changed

+145
-0
lines changed

2 files changed

+145
-0
lines changed

‎idepix-s2msi/src/main/java/org/esa/snap/idepix/s2msi/CoastalCloudDistinction.java‎

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
package org.esa.snap.idepix.s2msi;
2+
3+
import org.esa.snap.core.dataio.ProductSubsetDef;
4+
import org.esa.snap.core.datamodel.Band;
5+
import org.esa.snap.core.datamodel.CrsGeoCoding;
6+
import org.esa.snap.core.datamodel.GeneralFilterBand;
7+
import org.esa.snap.core.datamodel.GeoCoding;
8+
import org.esa.snap.core.datamodel.Kernel;
9+
import org.esa.snap.core.datamodel.Mask;
10+
import org.esa.snap.core.datamodel.Product;
11+
import org.esa.snap.core.datamodel.ProductData;
12+
import org.esa.snap.core.datamodel.VirtualBand;
13+
import org.esa.snap.core.gpf.OperatorException;
14+
import org.esa.snap.core.gpf.Tile;
15+
import org.opengis.referencing.operation.MathTransform;
16+
17+
import java.awt.Color;
18+
import java.awt.geom.AffineTransform;
19+
import java.io.IOException;
20+
21+
import static org.esa.snap.idepix.s2msi.util.S2IdepixConstants.*;
22+
23+
public class CoastalCloudDistinction {
24+
25+
private static final float COAST_BUFFER_SIZE = 1000f;
26+
27+
private final Band dilatedLandFlag;
28+
private final Band dilatedWaterFlag;
29+
private final Band idx1band;
30+
private final Band idx2band;
31+
private final int resolution;
32+
33+
public CoastalCloudDistinction(Product s2ClassifProduct) {
34+
final Product ccdHelperProduct = createCcdHelperProduct(s2ClassifProduct);
35+
resolution = determineSourceResolution(s2ClassifProduct);
36+
dilatedLandFlag = createDilatedLandFlag(ccdHelperProduct);
37+
dilatedWaterFlag = createDilatedWaterFlag(ccdHelperProduct);
38+
idx1band = getRatio2_11(ccdHelperProduct);
39+
idx2band = getRatio8_11(ccdHelperProduct);
40+
}
41+
42+
private Product createCcdHelperProduct(Product s2ClassifProduct) {
43+
try {
44+
return s2ClassifProduct.createSubset(new ProductSubsetDef(), "ccd_helper", "");
45+
} catch (IOException e) {
46+
throw new OperatorException("Cannot create helper product for coastal cloud distinction." ,e);
47+
}
48+
}
49+
50+
private Band createDilatedLandFlag(Product ucdHelperProduct) {
51+
return createDilatedIdepixFlag(ucdHelperProduct, "tmp_land_mask", IDEPIX_LAND_DESCR_TEXT, IDEPIX_LAND_NAME);
52+
}
53+
54+
private Band createDilatedWaterFlag(Product ucdHelperProduct) {
55+
return createDilatedIdepixFlag(ucdHelperProduct, "tmp_water_mask", IDEPIX_WATER_DESCR_TEXT, IDEPIX_WATER_NAME);
56+
}
57+
58+
private Band createDilatedIdepixFlag(Product ucdHelperProduct,
59+
String maskName,
60+
String descriptionName,
61+
String flagName) {
62+
final Mask flagMask = Mask.BandMathsType.create(maskName, descriptionName,
63+
ucdHelperProduct.getSceneRasterWidth(), ucdHelperProduct.getSceneRasterHeight(),
64+
String.format("%1$s.%2$s", IDEPIX_CLASSIF_FLAGS, flagName),
65+
new Color(178, 0, 178), 0.5f);
66+
flagMask.setOwner(ucdHelperProduct);
67+
68+
return computeDilate(flagMask);
69+
}
70+
71+
private Band computeDilate(Band band) {
72+
final int kernelSize = (int) Math.floor((2 * COAST_BUFFER_SIZE) / resolution);
73+
final Kernel kernel = new Kernel(kernelSize, kernelSize, new double[kernelSize * kernelSize]);
74+
final GeneralFilterBand filterBand = new GeneralFilterBand(
75+
"__dilate33_" + band.getName(), band, GeneralFilterBand.OpType.DILATION, kernel, 1
76+
);
77+
band.getProduct().addBand(filterBand);
78+
return filterBand;
79+
}
80+
81+
public void correctCloudFlag(int x, int y, Tile classifFlagTile, Tile targetFlagTile) {
82+
final boolean landFlag = classifFlagTile.getSampleBit(x, y, IDEPIX_LAND);
83+
final boolean waterFlag = classifFlagTile.getSampleBit(x, y, IDEPIX_WATER);
84+
final boolean cloudAmbiguousFlag = classifFlagTile.getSampleBit(x, y, IDEPIX_CLOUD_AMBIGUOUS);
85+
final boolean cloudFlag = classifFlagTile.getSampleBit(x, y, IDEPIX_CLOUD);
86+
final boolean cloudSureFlag = classifFlagTile.getSampleBit(x, y, IDEPIX_CLOUD_SURE);
87+
targetFlagTile.setSample(x, y, classifFlagTile.getSampleInt(x, y));
88+
final boolean isCoastal = (dilatedLandFlag.getSampleInt(x, y) == 255 && !landFlag) ||
89+
(dilatedWaterFlag.getSampleInt(x, y) == 255 && !waterFlag);
90+
if (isCoastal) {
91+
final float idx1 = idx1band.getSampleFloat(x, y);
92+
final float idx2 = idx2band.getSampleFloat(x, y);
93+
final boolean notCoast = idx1 > 0.7 || (idx1 < 1 && idx1 > 0.6 && idx2 > 0.9);
94+
final boolean cloudSureCoastalFlag = cloudSureFlag && notCoast;
95+
final boolean cloudAmbiguousCoastalFlag = cloudAmbiguousFlag && notCoast;
96+
targetFlagTile.setSample(x, y, IDEPIX_CLOUD_AMBIGUOUS, cloudAmbiguousCoastalFlag);
97+
targetFlagTile.setSample(x, y, IDEPIX_CLOUD_SURE, cloudSureCoastalFlag);
98+
targetFlagTile.setSample(x, y, IDEPIX_CLOUD, cloudAmbiguousCoastalFlag || cloudSureCoastalFlag);
99+
} else {
100+
targetFlagTile.setSample(x, y, IDEPIX_CLOUD_AMBIGUOUS, cloudAmbiguousFlag);
101+
targetFlagTile.setSample(x, y, IDEPIX_CLOUD_SURE, cloudSureFlag);
102+
targetFlagTile.setSample(x, y, IDEPIX_CLOUD, cloudFlag);
103+
}
104+
}
105+
106+
private Band getRatio2_11(Product product) {
107+
return createRatioBand(product.getBand("B2"), product.getBand("B11"), product);
108+
}
109+
110+
private Band getRatio8_11(Product product) {
111+
return createRatioBand(product.getBand("B8"), product.getBand("B11"), product);
112+
}
113+
114+
private int determineSourceResolution(Product product) throws OperatorException {
115+
final GeoCoding sceneGeoCoding = product.getSceneGeoCoding();
116+
if (sceneGeoCoding instanceof CrsGeoCoding) {
117+
final MathTransform imageToMapTransform = sceneGeoCoding.getImageToMapTransform();
118+
if (imageToMapTransform instanceof AffineTransform) {
119+
return (int) ((AffineTransform) imageToMapTransform).getScaleX();
120+
}
121+
}
122+
throw new OperatorException("Invalid product");
123+
}
124+
125+
private static Band createRatioBand(Band numerator, Band denominator, Product owner) {
126+
final String name = String.format("__ratio_%s_%s", numerator.getName(), denominator.getName());
127+
final String expression = String.format("%s / %s", numerator.getName(), denominator.getName());
128+
return createVirtualBand(name, expression, owner);
129+
}
130+
131+
private static VirtualBand createVirtualBand(String name, String expression, Product owner) {
132+
final VirtualBand virtBand = new VirtualBand(name,
133+
ProductData.TYPE_FLOAT32,
134+
owner.getSceneRasterWidth(),
135+
owner.getSceneRasterHeight(),
136+
expression);
137+
owner.addBand(virtBand);
138+
return virtBand;
139+
}
140+
141+
}

‎idepix-s2msi/src/main/java/org/esa/snap/idepix/s2msi/operators/S2IdepixCloudPostProcessOp.java‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@
1212
import org.esa.snap.core.gpf.annotations.SourceProduct;
1313
import org.esa.snap.core.util.ProductUtils;
1414
import org.esa.snap.core.util.RectangleExtender;
15+
import org.esa.snap.idepix.s2msi.CoastalCloudDistinction;
1516
import org.esa.snap.idepix.s2msi.UrbanCloudDistinction;
1617
import org.esa.snap.idepix.s2msi.util.S2IdepixUtils;
1718

@@ -46,6 +47,7 @@ public class S2IdepixCloudPostProcessOp extends Operator {
4647
private Band origClassifFlagBand;
4748

4849
private RectangleExtender rectCalculator;
50+
private CoastalCloudDistinction coastalCloudDistinction;
4951
private UrbanCloudDistinction urbanCloudDistinction;
5052

5153

@@ -63,6 +65,7 @@ public void initialize() throws OperatorException {
6365
origClassifFlagBand = classifiedProduct.getBand(IDEPIX_CLASSIF_FLAGS);
6466
ProductUtils.copyBand(IDEPIX_CLASSIF_FLAGS, classifiedProduct, cloudBufferProduct, false);
6567

68+
coastalCloudDistinction = new CoastalCloudDistinction(classifiedProduct);
6669
urbanCloudDistinction = new UrbanCloudDistinction(classifiedProduct);
6770

6871
setTargetProduct(cloudBufferProduct);
@@ -95,6 +98,7 @@ public void computeTile(Band targetBand, Tile targetTile, ProgressMonitor pm) th
9598
boolean isCloud;
9699
if (targetRectangle.contains(x, y)) {
97100
S2IdepixUtils.combineFlags(x, y, sourceFlagTile, targetTile);
101+
coastalCloudDistinction.correctCloudFlag(x, y, targetTile, targetTile);
98102
urbanCloudDistinction.correctCloudFlag(x, y, targetTile, targetTile);
99103
isCloud = isCloudPixel(targetTile, y, x);
100104
} else {

0 commit comments

Comments
 (0)