Skip to content

Commit a616b0f

Browse files
committed
feat: make rasterization orientation-agnostic by supporting positive scaleY
1 parent b8eeb95 commit a616b0f

9 files changed

Lines changed: 860 additions & 29 deletions

File tree

common/src/main/java/org/apache/sedona/common/raster/Rasterization.java

Lines changed: 65 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,9 @@
1818
*/
1919
package org.apache.sedona.common.raster;
2020

21+
import static org.apache.sedona.common.utils.RasterUtils.getRasterScaleY;
22+
import static org.apache.sedona.common.utils.RasterUtils.getRasterUpperLeftY;
23+
2124
import java.awt.image.WritableRaster;
2225
import java.math.BigDecimal;
2326
import java.math.RoundingMode;
@@ -69,6 +72,11 @@ protected static List<Object> rasterize(
6972
rasterized = coverageFactory.create("rasterized", params.writableRaster, rasterExtent);
7073
}
7174

75+
// If original raster is bottom-up, flip the rasterized result vertically
76+
if (params.bottomUp) {
77+
rasterized = RasterUtils.flipVerticallyGridSpace(rasterized);
78+
}
79+
7280
// Return results compatible with the original function
7381
List<Object> objects = new ArrayList<>();
7482
objects.add(params.writableRaster);
@@ -136,16 +144,22 @@ private static void rasterizePoint(
136144
int x = startX;
137145
int y = startY;
138146

139-
for (double worldY = geomExtent.getMinY();
140-
worldY < geomExtent.getMaxY();
141-
worldY += params.scaleY, y++) {
147+
for (double worldY = geomExtent.getMaxY();
148+
worldY > geomExtent.getMinY();
149+
worldY += params.scaleY, y--) {
142150
x = startX;
143151
for (double worldX = geomExtent.getMinX();
144152
worldX < geomExtent.getMaxX();
145153
worldX += params.scaleX, x++) {
146154

147-
// Flip y-axis (since raster Y starts from top-left)
148-
int yIndex = -y - 1;
155+
// Make zero indexed
156+
int yIndex = y - 1;
157+
158+
// Adjust for bottom-up rasters
159+
// Reverse the y index
160+
if (params.bottomUp) {
161+
yIndex = params.writableRaster.getHeight() - 1 - yIndex;
162+
}
149163

150164
// Create envelope for this pixel
151165
double cellMaxX = worldX + params.scaleX;
@@ -182,9 +196,9 @@ private static void rasterizeLineString(
182196
}
183197

184198
double x0 = (start.x - params.upperLeftX) / params.scaleX;
185-
double y0 = (params.upperLeftY - start.y) / params.scaleY;
199+
double y0 = (start.y - params.upperLeftY) / params.scaleY;
186200
double x1 = (end.x - params.upperLeftX) / params.scaleX;
187-
double y1 = (params.upperLeftY - end.y) / params.scaleY;
201+
double y1 = (end.y - params.upperLeftY) / params.scaleY;
188202

189203
// Apply Bresenham for this segment
190204
drawLineBresenham(params, x0, y0, x1, y1, value, 0.2);
@@ -221,6 +235,12 @@ private static void drawLineBresenham(
221235
int rasterX = (int) (Math.floor(x));
222236
int rasterY = (int) (Math.floor(y));
223237

238+
// Adjust for bottom-up rasters
239+
// Reverse the y index
240+
if (params.bottomUp) {
241+
rasterY = params.writableRaster.getHeight() - 1 - rasterY;
242+
}
243+
224244
// Only write if within raster bounds
225245
if (rasterX >= 0
226246
&& rasterX < params.writableRaster.getWidth()
@@ -338,9 +358,9 @@ static ReferencedEnvelope rasterizeGeomExtent(
338358

339359
// Using BigDecimal to avoid floating point errors
340360
double upperLeftX = metadata[0];
341-
double upperLeftY = metadata[1];
361+
double upperLeftY = getRasterUpperLeftY(metadata);
342362
double scaleX = metadata[4];
343-
double scaleY = metadata[5];
363+
double scaleY = getRasterScaleY(metadata);
344364

345365
// Compute the aligned min/max values
346366
double alignedMinX =
@@ -464,13 +484,14 @@ private static RasterizationParams calculateRasterizationParams(
464484
upperLeftY = geomExtent.getMaxY();
465485
} else {
466486
upperLeftX = metadata[0];
467-
upperLeftY = metadata[1];
487+
upperLeftY = getRasterUpperLeftY(metadata);
468488
}
469489

470490
WritableRaster writableRaster;
471491
if (useGeometryExtent) {
472492
int geomExtentWidth = (int) (Math.round(geomExtent.getWidth() / metadata[4]));
473-
int geomExtentHeight = (int) (Math.round(geomExtent.getHeight() / -metadata[5]));
493+
int geomExtentHeight =
494+
(int) (Math.round(geomExtent.getHeight() / -getRasterScaleY(metadata)));
474495

475496
writableRaster =
476497
RasterFactory.createBandedRaster(
@@ -483,19 +504,30 @@ private static RasterizationParams calculateRasterizationParams(
483504
RasterUtils.getDataTypeCode(pixelType), rasterWidth, rasterHeight, 1, null);
484505
}
485506

507+
boolean bottomUp = metadata[5] > 0;
508+
486509
return new RasterizationParams(
487-
writableRaster, pixelType, metadata[4], -metadata[5], upperLeftX, upperLeftY);
510+
writableRaster,
511+
pixelType,
512+
metadata[4],
513+
getRasterScaleY(metadata),
514+
upperLeftX,
515+
upperLeftY,
516+
bottomUp);
488517
}
489518

490519
private static void validateRasterMetadata(double[] rasterMetadata) {
491520
if (rasterMetadata[4] < 0) {
492-
throw new IllegalArgumentException("ScaleX cannot be negative");
493-
}
494-
if (rasterMetadata[5] > 0) {
495-
throw new IllegalArgumentException("ScaleY must be negative");
521+
throw new IllegalArgumentException(
522+
String.format(
523+
"Invalid raster metadata: scaleX is negative (%.2f). Right-to-left rasters are not supported.",
524+
rasterMetadata[4]));
496525
}
497526
if (rasterMetadata[6] != 0 || rasterMetadata[7] != 0) {
498-
throw new IllegalArgumentException("SkewX and SkewY must be zero");
527+
throw new IllegalArgumentException(
528+
String.format(
529+
"Invalid raster metadata: skewX is %.2f and skewY is %.2f. Both values must be zero.",
530+
rasterMetadata[6], rasterMetadata[7]));
499531
}
500532
}
501533

@@ -506,20 +538,23 @@ private static class RasterizationParams {
506538
double scaleY;
507539
double upperLeftX;
508540
double upperLeftY;
541+
boolean bottomUp;
509542

510543
RasterizationParams(
511544
WritableRaster writableRaster,
512545
String pixelType,
513546
double scaleX,
514547
double scaleY,
515548
double upperLeftX,
516-
double upperLeftY) {
549+
double upperLeftY,
550+
boolean bottomUp) {
517551
this.writableRaster = writableRaster;
518552
this.pixelType = pixelType;
519553
this.scaleX = scaleX;
520554
this.scaleY = scaleY;
521555
this.upperLeftX = upperLeftX;
522556
this.upperLeftY = upperLeftY;
557+
this.bottomUp = bottomUp;
523558
}
524559
}
525560

@@ -602,15 +637,15 @@ private static Map<Double, TreeSet<Double>> computeScanlineIntersections(
602637
// Using BigDecimal to avoid floating point errors
603638
double yStart =
604639
Math.round(
605-
(BigDecimal.valueOf(params.upperLeftY)
606-
.subtract(BigDecimal.valueOf(worldP1.y))
640+
(BigDecimal.valueOf(worldP1.y)
641+
.subtract(BigDecimal.valueOf(params.upperLeftY))
607642
.divide(BigDecimal.valueOf(params.scaleY), RoundingMode.CEILING))
608643
.doubleValue());
609644

610645
double yEnd =
611646
Math.round(
612-
(BigDecimal.valueOf(params.upperLeftY)
613-
.subtract(BigDecimal.valueOf(worldP2.y))
647+
(BigDecimal.valueOf(worldP2.y)
648+
.subtract(BigDecimal.valueOf(params.upperLeftY))
614649
.divide(BigDecimal.valueOf(params.scaleY), RoundingMode.FLOOR))
615650
.doubleValue());
616651

@@ -619,7 +654,7 @@ private static Map<Double, TreeSet<Double>> computeScanlineIntersections(
619654
yStart = Math.min((params.writableRaster.getHeight() - 0.5), Math.abs(yStart) - 0.5);
620655

621656
double p1X = (worldP1.x - params.upperLeftX) / params.scaleX;
622-
double p1Y = (params.upperLeftY - worldP1.y) / params.scaleY;
657+
double p1Y = (worldP1.y - params.upperLeftY) / params.scaleY;
623658

624659
if (worldP1.x == worldP2.x) {
625660
// Vertical line case: directly set xIntercept. Avoid divide by zero error when
@@ -695,6 +730,13 @@ private static void fillPolygon(
695730

696731
for (Map.Entry<Integer, List<int[]>> entry : scanlineFillRanges.entrySet()) {
697732
int y = entry.getKey();
733+
734+
// Adjust for bottom-up rasters
735+
// Reverse the y index
736+
if (params.bottomUp) {
737+
y = params.writableRaster.getHeight() - 1 - y;
738+
}
739+
698740
for (int[] range : entry.getValue()) {
699741
if (range.length == 1) {
700742
params.writableRaster.setSample(range[0], y, 0, value);

common/src/main/java/org/apache/sedona/common/utils/RasterUtils.java

Lines changed: 166 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -984,4 +984,170 @@ public static GridCoverage2D applyRasterMask(GridCoverage2D raster, GridCoverage
984984
false);
985985
return modifiedRaster;
986986
}
987+
988+
/**
989+
* This function calculates the upper-left Y world coordinate of the raster. It accounts for both
990+
* bottom-up and top-down raster conventions.
991+
*
992+
* <p>For bottom-up rasters (metadata[5] > 0), the upper-left Y world coordinate is directly taken
993+
* from metadata[1]. - For top-down rasters (metadata[5] < 0), the upper-left Y coordinate is
994+
* adjusted by adding the product of the raster height and the scale factor (metadata[5]).
995+
*
996+
* @param metadata The metadata array containing raster properties.
997+
* @return The upper-left Y world coordinate of the raster.
998+
*/
999+
public static double getRasterUpperLeftY(double[] metadata) {
1000+
if (metadata[5] < 0) {
1001+
return metadata[1];
1002+
} else {
1003+
int rasterHeight = (int) metadata[3];
1004+
return metadata[1] + (rasterHeight * metadata[5]);
1005+
}
1006+
}
1007+
1008+
/**
1009+
* This function calculates the Y scale factor of the raster. It ensures that the scale factor is
1010+
* always negative, regardless of the raster convention.
1011+
*
1012+
* <p>For bottom-up rasters (metadata[5] > 0), the scale factor is returned as-is. - For top-down
1013+
* rasters (metadata[5] < 0), the scale factor is negated to maintain consistency.
1014+
*
1015+
* @param metadata The metadata array containing raster properties.
1016+
* @return The Y scale factor of the raster.
1017+
*/
1018+
public static double getRasterScaleY(double[] metadata) {
1019+
if (metadata[5] < 0) {
1020+
return metadata[5];
1021+
} else {
1022+
return -metadata[5];
1023+
}
1024+
}
1025+
1026+
/**
1027+
* GRID SPACE: Grid space refers to the raster's index space (col, row). Pixel values remain at
1028+
* the same array indices, but the grid-to-CRS affine transform is rewritten.
1029+
*
1030+
* <p>EFFECT: - The raster is reinterpreted spatially (top-down ↔ bottom-up). - Pixel values are
1031+
* NOT moved in memory. - Pixel world coordinates DO change. - The world envelope (bounding box)
1032+
* is preserved.
1033+
*
1034+
* <p>Use this when you want to change how the raster is oriented in CRS space without touching
1035+
* the pixel data.
1036+
*/
1037+
public static GridCoverage2D flipVerticallyGridSpace(GridCoverage2D raster) {
1038+
1039+
GridGeometry2D gridGeometry = raster.getGridGeometry();
1040+
1041+
MathTransform mt = gridGeometry.getGridToCRS2D();
1042+
if (!(mt instanceof AffineTransform2D)) {
1043+
throw new UnsupportedOperationException(
1044+
"flipVertically only supports AffineTransform2D grid-to-CRS transforms");
1045+
}
1046+
1047+
AffineTransform2D affine = (AffineTransform2D) mt;
1048+
1049+
RenderedImage image = raster.getRenderedImage();
1050+
int height = image.getHeight();
1051+
1052+
/*
1053+
* Affine transform definition (GeoTools / GDAL compatible):
1054+
*
1055+
* worldX = scaleX * col + shearX * row + translateX
1056+
* worldY = shearY * col + scaleY * row + translateY
1057+
*/
1058+
double scaleX = affine.getScaleX();
1059+
double shearX = affine.getShearX();
1060+
double shearY = affine.getShearY();
1061+
double scaleY = affine.getScaleY();
1062+
double translateX = affine.getTranslateX();
1063+
double translateY = affine.getTranslateY();
1064+
1065+
/*
1066+
* Apply grid-space vertical flip:
1067+
*
1068+
* row = (H - 1) - row'
1069+
*/
1070+
double newShearX = -shearX;
1071+
double newScaleY = -scaleY;
1072+
1073+
double newTranslateX = translateX + shearX * (height - 1);
1074+
double newTranslateY = translateY + scaleY * (height - 1);
1075+
1076+
AffineTransform2D flippedAffine =
1077+
new AffineTransform2D(scaleX, shearY, newShearX, newScaleY, newTranslateX, newTranslateY);
1078+
1079+
GridGeometry2D newGridGeometry =
1080+
new GridGeometry2D(
1081+
gridGeometry.getGridRange(),
1082+
flippedAffine,
1083+
gridGeometry.getCoordinateReferenceSystem());
1084+
1085+
return RasterUtils.clone(
1086+
raster.getRenderedImage(),
1087+
newGridGeometry,
1088+
raster.getSampleDimensions(),
1089+
raster,
1090+
null,
1091+
true);
1092+
}
1093+
1094+
/**
1095+
* PIXEL SPACE: Pixel space refers to pixel values anchored at their world coordinates. Pixel
1096+
* values are physically rearranged so that each value remains at the same CRS location.
1097+
*
1098+
* <p>EFFECT: - Raster rows are reversed in memory. - The affine transform is updated accordingly.
1099+
* - Pixel world coordinates are preserved. - The world envelope (bounding box) is preserved.
1100+
*
1101+
* <p>Use this when pixel values must remain fixed at their geographic locations.
1102+
*/
1103+
public static GridCoverage2D flipVerticallyPixelSpace(GridCoverage2D raster) {
1104+
GridGeometry2D gridGeometry = raster.getGridGeometry();
1105+
1106+
MathTransform mt = gridGeometry.getGridToCRS2D();
1107+
if (!(mt instanceof AffineTransform2D)) {
1108+
throw new UnsupportedOperationException(
1109+
"flipVertically only supports AffineTransform2D grid-to-CRS transforms");
1110+
}
1111+
1112+
AffineTransform2D affine = (AffineTransform2D) mt;
1113+
1114+
RenderedImage image = raster.getRenderedImage();
1115+
int height = image.getHeight();
1116+
1117+
// Flip the affine transformation
1118+
double scaleX = affine.getScaleX();
1119+
double shearX = affine.getShearX();
1120+
double shearY = affine.getShearY();
1121+
double scaleY = affine.getScaleY();
1122+
double translateX = affine.getTranslateX();
1123+
double translateY = affine.getTranslateY();
1124+
1125+
double newShearX = -shearX;
1126+
double newScaleY = -scaleY;
1127+
double newTranslateX = translateX + shearX * (height - 1);
1128+
double newTranslateY = translateY + scaleY * (height - 1);
1129+
1130+
AffineTransform2D flippedAffine =
1131+
new AffineTransform2D(scaleX, shearY, newShearX, newScaleY, newTranslateX, newTranslateY);
1132+
1133+
// Create a new WritableRaster to preserve pixel values
1134+
Raster originalRaster = RasterUtils.getRaster(image);
1135+
WritableRaster flippedRaster = originalRaster.createCompatibleWritableRaster();
1136+
1137+
// Copy pixel values to the flipped raster
1138+
for (int y = 0; y < height; y++) {
1139+
for (int x = 0; x < image.getWidth(); x++) {
1140+
flippedRaster.setPixel(x, height - y - 1, originalRaster.getPixel(x, y, (double[]) null));
1141+
}
1142+
}
1143+
1144+
GridGeometry2D newGridGeometry =
1145+
new GridGeometry2D(
1146+
gridGeometry.getGridRange(),
1147+
flippedAffine,
1148+
gridGeometry.getCoordinateReferenceSystem());
1149+
1150+
return RasterUtils.clone(
1151+
flippedRaster, newGridGeometry, raster.getSampleDimensions(), raster, null, true);
1152+
}
9871153
}

0 commit comments

Comments
 (0)