Kontinuation commented on code in PR #1788: URL: https://github.com/apache/sedona/pull/1788#discussion_r1967047601
########## common/src/main/java/org/apache/sedona/common/raster/Rasterization.java: ########## @@ -0,0 +1,598 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ +package org.apache.sedona.common.raster; + +import java.awt.image.WritableRaster; +import java.util.*; +import javax.media.jai.RasterFactory; +import org.apache.sedona.common.Functions; +import org.apache.sedona.common.utils.RasterUtils; +import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.coverage.grid.GridCoverageFactory; +import org.geotools.geometry.Envelope2D; +import org.geotools.geometry.jts.JTS; +import org.locationtech.jts.geom.*; +import org.opengis.geometry.BoundingBox; +import org.opengis.referencing.FactoryException; + +public class Rasterization { + protected static List<Object> rasterize( + Geometry geom, + GridCoverage2D raster, + String pixelType, + double value, + boolean useGeometryExtent, + boolean allTouched) + throws FactoryException { + + // Validate the input geometry and raster metadata + double[] metadata = RasterAccessors.metadata(raster); + validateRasterMetadata(metadata); + if (!RasterPredicates.rsIntersects(raster, geom)) { + throw new IllegalArgumentException("Geometry does not intersect Raster."); + } + + Envelope2D rasterExtent = raster.getEnvelope2D(); + + Envelope2D geomExtent = rasterizeGeomExtent(geom, raster, metadata, allTouched); + + RasterizationParams params = + calculateRasterizationParams(raster, useGeometryExtent, metadata, geomExtent, pixelType); + + rasterizeGeometry(raster, metadata, geom, params, geomExtent, value, allTouched); + + // Create a GridCoverage2D for the rasterized result + GridCoverageFactory coverageFactory = new GridCoverageFactory(); + GridCoverage2D rasterized; + + if (useGeometryExtent) { + rasterized = coverageFactory.create("rasterized", params.writableRaster, geomExtent); + } else { + rasterized = coverageFactory.create("rasterized", params.writableRaster, rasterExtent); + } + + // Return results compatible with the original function + List<Object> objects = new ArrayList<>(); + objects.add(params.writableRaster); + objects.add(rasterized); + + return objects; + } + + private static void rasterizeGeometry( + GridCoverage2D raster, + double[] metadata, + Geometry geom, + RasterizationParams params, + Envelope2D geomExtent, + double value, + boolean allTouched) + throws FactoryException { + + switch (geom.getGeometryType()) { + case "GeometryCollection": + case "MultiPolygon": + case "MultiPoint": + rasterizeGeometryCollection(raster, metadata, geom, params, value, allTouched); + break; + case "Point": + rasterizePoint(geom, params, geomExtent, value); + break; + case "LineString": + case "MultiLineString": + case "LinearRing": + rasterizeLineString(geom, params, value, geomExtent); + break; + default: + rasterizePolygon(geom, params, geomExtent, value, allTouched); + break; + } + } + + private static void rasterizeGeometryCollection( + GridCoverage2D raster, + double[] metadata, + Geometry geom, + RasterizationParams params, + double value, + boolean allTouched) + throws FactoryException { + + for (int i = 0; i < geom.getNumGeometries(); i++) { + Geometry subGeom = geom.getGeometryN(i); + Envelope2D subGeomExtent = rasterizeGeomExtent(subGeom, raster, metadata, allTouched); + rasterizeGeometry(raster, metadata, subGeom, params, subGeomExtent, value, allTouched); + } + } + + private static void rasterizePoint( + Geometry geom, RasterizationParams params, Envelope2D geomExtent, double value) { + + int startX = (int) Math.round(((geomExtent.getMinX() - params.upperLeftX) / params.scaleX)); + int startY = (int) Math.round(((geomExtent.getMinY() - params.upperLeftY) / params.scaleY)); + int x = startX; + int y = startY; + + for (double worldY = geomExtent.getMinY(); + worldY < geomExtent.getMaxY(); + worldY += params.scaleY, y++) { + x = startX; + for (double worldX = geomExtent.getMinX(); + worldX < geomExtent.getMaxX(); + worldX += params.scaleX, x++) { + + // Flip y-axis (since raster Y starts from top-left) + int yIndex = -y - 1; + + // Create envelope for this pixel + double cellMaxX = worldX + params.scaleX; + double cellMaxY = worldY + params.scaleY; + + boolean intersects = + geom.intersects(JTS.toGeometry(new Envelope(worldX, cellMaxX, worldY, cellMaxY))); + + if (intersects) { + params.writableRaster.setSample(x, yIndex, 0, value); + } + } + } + } + + private static void rasterizeLineString( + Geometry geom, RasterizationParams params, double value, Envelope2D geomExtent) { + for (int i = 0; i < geom.getNumGeometries(); i++) { + LineString line = (LineString) geom.getGeometryN(i); + Coordinate[] coords = line.getCoordinates(); + + for (int j = 0; j < coords.length - 1; j++) { + // Extract start and end points for the segment + LineSegment clippedSegment = + clipSegmentToRasterBounds(coords[j], coords[j + 1], geomExtent); + + Coordinate start; + Coordinate end; + if (clippedSegment != null) { + start = new Coordinate(clippedSegment.p0.x, clippedSegment.p0.y); + end = new Coordinate(clippedSegment.p1.x, clippedSegment.p1.y); + } else { + continue; // Skip case where segment is completely outside geomExtent + } + + double x0 = (start.x - params.upperLeftX) / params.scaleX; + double y0 = (params.upperLeftY - start.y) / params.scaleY; + double x1 = (end.x - params.upperLeftX) / params.scaleX; + double y1 = (params.upperLeftY - end.y) / params.scaleY; + + // Apply Bresenham for this segment + drawLineBresenham(params, x0, y0, x1, y1, value, 0.2); + } + } + } + + // Modified Bresenham with Fractional Steps + private static void drawLineBresenham( + RasterizationParams params, + double x0, + double y0, + double x1, + double y1, + double value, + double stepSize) { Review Comment: I found that GDAL uses a [more complex approach](https://github.com/OSGeo/gdal/blob/v3.10.2/alg/llrasterize.cpp#L651-L710) for allTouched mode: it figures out which step to take (right or up/down) to step into the next pixel. This could be more reliable than using a fractional step of 0.2. I am OK with the current approach, but I think we can switch to GDAL's approach in the future. -- This is an automated message from the Apache Git Service. To respond to the message, please log on to GitHub and use the URL above to go to the specific comment. To unsubscribe, e-mail: [email protected] For queries about this service, please contact Infrastructure at: [email protected]
