This is an automated email from the ASF dual-hosted git repository. jiayu pushed a commit to branch getband in repository https://gitbox.apache.org/repos/asf/sedona.git
commit f76b2071ef1eb15dbf9a21d951a4370ef94bf575 Author: Jia Yu <[email protected]> AuthorDate: Mon Jun 12 03:57:17 2023 -0700 Add functions --- .../apache/sedona/common/raster/Constructors.java | 38 ------ .../org/apache/sedona/common/raster/Functions.java | 125 ------------------ .../apache/sedona/common/raster/MapAlgebra.java | 147 +++++++++++++++++++++ .../sedona/common/raster/PixelFunctions.java | 86 ++++++++++++ .../sedona/common/raster/RasterAccessors.java | 102 ++++++++++++++ .../sedona/common/raster/RasterConstructors.java | 94 +++++++++++++ .../apache/sedona/common/raster/RasterEditors.java | 44 ++++++ .../raster/{Outputs.java => RasterOutputs.java} | 2 +- .../org/apache/sedona/common/raster/Serde.java | 7 +- .../apache/sedona/common/raster/FunctionsTest.java | 65 ++++----- .../sedona/common/raster/MapAlgebraTest.java | 147 +++++++++++++++++++++ .../sedona/common/raster/RasterAccessorsTest.java | 101 ++++++++++++++ ...uctorsTest.java => RasterConstructorsTest.java} | 39 +++++- .../sedona/common/raster/RasterTestBase.java | 12 +- .../org/apache/sedona/common/raster/SerdeTest.java | 3 +- docs/api/sql/Raster-loader.md | 49 +++++++ docs/api/sql/Raster-operators.md | 101 +++++++++++++- docs/api/sql/Raster-writer.md | 3 + .../scala/org/apache/sedona/sql/UDF/Catalog.scala | 4 + .../raster/{Functions.scala => MapAlgebra.scala} | 144 ++++---------------- .../expressions/raster/PixelFunctions.scala | 85 ++++++++++++ .../expressions/raster/RasterAccessors.scala | 121 +++++++++++++++++ ...Constructors.scala => RasterConstructors.scala} | 43 +++++- .../{Constructors.scala => RasterEditors.scala} | 51 ++----- .../raster/{Outputs.scala => RasterOutputs.scala} | 12 +- .../org/apache/sedona/sql/rasteralgebraTest.scala | 53 ++++++++ .../org/apache/sedona/sql/serdeAwareTest.scala | 2 +- 27 files changed, 1287 insertions(+), 393 deletions(-) diff --git a/common/src/main/java/org/apache/sedona/common/raster/Constructors.java b/common/src/main/java/org/apache/sedona/common/raster/Constructors.java deleted file mode 100644 index e39a80b7..00000000 --- a/common/src/main/java/org/apache/sedona/common/raster/Constructors.java +++ /dev/null @@ -1,38 +0,0 @@ -/** - * Licensed 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 - * <p> - * http://www.apache.org/licenses/LICENSE-2.0 - * <p> - * 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 org.geotools.coverage.grid.GridCoverage2D; -import org.geotools.gce.arcgrid.ArcGridReader; -import org.geotools.gce.geotiff.GeoTiffReader; - -import java.io.ByteArrayInputStream; -import java.io.IOException; - -public class Constructors { - - public static GridCoverage2D fromArcInfoAsciiGrid(byte[] bytes) throws IOException { - try (ByteArrayInputStream inputStream = new ByteArrayInputStream(bytes)) { - ArcGridReader reader = new ArcGridReader(inputStream); - return reader.read(null); - } - } - - public static GridCoverage2D fromGeoTiff(byte[] bytes) throws IOException { - try (ByteArrayInputStream inputStream = new ByteArrayInputStream(bytes)) { - GeoTiffReader geoTiffReader = new GeoTiffReader(inputStream); - return geoTiffReader.read(null); - } - } -} diff --git a/common/src/main/java/org/apache/sedona/common/raster/Functions.java b/common/src/main/java/org/apache/sedona/common/raster/Functions.java deleted file mode 100644 index 809f9b33..00000000 --- a/common/src/main/java/org/apache/sedona/common/raster/Functions.java +++ /dev/null @@ -1,125 +0,0 @@ -/** - * Licensed 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 - * <p> - * http://www.apache.org/licenses/LICENSE-2.0 - * <p> - * 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 org.geotools.coverage.CoverageFactoryFinder; -import org.geotools.coverage.grid.GridCoordinates2D; -import org.geotools.coverage.grid.GridCoverage2D; -import org.geotools.coverage.grid.GridCoverageFactory; -import org.geotools.coverage.grid.GridGeometry2D; -import org.geotools.gce.geotiff.GeoTiffWriter; -import org.geotools.geometry.DirectPosition2D; -import org.geotools.geometry.Envelope2D; -import org.geotools.geometry.jts.ReferencedEnvelope; -import org.geotools.referencing.CRS; -import org.geotools.referencing.crs.DefaultEngineeringCRS; -import org.locationtech.jts.geom.*; -import org.opengis.referencing.FactoryException; -import org.opengis.referencing.crs.CoordinateReferenceSystem; -import org.opengis.referencing.operation.TransformException; - -import java.awt.image.Raster; -import java.util.ArrayList; -import java.util.Collections; -import java.util.List; -import java.util.Optional; -import java.util.function.DoublePredicate; -import java.util.stream.Collectors; -import java.util.stream.DoubleStream; - -public class Functions { - - public static Geometry envelope(GridCoverage2D raster) throws FactoryException { - Envelope2D envelope2D = raster.getEnvelope2D(); - - Envelope envelope = new Envelope(envelope2D.getMinX(), envelope2D.getMaxX(), envelope2D.getMinY(), envelope2D.getMaxY()); - int srid = srid(raster); - return new GeometryFactory(new PrecisionModel(), srid).toGeometry(envelope); - } - - public static int numBands(GridCoverage2D raster) { - return raster.getNumSampleDimensions(); - } - - public static GridCoverage2D setSrid(GridCoverage2D raster, int srid) throws FactoryException { - CoordinateReferenceSystem crs; - if (srid == 0) { - crs = DefaultEngineeringCRS.CARTESIAN_2D; - } else { - crs = CRS.decode("EPSG:" + srid); - } - ReferencedEnvelope referencedEnvelope = new ReferencedEnvelope(raster.getEnvelope2D(), crs); - GridCoverageFactory gridCoverageFactory = CoverageFactoryFinder.getGridCoverageFactory(null); - return gridCoverageFactory.create(raster.getName().toString(), raster.getRenderedImage(), referencedEnvelope); - } - - public static int srid(GridCoverage2D raster) throws FactoryException { - CoordinateReferenceSystem crs = raster.getCoordinateReferenceSystem(); - if (crs instanceof DefaultEngineeringCRS) { - // GeoTools defaults to internal non-standard epsg codes, like 404000, if crs is missing. - // We need to check for this case and return 0 instead. - if (((DefaultEngineeringCRS) crs).isWildcard()) { - return 0; - } - } - return Optional.ofNullable(CRS.lookupEpsgCode(crs, true)).orElse(0); - } - - public static Double value(GridCoverage2D rasterGeom, Geometry geometry, int band) throws TransformException { - return values(rasterGeom, Collections.singletonList(geometry), band).get(0); - } - - public static List<Double> values(GridCoverage2D rasterGeom, List<Geometry> geometries, int band) throws TransformException { - int numBands = rasterGeom.getNumSampleDimensions(); - if (band < 1 || band > numBands) { - // Invalid band index. Return nulls. - return geometries.stream().map(geom -> (Double) null).collect(Collectors.toList()); - } - Raster raster = rasterGeom.getRenderedImage().getData(); - GridGeometry2D gridGeometry = rasterGeom.getGridGeometry(); - double[] noDataValues = rasterGeom.getSampleDimension(band - 1).getNoDataValues(); - DoublePredicate isNoData = d -> noDataValues != null && DoubleStream.of(noDataValues).anyMatch(noDataValue -> Double.compare(noDataValue, d) == 0); - double[] pixelBuffer = new double[numBands]; - - List<Double> result = new ArrayList<>(geometries.size()); - for (Geometry geom : geometries) { - if (geom == null) { - result.add(null); - } else { - Point point = ensurePoint(geom); - DirectPosition2D directPosition2D = new DirectPosition2D(point.getX(), point.getY()); - GridCoordinates2D gridCoordinates2D = gridGeometry.worldToGrid(directPosition2D); - try { - double pixel = raster.getPixel(gridCoordinates2D.x, gridCoordinates2D.y, pixelBuffer)[band - 1]; - if (isNoData.test(pixel)) { - result.add(null); - } else { - result.add(pixel); - } - } catch (ArrayIndexOutOfBoundsException exc) { - // Points outside the extent should return null - result.add(null); - } - } - } - return result; - } - - private static Point ensurePoint(Geometry geometry) { - if (geometry instanceof Point) { - return (Point) geometry; - } - throw new IllegalArgumentException("Attempting to get the value of a pixel with a non-point geometry."); - } -} diff --git a/common/src/main/java/org/apache/sedona/common/raster/MapAlgebra.java b/common/src/main/java/org/apache/sedona/common/raster/MapAlgebra.java new file mode 100644 index 00000000..0a24cce2 --- /dev/null +++ b/common/src/main/java/org/apache/sedona/common/raster/MapAlgebra.java @@ -0,0 +1,147 @@ +/* + * 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 org.geotools.coverage.CoverageFactoryFinder; +import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.coverage.grid.GridCoverageFactory; + +import javax.media.jai.RasterFactory; + +import java.awt.Point; +import java.awt.image.Raster; +import java.awt.image.RenderedImage; +import java.awt.image.WritableRaster; + +public class MapAlgebra +{ + /** + * Returns the values of the given band as a 1D array. + * @param rasterGeom + * @param bandIndex starts at 1 + * @return + */ + public static double[] bandAsArray(GridCoverage2D rasterGeom, int bandIndex) { + int numBands = rasterGeom.getNumSampleDimensions(); + if (bandIndex < 1 || bandIndex > numBands) { + // Invalid band index. Return nulls. + return null; + } + + Raster raster = rasterGeom.getRenderedImage().getData(); + // Get the width and height of the raster + int width = raster.getWidth(); + int height = raster.getHeight(); + + // Array to hold the band values + double[] bandValues = new double[width * height]; + + // Iterate over each pixel + for (int y = 0; y < height; y++) { + for (int x = 0; x < width; x++) { + // Get the value and store it in the 1D array + bandValues[y * width + x] = raster.getSample(x, y, bandIndex - 1); + } + } + return bandValues; + } + + /** + * Adds a new band to the given raster, using the given array as the band values. + * @param rasterGeom + * @param bandValues + * @param bandIndex starts at 1, and no larger than numBands + 1 + * @return + */ + public static GridCoverage2D addBandFromArray(GridCoverage2D rasterGeom, double[] bandValues, int bandIndex) { + int numBands = rasterGeom.getNumSampleDimensions(); + // Allow the band index to be one larger than the number of bands, which will append the band to the end + if (bandIndex < 1 || bandIndex > numBands + 1) { + throw new IllegalArgumentException("Band index is out of bounds. Must be between 1 and " + (numBands + 1) + ")"); + } + + if (bandIndex == numBands + 1) { + return copyRasterAndAppendBand(rasterGeom, bandValues); + } + else { + return copyRasterAndReplaceBand(rasterGeom, bandIndex, bandValues); + } + } + + + /** + * Adds a new band to the given raster, using the given array as the band values. The new band is appended to the end. + * @param rasterGeom + * @param bandValues + * @return + */ + public static GridCoverage2D addBandFromArray(GridCoverage2D rasterGeom, double[] bandValues) { + return addBandFromArray(rasterGeom, bandValues, rasterGeom.getNumSampleDimensions() + 1); + } + + /** + * This is an experimental method as it does not copy the original raster properties (e.g. color model, sample model, etc.) + * TODO: Copy the original raster properties + * @param gridCoverage2D + * @param bandValues + * @return + */ + private static GridCoverage2D copyRasterAndAppendBand(GridCoverage2D gridCoverage2D, double[] bandValues) { + // Get the original image and its properties + RenderedImage originalImage = gridCoverage2D.getRenderedImage(); + Raster raster = originalImage.getData(); + Point location = gridCoverage2D.getRenderedImage().getData().getBounds().getLocation(); + WritableRaster wr = RasterFactory.createBandedRaster(raster.getDataBuffer().getDataType(), originalImage.getWidth(), originalImage.getHeight(), gridCoverage2D.getNumSampleDimensions() + 1, location); + // Copy the raster data and append the new band values + for (int i = 0; i < raster.getWidth(); i++) { + for (int j = 0; j < raster.getHeight(); j++) { + double[] pixels = raster.getPixel(i, j, (double[]) null); + double[] copiedPixels = new double[pixels.length + 1]; + System.arraycopy(pixels, 0, copiedPixels, 0, pixels.length); + copiedPixels[pixels.length] = bandValues[j * raster.getWidth() + i]; + wr.setPixel(i, j, copiedPixels); + } + } + // Create a new GridCoverage2D with the copied image + GridCoverageFactory gridCoverageFactory = CoverageFactoryFinder.getGridCoverageFactory(null); + return gridCoverageFactory.create(gridCoverage2D.getName(), wr, gridCoverage2D.getEnvelope()); + } + + private static GridCoverage2D copyRasterAndReplaceBand(GridCoverage2D gridCoverage2D, int bandIndex, double[] bandValues) { + // Do not allow the band index to be out of bounds + if (bandIndex < 1 || bandIndex > gridCoverage2D.getNumSampleDimensions()) { + throw new IllegalArgumentException("Band index is out of bounds. Must be between 1 and " + gridCoverage2D.getNumSampleDimensions() + ")"); + } + // Get the original image and its properties + RenderedImage originalImage = gridCoverage2D.getRenderedImage(); + Raster raster = originalImage.getData(); + WritableRaster wr = raster.createCompatibleWritableRaster(); + // Copy the raster data and replace the band values + for (int i = 0; i < raster.getWidth(); i++) { + for (int j = 0; j < raster.getHeight(); j++) { + double[] bands = raster.getPixel(i, j, (double[]) null); + bands[bandIndex - 1] = bandValues[j * raster.getWidth() + i]; + wr.setPixel(i, j, bands); + } + } + // Create a new GridCoverage2D with the copied image + GridCoverageFactory gridCoverageFactory = CoverageFactoryFinder.getGridCoverageFactory(null); + return gridCoverageFactory.create(gridCoverage2D.getName(), wr, gridCoverage2D.getEnvelope()); + } +} diff --git a/common/src/main/java/org/apache/sedona/common/raster/PixelFunctions.java b/common/src/main/java/org/apache/sedona/common/raster/PixelFunctions.java new file mode 100644 index 00000000..6d6ca6c3 --- /dev/null +++ b/common/src/main/java/org/apache/sedona/common/raster/PixelFunctions.java @@ -0,0 +1,86 @@ +/* + * 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 org.geotools.coverage.grid.GridCoordinates2D; +import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.coverage.grid.GridGeometry2D; +import org.geotools.geometry.DirectPosition2D; +import org.locationtech.jts.geom.Geometry; +import org.locationtech.jts.geom.Point; +import org.opengis.referencing.operation.TransformException; + +import java.awt.image.Raster; +import java.util.ArrayList; +import java.util.Collections; +import java.util.List; +import java.util.function.DoublePredicate; +import java.util.stream.Collectors; +import java.util.stream.DoubleStream; + +public class PixelFunctions +{ + public static Double value(GridCoverage2D rasterGeom, Geometry geometry, int band) throws TransformException + { + return values(rasterGeom, Collections.singletonList(geometry), band).get(0); + } + + public static List<Double> values(GridCoverage2D rasterGeom, List<Geometry> geometries, int band) throws TransformException { + int numBands = rasterGeom.getNumSampleDimensions(); + if (band < 1 || band > numBands) { + // Invalid band index. Return nulls. + return geometries.stream().map(geom -> (Double) null).collect(Collectors.toList()); + } + Raster raster = rasterGeom.getRenderedImage().getData(); + GridGeometry2D gridGeometry = rasterGeom.getGridGeometry(); + double[] noDataValues = rasterGeom.getSampleDimension(band - 1).getNoDataValues(); + DoublePredicate isNoData = d -> noDataValues != null && DoubleStream.of(noDataValues).anyMatch(noDataValue -> Double.compare(noDataValue, d) == 0); + double[] pixelBuffer = new double[numBands]; + + List<Double> result = new ArrayList<>(geometries.size()); + for (Geometry geom : geometries) { + if (geom == null) { + result.add(null); + } else { + Point point = ensurePoint(geom); + DirectPosition2D directPosition2D = new DirectPosition2D(point.getX(), point.getY()); + GridCoordinates2D gridCoordinates2D = gridGeometry.worldToGrid(directPosition2D); + try { + double pixel = raster.getPixel(gridCoordinates2D.x, gridCoordinates2D.y, pixelBuffer)[band - 1]; + if (isNoData.test(pixel)) { + result.add(null); + } else { + result.add(pixel); + } + } catch (ArrayIndexOutOfBoundsException exc) { + // Points outside the extent should return null + result.add(null); + } + } + } + return result; + } + + private static Point ensurePoint(Geometry geometry) { + if (geometry instanceof Point) { + return (Point) geometry; + } + throw new IllegalArgumentException("Attempting to get the value of a pixel with a non-point geometry."); + } +} diff --git a/common/src/main/java/org/apache/sedona/common/raster/RasterAccessors.java b/common/src/main/java/org/apache/sedona/common/raster/RasterAccessors.java new file mode 100644 index 00000000..9f37aea6 --- /dev/null +++ b/common/src/main/java/org/apache/sedona/common/raster/RasterAccessors.java @@ -0,0 +1,102 @@ +/* + * 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 org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.geometry.Envelope2D; +import org.geotools.referencing.CRS; +import org.geotools.referencing.crs.DefaultEngineeringCRS; +import org.geotools.referencing.operation.transform.AffineTransform2D; +import org.locationtech.jts.geom.Envelope; +import org.locationtech.jts.geom.Geometry; +import org.locationtech.jts.geom.GeometryFactory; +import org.locationtech.jts.geom.PrecisionModel; +import org.opengis.referencing.FactoryException; +import org.opengis.referencing.crs.CoordinateReferenceSystem; +import org.opengis.referencing.operation.MathTransform; + +import java.awt.image.RenderedImage; +import java.util.Optional; + +public class RasterAccessors +{ + public static int srid(GridCoverage2D raster) throws FactoryException + { + CoordinateReferenceSystem crs = raster.getCoordinateReferenceSystem(); + if (crs instanceof DefaultEngineeringCRS) { + // GeoTools defaults to internal non-standard epsg codes, like 404000, if crs is missing. + // We need to check for this case and return 0 instead. + if (((DefaultEngineeringCRS) crs).isWildcard()) { + return 0; + } + } + return Optional.ofNullable(CRS.lookupEpsgCode(crs, true)).orElse(0); + } + + public static int numBands(GridCoverage2D raster) { + return raster.getNumSampleDimensions(); + } + + public static Geometry envelope(GridCoverage2D raster) throws FactoryException { + Envelope2D envelope2D = raster.getEnvelope2D(); + + Envelope envelope = new Envelope(envelope2D.getMinX(), envelope2D.getMaxX(), envelope2D.getMinY(), envelope2D.getMaxY()); + int srid = srid(raster); + return new GeometryFactory(new PrecisionModel(), srid).toGeometry(envelope); + } + + /** + * Returns the metadata of a raster as an array of doubles. + * @param raster + * @return double[] with the following values: + * 0: minX: upper left x + * 1: maxY: upper left y + * 2: width: number of pixels on x axis + * 3: height: number of pixels on y axis + * 4: scaleX: pixel width + * 5: scaleY: pixel height + * 6: shearX: skew on x axis + * 7: shearY: skew on y axis + * 8: srid + * 9: numBands + * @throws FactoryException + */ + public static double[] metadata(GridCoverage2D raster) + throws FactoryException + { + RenderedImage image = raster.getRenderedImage(); + // Georeference metadata + Envelope2D envelope2D = raster.getEnvelope2D(); + MathTransform gridToCRS = raster.getGridGeometry().getGridToCRS2D(); + if (gridToCRS instanceof AffineTransform2D) { + AffineTransform2D affine = (AffineTransform2D) gridToCRS; + + // Get the affine parameters + double scaleX = affine.getScaleX(); + double scaleY = affine.getScaleY(); + double shearX = affine.getShearX(); + double shearY = affine.getShearY(); + return new double[] {envelope2D.getMinX(), envelope2D.getMaxY(), image.getWidth(), image.getHeight(), scaleX, scaleY, shearX, shearY, srid(raster), raster.getNumSampleDimensions()}; + } + else { + // Handle the case where gridToCRS is not an AffineTransform2D + throw new UnsupportedOperationException("Only AffineTransform2D is supported"); + } + } +} diff --git a/common/src/main/java/org/apache/sedona/common/raster/RasterConstructors.java b/common/src/main/java/org/apache/sedona/common/raster/RasterConstructors.java new file mode 100644 index 00000000..8455c98e --- /dev/null +++ b/common/src/main/java/org/apache/sedona/common/raster/RasterConstructors.java @@ -0,0 +1,94 @@ +/** + * Licensed 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 + * <p> + * http://www.apache.org/licenses/LICENSE-2.0 + * <p> + * 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 org.geotools.coverage.CoverageFactoryFinder; +import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.coverage.grid.GridCoverageFactory; +import org.geotools.gce.arcgrid.ArcGridReader; +import org.geotools.gce.geotiff.GeoTiffReader; +import org.geotools.geometry.jts.ReferencedEnvelope; +import org.geotools.referencing.crs.DefaultEngineeringCRS; + +import javax.media.jai.RasterFactory; + +import java.awt.image.DataBuffer; +import java.awt.image.WritableRaster; +import java.io.ByteArrayInputStream; +import java.io.IOException; + +public class RasterConstructors +{ + + public static GridCoverage2D fromArcInfoAsciiGrid(byte[] bytes) throws IOException { + try (ByteArrayInputStream inputStream = new ByteArrayInputStream(bytes)) { + ArcGridReader reader = new ArcGridReader(inputStream); + return reader.read(null); + } + } + + public static GridCoverage2D fromGeoTiff(byte[] bytes) throws IOException { + try (ByteArrayInputStream inputStream = new ByteArrayInputStream(bytes)) { + GeoTiffReader geoTiffReader = new GeoTiffReader(inputStream); + return geoTiffReader.read(null); + } + } + + /** + * Create a new empty raster with the given number of empty bands + * The bounding envelope is defined by the upper left corner and the scale + * The math formula of the envelope is: minX = upperLeftX = lowerLeftX (minX), minY (lowerLeftY) = upperLeftY - height * pixelSizeInDegrees + * The raster is defined by the width and height + * The affine transform is defined by the skewX and skewY + * The upper left corner is defined by the upperLeftX and upperLeftY + * The scale is defined by the scaleX and scaleY + * @param widthInPixel + * @param heightInPixel + * @param upperLeftX the upper left corner of the raster. Note that: the minX of the envelope is equal to the upperLeftX + * @param upperLeftY the upper left corner of the raster. Note that: the minY of the envelope is equal to the upperLeftY - height * scaleY + * @param pixelSizeInDegrees the size of the pixel in degrees + * @param numBand the number of bands + * @return + */ + public static GridCoverage2D makeEmptyRaster(int widthInPixel, int heightInPixel, double upperLeftX, double upperLeftY, double pixelSizeInDegrees, int numBand) + { + // Create a new empty raster + WritableRaster raster = RasterFactory.createBandedRaster(DataBuffer.TYPE_DOUBLE, widthInPixel, heightInPixel, numBand, null); + + ReferencedEnvelope referencedEnvelope = new ReferencedEnvelope(upperLeftX, upperLeftX + widthInPixel * pixelSizeInDegrees, upperLeftY - heightInPixel * pixelSizeInDegrees, upperLeftY, DefaultEngineeringCRS.GENERIC_2D); + // Create a new coverage + GridCoverageFactory gridCoverageFactory = CoverageFactoryFinder.getGridCoverageFactory(null); + return gridCoverageFactory.create("genericCoverage", raster, referencedEnvelope); + } + + /** + * Create a new empty raster with 1 empty band + * The bounding envelope is defined by the upper left corner and the scale + * The math formula of the envelope is: minX = upperLeftX = lowerLeftX (minX), minY (lowerLeftY) = upperLeftY - height * pixelSizeInDegrees + * The raster is defined by the width and height + * The affine transform is defined by the skewX and skewY + * The upper left corner is defined by the upperLeftX and upperLeftY + * The scale is defined by the scaleX and scaleY + * @param widthInPixel + * @param heightInPixel + * @param upperLeftX the upper left corner of the raster. Note that: the minX of the envelope is equal to the upperLeftX + * @param upperLeftY the upper left corner of the raster. Note that: the minY of the envelope is equal to the upperLeftY - height * scaleY + * @param pixelSizeInDegrees the size of the pixel in degrees + * @return + */ + public static GridCoverage2D makeEmptyRaster(int widthInPixel, int heightInPixel, double upperLeftX, double upperLeftY, double pixelSizeInDegrees) + { + return makeEmptyRaster(widthInPixel, heightInPixel, upperLeftX, upperLeftY, pixelSizeInDegrees, 1); + } +} diff --git a/common/src/main/java/org/apache/sedona/common/raster/RasterEditors.java b/common/src/main/java/org/apache/sedona/common/raster/RasterEditors.java new file mode 100644 index 00000000..32cdc3f5 --- /dev/null +++ b/common/src/main/java/org/apache/sedona/common/raster/RasterEditors.java @@ -0,0 +1,44 @@ +/* + * 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 org.geotools.coverage.CoverageFactoryFinder; +import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.coverage.grid.GridCoverageFactory; +import org.geotools.geometry.jts.ReferencedEnvelope; +import org.geotools.referencing.CRS; +import org.geotools.referencing.crs.DefaultEngineeringCRS; +import org.opengis.referencing.FactoryException; +import org.opengis.referencing.crs.CoordinateReferenceSystem; + +public class RasterEditors +{ + public static GridCoverage2D setSrid(GridCoverage2D raster, int srid) throws FactoryException + { + CoordinateReferenceSystem crs; + if (srid == 0) { + crs = DefaultEngineeringCRS.CARTESIAN_2D; + } else { + crs = CRS.decode("EPSG:" + srid); + } + ReferencedEnvelope referencedEnvelope = new ReferencedEnvelope(raster.getEnvelope2D(), crs); + GridCoverageFactory gridCoverageFactory = CoverageFactoryFinder.getGridCoverageFactory(null); + return gridCoverageFactory.create(raster.getName().toString(), raster.getRenderedImage(), referencedEnvelope); + } +} diff --git a/common/src/main/java/org/apache/sedona/common/raster/Outputs.java b/common/src/main/java/org/apache/sedona/common/raster/RasterOutputs.java similarity index 99% rename from common/src/main/java/org/apache/sedona/common/raster/Outputs.java rename to common/src/main/java/org/apache/sedona/common/raster/RasterOutputs.java index 9d6f5bed..8fbf8217 100644 --- a/common/src/main/java/org/apache/sedona/common/raster/Outputs.java +++ b/common/src/main/java/org/apache/sedona/common/raster/RasterOutputs.java @@ -33,7 +33,7 @@ import javax.imageio.ImageWriteParam; import java.io.ByteArrayOutputStream; import java.io.IOException; -public class Outputs +public class RasterOutputs { public static byte[] asGeoTiff(GridCoverage2D raster, String compressionType, float compressionQuality) { ByteArrayOutputStream out = new ByteArrayOutputStream(); diff --git a/common/src/main/java/org/apache/sedona/common/raster/Serde.java b/common/src/main/java/org/apache/sedona/common/raster/Serde.java index 1762a543..2a9e4d7e 100644 --- a/common/src/main/java/org/apache/sedona/common/raster/Serde.java +++ b/common/src/main/java/org/apache/sedona/common/raster/Serde.java @@ -17,8 +17,13 @@ import org.geotools.coverage.grid.GridCoverage2D; import org.geotools.coverage.grid.GridCoverageFactory; import javax.media.jai.remote.SerializableRenderedImage; + import java.awt.image.RenderedImage; -import java.io.*; +import java.io.ByteArrayInputStream; +import java.io.ByteArrayOutputStream; +import java.io.IOException; +import java.io.ObjectInputStream; +import java.io.ObjectOutputStream; public class Serde { diff --git a/common/src/test/java/org/apache/sedona/common/raster/FunctionsTest.java b/common/src/test/java/org/apache/sedona/common/raster/FunctionsTest.java index 90e0d68f..287d6456 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/FunctionsTest.java +++ b/common/src/test/java/org/apache/sedona/common/raster/FunctionsTest.java @@ -26,64 +26,45 @@ import java.util.Arrays; import java.util.List; import java.util.Objects; -import static org.junit.Assert.*; +import static org.junit.Assert.assertEquals; +import static org.junit.Assert.assertNotNull; +import static org.junit.Assert.assertNull; +import static org.junit.Assert.assertTrue; public class FunctionsTest extends RasterTestBase { - @Test - public void envelope() throws FactoryException { - Geometry envelope = Functions.envelope(oneBandRaster); - assertEquals(3600.0d, envelope.getArea(), 0.1d); - assertEquals(378922.0d + 30.0d, envelope.getCentroid().getX(), 0.1d); - assertEquals(4072345.0d + 30.0d, envelope.getCentroid().getY(), 0.1d); - - assertEquals(4326, Functions.envelope(multiBandRaster).getSRID()); - } - - @Test - public void testNumBands() { - assertEquals(1, Functions.numBands(oneBandRaster)); - assertEquals(4, Functions.numBands(multiBandRaster)); - } - @Test public void testSetSrid() throws FactoryException { - assertEquals(0, Functions.srid(oneBandRaster)); - assertEquals(4326, Functions.srid(multiBandRaster)); + assertEquals(0, RasterAccessors.srid(oneBandRaster)); + assertEquals(4326, RasterAccessors.srid(multiBandRaster)); - GridCoverage2D oneBandRasterWithUpdatedSrid = Functions.setSrid(oneBandRaster, 4326); - assertEquals(4326, Functions.srid(oneBandRasterWithUpdatedSrid)); - assertEquals(4326, Functions.envelope(oneBandRasterWithUpdatedSrid).getSRID()); - assertTrue(Functions.envelope(oneBandRasterWithUpdatedSrid).equalsTopo(Functions.envelope(oneBandRaster))); + GridCoverage2D oneBandRasterWithUpdatedSrid = RasterEditors.setSrid(oneBandRaster, 4326); + assertEquals(4326, RasterAccessors.srid(oneBandRasterWithUpdatedSrid)); + assertEquals(4326, RasterAccessors.envelope(oneBandRasterWithUpdatedSrid).getSRID()); + assertTrue(RasterAccessors.envelope(oneBandRasterWithUpdatedSrid).equalsTopo(RasterAccessors.envelope(oneBandRaster))); - GridCoverage2D multiBandRasterWithUpdatedSrid = Functions.setSrid(multiBandRaster, 0); - assertEquals(0 , Functions.srid(multiBandRasterWithUpdatedSrid)); - } - - @Test - public void testSrid() throws FactoryException { - assertEquals(0, Functions.srid(oneBandRaster)); - assertEquals(4326, Functions.srid(multiBandRaster)); + GridCoverage2D multiBandRasterWithUpdatedSrid = RasterEditors.setSrid(multiBandRaster, 0); + assertEquals(0 , RasterAccessors.srid(multiBandRasterWithUpdatedSrid)); } @Test public void value() throws TransformException { - assertNull("Points outside of the envelope should return null.", Functions.value(oneBandRaster, point(1, 1), 1)); - assertNull("Invalid band should return null.", Functions.value(oneBandRaster, point(378923, 4072346), 0)); - assertNull("Invalid band should return null.", Functions.value(oneBandRaster, point(378923, 4072346), 2)); + assertNull("Points outside of the envelope should return null.", PixelFunctions.value(oneBandRaster, point(1, 1), 1)); + assertNull("Invalid band should return null.", PixelFunctions.value(oneBandRaster, point(378923, 4072346), 0)); + assertNull("Invalid band should return null.", PixelFunctions.value(oneBandRaster, point(378923, 4072346), 2)); - Double value = Functions.value(oneBandRaster, point(378923, 4072346), 1); + Double value = PixelFunctions.value(oneBandRaster, point(378923, 4072346), 1); assertNotNull(value); assertEquals(2.0d, value, 0.1d); - assertNull("Null should be returned for no data values.", Functions.value(oneBandRaster, point(378923, 4072376), 1)); + assertNull("Null should be returned for no data values.", PixelFunctions.value(oneBandRaster, point(378923, 4072376), 1)); } @Test public void valueWithMultibandRaster() throws TransformException { // Multiband raster - assertEquals(9d, Functions.value(multiBandRaster, point(4.5d,4.5d), 3), 0.1d); - assertEquals(255d, Functions.value(multiBandRaster, point(4.5d,4.5d), 4), 0.1d); + assertEquals(9d, PixelFunctions.value(multiBandRaster, point(4.5d,4.5d), 3), 0.1d); + assertEquals(255d, PixelFunctions.value(multiBandRaster, point(4.5d,4.5d), 4), 0.1d); } @Test @@ -91,19 +72,19 @@ public class FunctionsTest extends RasterTestBase { // The function 'value' is implemented using 'values'. // These test only cover bits not already covered by tests for 'value' List<Geometry> points = Arrays.asList(new Geometry[]{point(378923, 4072346), point(378924, 4072346)}); - List<Double> values = Functions.values(oneBandRaster, points, 1); + List<Double> values = PixelFunctions.values(oneBandRaster, points, 1); assertEquals(2, values.size()); assertTrue(values.stream().allMatch(Objects::nonNull)); - values = Functions.values(oneBandRaster, points, 0); + values = PixelFunctions.values(oneBandRaster, points, 0); assertEquals(2, values.size()); assertTrue("All values should be null for invalid band index.", values.stream().allMatch(Objects::isNull)); - values = Functions.values(oneBandRaster, points, 2); + values = PixelFunctions.values(oneBandRaster, points, 2); assertEquals(2, values.size()); assertTrue("All values should be null for invalid band index.", values.stream().allMatch(Objects::isNull)); - values = Functions.values(oneBandRaster, Arrays.asList(new Geometry[]{point(378923, 4072346), null}), 1); + values = PixelFunctions.values(oneBandRaster, Arrays.asList(new Geometry[]{point(378923, 4072346), null}), 1); assertEquals(2, values.size()); assertNull("Null geometries should return null values.", values.get(1)); } diff --git a/common/src/test/java/org/apache/sedona/common/raster/MapAlgebraTest.java b/common/src/test/java/org/apache/sedona/common/raster/MapAlgebraTest.java new file mode 100644 index 00000000..085f1439 --- /dev/null +++ b/common/src/test/java/org/apache/sedona/common/raster/MapAlgebraTest.java @@ -0,0 +1,147 @@ +/* + * 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 org.geotools.coverage.grid.GridCoverage2D; +import org.junit.Test; +import org.opengis.referencing.FactoryException; + +import java.io.IOException; + +import static org.junit.Assert.assertEquals; +import static org.junit.Assert.assertNull; + +public class MapAlgebraTest extends RasterTestBase +{ + @Test + public void testAddBandAsArrayAppend() + throws FactoryException + { + GridCoverage2D raster = createEmptyRaster(1); + double[] band1 = new double[raster.getRenderedImage().getWidth() * raster.getRenderedImage().getHeight()]; + for (int i = 0; i < band1.length; i++) { + band1[i] = i; + } + double[] band2 = new double[raster.getRenderedImage().getWidth() * raster.getRenderedImage().getHeight()]; + for (int i = 0; i < band2.length; i++) { + band2[i] = i * 2; + } + // Replace the first band + GridCoverage2D rasterWithBand1 = MapAlgebra.addBandFromArray(raster, band1, 1); + assertEquals(1, RasterAccessors.numBands(rasterWithBand1)); + assertEquals(raster.getEnvelope(), rasterWithBand1.getEnvelope()); + assertEquals(raster.getCoordinateReferenceSystem2D(), rasterWithBand1.getCoordinateReferenceSystem2D()); + assertEquals(RasterAccessors.srid(raster), RasterAccessors.srid(rasterWithBand1)); + + // Append a new band + GridCoverage2D rasterWithBand2 = MapAlgebra.addBandFromArray(rasterWithBand1, band2); + assertEquals(2, RasterAccessors.numBands(rasterWithBand2)); + assertEquals(raster.getEnvelope(), rasterWithBand2.getEnvelope()); + assertEquals(raster.getCoordinateReferenceSystem2D(), rasterWithBand2.getCoordinateReferenceSystem2D()); + assertEquals(RasterAccessors.srid(raster), RasterAccessors.srid(rasterWithBand2)); + + // Check the value of the first band when use the raster with only one band + double[] firstBand = MapAlgebra.bandAsArray(rasterWithBand1, 1); + for (int i = 0; i < firstBand.length; i++) { + assertEquals(i, firstBand[i], 0.1); + } + // Check the value of the first band when use the raster with two bands + firstBand = MapAlgebra.bandAsArray(rasterWithBand2, 1); + for (int i = 0; i < firstBand.length; i++) { + assertEquals(i, firstBand[i], 0.1); + } + // Check the value of the second band + double[] secondBand = MapAlgebra.bandAsArray(rasterWithBand2, 2); + for (int i = 0; i < secondBand.length; i++) { + assertEquals(i * 2, secondBand[i], 0.1); + } + } + + @Test + public void testAddBandAsArrayReplace() + throws FactoryException + { + GridCoverage2D raster = createEmptyRaster(2); + double[] band1 = new double[raster.getRenderedImage().getWidth() * raster.getRenderedImage().getHeight()]; + for (int i = 0; i < band1.length; i++) { + band1[i] = i; + } + double[] band2 = new double[raster.getRenderedImage().getWidth() * raster.getRenderedImage().getHeight()]; + for (int i = 0; i < band2.length; i++) { + band2[i] = i * 2; + } + // Replace the first band + GridCoverage2D rasterWithBand1 = MapAlgebra.addBandFromArray(raster, band1, 1); + assertEquals(2, RasterAccessors.numBands(rasterWithBand1)); + assertEquals(raster.getEnvelope(), rasterWithBand1.getEnvelope()); + assertEquals(raster.getCoordinateReferenceSystem2D(), rasterWithBand1.getCoordinateReferenceSystem2D()); + assertEquals(RasterAccessors.srid(raster), RasterAccessors.srid(rasterWithBand1)); + + // Replace the second band + GridCoverage2D rasterWithBand2 = MapAlgebra.addBandFromArray(rasterWithBand1, band2, 2); + assertEquals(2, RasterAccessors.numBands(rasterWithBand2)); + assertEquals(raster.getEnvelope(), rasterWithBand2.getEnvelope()); + assertEquals(raster.getCoordinateReferenceSystem2D(), rasterWithBand2.getCoordinateReferenceSystem2D()); + assertEquals(RasterAccessors.srid(raster), RasterAccessors.srid(rasterWithBand2)); + + // Check the value of the first band when use the raster with only one band + double[] firstBand = MapAlgebra.bandAsArray(rasterWithBand1, 1); + for (int i = 0; i < firstBand.length; i++) { + assertEquals(i, firstBand[i], 0.1); + } + // Check the value of the first band when use the raster with two bands + firstBand = MapAlgebra.bandAsArray(rasterWithBand2, 1); + for (int i = 0; i < firstBand.length; i++) { + assertEquals(i, firstBand[i], 0.1); + } + // Check the value of the second band + double[] secondBand = MapAlgebra.bandAsArray(rasterWithBand2, 2); + for (int i = 0; i < secondBand.length; i++) { + assertEquals(i * 2, secondBand[i], 0.1); + } + } + + @Test + public void testBandAsArray() { + int widthInPixel = 10; + int heightInPixel = 10; + double upperLeftX = 0; + double upperLeftY = 0; + double cellSize = 1; + GridCoverage2D raster = RasterConstructors.makeEmptyRaster(widthInPixel, heightInPixel, upperLeftX, upperLeftY, cellSize); + // Out of bound index should return null + double[] band = MapAlgebra.bandAsArray(raster, 0); + assertNull(band); + band = MapAlgebra.bandAsArray(raster, 1); + assertEquals(widthInPixel * heightInPixel, band.length); + for (int i = 0; i < band.length; i++) { + // The default value is 0.0 + assertEquals(0.0, band[i], 0.1); + } + // Now set the value of the first band and check again + for (int i = 0; i < band.length; i++) { + band[i] = i; + } + double[] bandNew = MapAlgebra.bandAsArray(MapAlgebra.addBandFromArray(raster, band, 1), 1); + assertEquals(band.length, bandNew.length); + for (int i = 0; i < band.length; i++) { + assertEquals(band[i], bandNew[i], 0.1); + } + } +} diff --git a/common/src/test/java/org/apache/sedona/common/raster/RasterAccessorsTest.java b/common/src/test/java/org/apache/sedona/common/raster/RasterAccessorsTest.java new file mode 100644 index 00000000..91ffff80 --- /dev/null +++ b/common/src/test/java/org/apache/sedona/common/raster/RasterAccessorsTest.java @@ -0,0 +1,101 @@ +/* + * 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 org.geotools.coverage.grid.GridCoverage2D; +import org.junit.Test; +import org.locationtech.jts.geom.Geometry; +import org.opengis.referencing.FactoryException; + +import static org.junit.Assert.assertEquals; + +public class RasterAccessorsTest extends RasterTestBase +{ + @Test + public void envelope() throws FactoryException + { + Geometry envelope = RasterAccessors.envelope(oneBandRaster); + assertEquals(3600.0d, envelope.getArea(), 0.1d); + assertEquals(378922.0d + 30.0d, envelope.getCentroid().getX(), 0.1d); + assertEquals(4072345.0d + 30.0d, envelope.getCentroid().getY(), 0.1d); + + assertEquals(4326, RasterAccessors.envelope(multiBandRaster).getSRID()); + } + + @Test + public void testNumBands() { + assertEquals(1, RasterAccessors.numBands(oneBandRaster)); + assertEquals(4, RasterAccessors.numBands(multiBandRaster)); + } + + @Test + public void testSrid() throws FactoryException { + assertEquals(0, RasterAccessors.srid(oneBandRaster)); + assertEquals(4326, RasterAccessors.srid(multiBandRaster)); + } + + @Test + public void testMetaData() + throws FactoryException + { + double upperLeftX = 1; + double upperLeftY = 2; + int widthInPixel = 3; + int heightInPixel = 4; + double pixelSize = 5; + + GridCoverage2D gridCoverage2D = RasterConstructors.makeEmptyRaster(widthInPixel, heightInPixel, upperLeftX, upperLeftY, pixelSize); + double[] metadata = RasterAccessors.metadata(gridCoverage2D); + assertEquals(upperLeftX, metadata[0], 0.1d); + assertEquals(upperLeftY, metadata[1], 0.1d); + assertEquals(widthInPixel, metadata[2], 0.1d); + assertEquals(heightInPixel, metadata[3], 0.1d); + assertEquals(pixelSize, metadata[4], 0.1d); + assertEquals(-1 * pixelSize, metadata[5], 0.1d); + assertEquals(0, metadata[6], 0.1d); + assertEquals(0, metadata[7], 0.1d); + assertEquals(0, metadata[8], 0.1d); + assertEquals(1, metadata[9], 0.1d); + assertEquals(10, metadata.length); + + upperLeftX = 5; + upperLeftY = 6; + widthInPixel = 7; + heightInPixel = 8; + pixelSize = 9; + int numBands = 10; + + gridCoverage2D = RasterConstructors.makeEmptyRaster(widthInPixel, heightInPixel, upperLeftX, upperLeftY, pixelSize, numBands); + + metadata = RasterAccessors.metadata(gridCoverage2D); + + assertEquals(upperLeftX, metadata[0], 0.1d); + assertEquals(upperLeftY, metadata[1], 0.1d); + assertEquals(widthInPixel, metadata[2], 0.1d); + assertEquals(heightInPixel, metadata[3], 0.1d); + assertEquals(pixelSize, metadata[4], 0.1d); + assertEquals(-1 * pixelSize, metadata[5], 0.1d); + assertEquals(0, metadata[6], 0.1d); + assertEquals(0, metadata[7], 0.1d); + assertEquals(0, metadata[8], 0.1d); + assertEquals(numBands, metadata[9], 0.1d); + + assertEquals(10, metadata.length); + } +} diff --git a/common/src/test/java/org/apache/sedona/common/raster/ConstructorsTest.java b/common/src/test/java/org/apache/sedona/common/raster/RasterConstructorsTest.java similarity index 51% rename from common/src/test/java/org/apache/sedona/common/raster/ConstructorsTest.java rename to common/src/test/java/org/apache/sedona/common/raster/RasterConstructorsTest.java index f2ec505b..7f511844 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/ConstructorsTest.java +++ b/common/src/test/java/org/apache/sedona/common/raster/RasterConstructorsTest.java @@ -21,15 +21,16 @@ import org.opengis.referencing.FactoryException; import java.io.IOException; import java.nio.charset.StandardCharsets; -import static org.junit.Assert.*; +import static org.junit.Assert.assertEquals; -public class ConstructorsTest extends RasterTestBase { +public class RasterConstructorsTest + extends RasterTestBase { @Test public void fromArcInfoAsciiGrid() throws IOException, FactoryException { - GridCoverage2D gridCoverage2D = Constructors.fromArcInfoAsciiGrid(arc.getBytes(StandardCharsets.UTF_8)); + GridCoverage2D gridCoverage2D = RasterConstructors.fromArcInfoAsciiGrid(arc.getBytes(StandardCharsets.UTF_8)); - Geometry envelope = Functions.envelope(gridCoverage2D); + Geometry envelope = RasterAccessors.envelope(gridCoverage2D); assertEquals(3600, envelope.getArea(), 0.1); assertEquals(378922d + 30, envelope.getCentroid().getX(), 0.1); assertEquals(4072345d + 30, envelope.getCentroid().getY(), 0.1); @@ -41,9 +42,9 @@ public class ConstructorsTest extends RasterTestBase { @Test public void fromGeoTiff() throws IOException, FactoryException { - GridCoverage2D gridCoverage2D = Constructors.fromGeoTiff(geoTiff); + GridCoverage2D gridCoverage2D = RasterConstructors.fromGeoTiff(geoTiff); - Geometry envelope = Functions.envelope(gridCoverage2D); + Geometry envelope = RasterAccessors.envelope(gridCoverage2D); assertEquals(100, envelope.getArea(), 0.1); assertEquals(5, envelope.getCentroid().getX(), 0.1); assertEquals(5, envelope.getCentroid().getY(), 0.1); @@ -52,4 +53,30 @@ public class ConstructorsTest extends RasterTestBase { assertEquals(10d, gridCoverage2D.getRenderedImage().getData().getPixel(5, 5, (double[])null)[0], 0.1); assertEquals(4, gridCoverage2D.getNumSampleDimensions()); } + + @Test + public void makeEmptyRaster() throws FactoryException { + double upperLeftX = 0; + double upperLeftY = 0; + int widthInPixel = 1; + int heightInPixel = 2; + double pixelSize = 2; + + GridCoverage2D gridCoverage2D = RasterConstructors.makeEmptyRaster(widthInPixel, heightInPixel, upperLeftX, upperLeftY, pixelSize); + Geometry envelope = RasterAccessors.envelope(gridCoverage2D); + assertEquals(upperLeftX, envelope.getEnvelopeInternal().getMinX(), 0.001); + assertEquals(upperLeftX + widthInPixel * pixelSize, envelope.getEnvelopeInternal().getMaxX(), 0.001); + assertEquals(upperLeftY - heightInPixel * pixelSize, envelope.getEnvelopeInternal().getMinY(), 0.001); + assertEquals(upperLeftY, envelope.getEnvelopeInternal().getMaxY(), 0.001); + + assertEquals("POLYGON ((0 -4, 0 0, 2 0, 2 -4, 0 -4))", envelope.toString()); + double expectedWidthInDegree = pixelSize * widthInPixel; + double expectedHeightInDegree = pixelSize * heightInPixel; + + assertEquals(expectedWidthInDegree * expectedHeightInDegree, envelope.getArea(), 0.001); + assertEquals(heightInPixel, gridCoverage2D.getRenderedImage().getTileHeight()); + assertEquals(widthInPixel, gridCoverage2D.getRenderedImage().getTileWidth()); + assertEquals(0d, gridCoverage2D.getRenderedImage().getData().getPixel(0, 0, (double[])null)[0], 0.001); + assertEquals(1, gridCoverage2D.getNumSampleDimensions()); + } } \ No newline at end of file diff --git a/common/src/test/java/org/apache/sedona/common/raster/RasterTestBase.java b/common/src/test/java/org/apache/sedona/common/raster/RasterTestBase.java index ba19e1f1..69d08057 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/RasterTestBase.java +++ b/common/src/test/java/org/apache/sedona/common/raster/RasterTestBase.java @@ -21,7 +21,7 @@ import org.geotools.referencing.crs.DefaultGeographicCRS; import org.junit.Before; import org.opengis.parameter.GeneralParameterValue; -import java.awt.*; +import java.awt.Color; import java.awt.image.BufferedImage; import java.io.ByteArrayOutputStream; import java.io.IOException; @@ -35,13 +35,21 @@ public class RasterTestBase { @Before public void setup() throws IOException { - oneBandRaster = Constructors.fromArcInfoAsciiGrid(arc.getBytes(StandardCharsets.UTF_8)); + oneBandRaster = RasterConstructors.fromArcInfoAsciiGrid(arc.getBytes(StandardCharsets.UTF_8)); multiBandRaster = createMultibandRaster(); ByteArrayOutputStream bos = new ByteArrayOutputStream(); new GeoTiffWriter(bos).write(multiBandRaster, new GeneralParameterValue[]{}); geoTiff = bos.toByteArray(); } + GridCoverage2D createEmptyRaster(int numBands) { + int widthInPixel = 4; + int heightInPixel = 5; + double upperLeftX = 101; + double upperLeftY = 102; + double cellSize = 2; + return RasterConstructors.makeEmptyRaster(widthInPixel, heightInPixel, upperLeftX, upperLeftY, cellSize, numBands); + } GridCoverage2D createMultibandRaster() throws IOException { GridCoverageFactory factory = new GridCoverageFactory(); BufferedImage image = new BufferedImage(10, 10, BufferedImage.TYPE_INT_ARGB); diff --git a/common/src/test/java/org/apache/sedona/common/raster/SerdeTest.java b/common/src/test/java/org/apache/sedona/common/raster/SerdeTest.java index 4500e11e..1a13534b 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/SerdeTest.java +++ b/common/src/test/java/org/apache/sedona/common/raster/SerdeTest.java @@ -18,7 +18,8 @@ import org.junit.Test; import java.io.IOException; -import static org.junit.Assert.*; +import static org.junit.Assert.assertEquals; +import static org.junit.Assert.assertNotNull; public class SerdeTest extends RasterTestBase { diff --git a/docs/api/sql/Raster-loader.md b/docs/api/sql/Raster-loader.md index 4ed58b52..520b6fae 100644 --- a/docs/api/sql/Raster-loader.md +++ b/docs/api/sql/Raster-loader.md @@ -47,8 +47,57 @@ var df = spark.read.format("binaryFile").load("/some/path/*.tiff") df = df.withColumn("raster", f.expr("RS_FromGeoTiff(content)")) ``` +### RS_MakeEmptyRaster + +Introduction: Returns an empty raster geometry.Every band in the raster is initialized to `0.0`. + +Format: `RS_MakeEmptyRaster(width: Int, height: Int, upperleftX: Double, upperleftY: Double, cellSize:Double, numBands:Int)` + +Since: `v1.4.1` + +* Width: The width of the raster in pixels. +* Height: The height of the raster in pixels. +* UpperleftX: The X coordinate of the upper left corner of the raster, in terms of the CRS units. +* UpperleftY: The Y coordinate of the upper left corner of the raster, in terms of the CRS units. +* CellSize: The size of the cells in the raster, in terms of the CRS units. +* NumBands (optional): The number of bands in the raster. If not specified, the raster will have a single band. + +SQL example 1 (with 2 bands): + +```sql +SELECT RS_MakeEmptyRaster(10, 10, 0.0, 0.0, 1.0, 2) as raster +``` + +Output: +``` ++--------------------------------------------+ +|rs_makeemptyraster(10, 10, 0.0, 0.0, 1.0, 2)| ++--------------------------------------------+ +| GridCoverage2D["g...| ++--------------------------------------------+ +``` + +SQL example 2 (number of bands not specified, default to 1 band): + +```sql +SELECT RS_MakeEmptyRaster(10, 10, 0.0, 0.0, 1.0) as raster +``` + +Output: +``` ++--------------------------------------------+ +|rs_makeemptyraster(10, 10, 0.0, 0.0, 1.0, 2)| ++--------------------------------------------+ +| GridCoverage2D["g...| ++--------------------------------------------+ +``` + + ## Load GeoTiff to Array[Double] format +!!!warning + This function has been deprecated since v1.4.1. Please use `RS_FromGeoTiff` instead and `binaryFile` data source to read GeoTiff files. + The `geotiff` loader of Sedona is a Spark built-in data source. It can read a single geotiff image or a number of geotiff images into a DataFrame. Each geotiff is a row in the resulting DataFrame and stored in an array of Double type format. Since: `v1.1.0` diff --git a/docs/api/sql/Raster-operators.md b/docs/api/sql/Raster-operators.md index 0b49daf0..fcaddaa3 100644 --- a/docs/api/sql/Raster-operators.md +++ b/docs/api/sql/Raster-operators.md @@ -1,4 +1,4 @@ -## RasterUDT based operators +## Raster based operators ### RS_Envelope @@ -17,6 +17,39 @@ Output: POLYGON((0 0,20 0,20 60,0 60,0 0)) ``` +### RS_MetaData + +Introduction: Returns the metadata of the raster as an array of double. The array contains the following values: + +- 0: upper left x coordinate of the raster, in terms of CRS units (the minimum x coordinate) +- 1: upper left y coordinate of the raster, in terms of CRS units (the maximum y coordinate) +- 2: width of the raster, in terms of pixels +- 3: height of the raster, in terms of pixels +- 4: width of a pixel, in terms of CRS units (scaleX) +- 5: height of a pixel, in terms of CRS units (scaleY) +- 6: skew in x direction (rotation x) +- 7: skew in y direction (rotation y) +- 8: srid of the raster +- 9: number of bands + +Format: `RS_MetaData (raster: Raster)` + +Since: `v1.4.1` + +SQL example: +```sql +SELECT RS_MetaData(raster) FROM raster_table +``` + +Output: +``` ++-----------------------------------------------------------------------------------------------------------------------+ +|rs_metadata(raster) | ++-----------------------------------------------------------------------------------------------------------------------+ +|[-1.3095817809482181E7, 4021262.7487925636, 512.0, 517.0, 72.32861272132695, -72.32861272132695, 0.0, 0.0, 3857.0, 1.0]| ++-----------------------------------------------------------------------------------------------------------------------+ +``` + ### RS_NumBands Introduction: Returns the number of the bands in the raster. @@ -35,7 +68,7 @@ Output: 4 ``` -## RS_SetSRID +### RS_SetSRID Introduction: Sets the spatial reference system identifier (SRID) of the raster geometry. @@ -143,7 +176,67 @@ Output: ``` -## Array[Double] based operators +## Raster to Map Algebra operators + +To bridge the gap between the raster and map algebra worlds, the following operators are provided. These operators convert a raster to a map algebra object. The map algebra object can then be used with the map algebra operators described in the next section. + +### RS_BandAsArray + +Introduction: Extract a band from a raster as an array of doubles. + +Format: `RS_BandAsArray (raster: Raster, bandIndex: Int)`. + +Since: `v1.4.1` + +BandIndex is 1-based and must be between 1 and RS_NumBands(raster). It returns null if the bandIndex is out of range or the raster is null. + +SQL example: +```sql +SELECT RS_BandAsArray(raster, 1) FROM raster_table +``` + +Output: + +``` ++--------------------+ +| band| ++--------------------+ +|[0.0, 0.0, 0.0, 0...| ++--------------------+ +``` + +### RS_AddBandFromArray + +Introduction: Add a band to a raster from an array of doubles. + +Format: `RS_AddBandFromArray (raster: Raster, band: Array[Double], bandIndex:Int)`. + +Since: `v1.4.1` + +The bandIndex is 1-based and must be between 1 and RS_NumBands(raster) + 1. It throws an exception if the bandIndex is out of range or the raster is null. + +When the bandIndex is RS_NumBands(raster) + 1, it appends the band to the end of the raster. Otherwise, it replaces the existing band at the bandIndex. + +Note that: `bandIndex == RS_NumBands(raster) + 1` is an experimental feature and might not lead to the loss of raster metadata and properties such as color models. + +SQL example: + +```sql +SELECT RS_AddBandFromArray(raster, RS_MultiplyFactor(RS_BandAsArray(RS_FromGeoTiff(content), 1), 2), 1) AS raster FROM raster_table +``` + +Output: +``` ++--------------------+ +| raster| ++--------------------+ +|GridCoverage2D["g...| ++--------------------+ +``` + +## Map Algebra operators + +Map algebra operators work on a single band of a raster. Each band is represented as an array of doubles. The operators return an array of doubles. ### RS_Add @@ -168,6 +261,8 @@ Format: `RS_Append(data: Array[Double], newBand: Array[Double], nBands: Int)` Since: `v1.2.1` +Deprecated since: `v1.4.1` + Spark SQL example: ```scala diff --git a/docs/api/sql/Raster-writer.md b/docs/api/sql/Raster-writer.md index eded6509..0ddcee04 100644 --- a/docs/api/sql/Raster-writer.md +++ b/docs/api/sql/Raster-writer.md @@ -167,6 +167,9 @@ The newly created DataFrame can be written to disk again but must be under a dif ## Write Array[Double] to GeoTiff files +!!!warning + This function has been deprecated since v1.4.1. Please use `RS_AsGeoTiff` instead and `raster` data source to write GeoTiff files. + Introduction: You can write a GeoTiff dataframe as GeoTiff images using the spark `write` feature with the format `geotiff`. The geotiff raster column needs to be an array of double type data. Since: `v1.2.1` diff --git a/sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala b/sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala index d50af2ab..2ac0dd34 100644 --- a/sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala +++ b/sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala @@ -174,10 +174,14 @@ object Catalog { function[RS_Array](), function[RS_Normalize](), function[RS_Append](), + function[RS_AddBandFromArray](), + function[RS_BandAsArray](), function[RS_FromArcInfoAsciiGrid](), function[RS_FromGeoTiff](), + function[RS_MakeEmptyRaster](1), function[RS_Envelope](), function[RS_NumBands](), + function[RS_Metadata](), function[RS_SetSRID](), function[RS_SRID](), function[RS_Value](1), diff --git a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/Functions.scala b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/MapAlgebra.scala similarity index 86% rename from sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/Functions.scala rename to sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/MapAlgebra.scala index d904bf9d..e08e9bff 100644 --- a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/Functions.scala +++ b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/MapAlgebra.scala @@ -16,24 +16,19 @@ * specific language governing permissions and limitations * under the License. */ - package org.apache.spark.sql.sedona_sql.expressions.raster -import org.apache.sedona.common.geometrySerde.GeometrySerializer -import org.apache.sedona.common.raster.{Functions, Serde} +import org.apache.sedona.common.raster.{MapAlgebra, Serde} import org.apache.spark.sql.catalyst.InternalRow import org.apache.spark.sql.catalyst.expressions.codegen.CodegenFallback import org.apache.spark.sql.catalyst.expressions.{ExpectsInputTypes, Expression} import org.apache.spark.sql.catalyst.util.{ArrayData, GenericArrayData} -import org.apache.spark.sql.sedona_sql.UDT.{GeometryUDT, RasterUDT} +import org.apache.spark.sql.sedona_sql.UDT.RasterUDT +import org.apache.spark.sql.sedona_sql.expressions.raster.implicits.RasterInputExpressionEnhancer import org.apache.spark.sql.sedona_sql.expressions.{SerdeAware, UserDataGeneratator} -import org.apache.spark.sql.sedona_sql.expressions.implicits._ -import org.apache.spark.sql.sedona_sql.expressions.raster.implicits._ import org.apache.spark.sql.types._ import org.geotools.coverage.grid.GridCoverage2D - - /// Calculate Normalized Difference between two bands case class RS_NormalizedDifference(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback with UserDataGeneratator { @@ -328,7 +323,7 @@ case class RS_Count(inputExpressions: Seq[Expression]) override def nullable: Boolean = false - override def eval(inputRow: InternalRow): Any = { + override def eval(inputRow: InternalRow): Any = { val band = inputExpressions(0).eval(inputRow).asInstanceOf[ArrayData].toDoubleArray() val target = inputExpressions(1).eval(inputRow).asInstanceOf[Decimal].toDouble findCount(band, target) @@ -810,42 +805,16 @@ case class RS_Append(inputExpressions: Seq[Expression]) } } -case class RS_Envelope(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback with ExpectsInputTypes { - override def nullable: Boolean = true - - override def eval(input: InternalRow): Any = { - val raster = inputExpressions(0).toRaster(input) - if (raster == null) { - null - } else { - Functions.envelope(raster).toGenericArrayData - } - } - - override def dataType: DataType = GeometryUDT - - override def children: Seq[Expression] = inputExpressions - - protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { - copy(inputExpressions = newChildren) - } - - override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT) -} +case class RS_AddBandFromArray(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback + with ExpectsInputTypes with SerdeAware { -case class RS_NumBands(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback with ExpectsInputTypes { override def nullable: Boolean = true override def eval(input: InternalRow): Any = { - val raster = inputExpressions(0).toRaster(input) - if (raster == null) { - null - } else { - Functions.numBands(raster) - } + Option(evalWithoutSerialization(input)).map(Serde.serialize).orNull } - override def dataType: DataType = IntegerType + override def dataType: DataType = RasterUDT override def children: Seq[Expression] = inputExpressions @@ -853,111 +822,44 @@ case class RS_NumBands(inputExpressions: Seq[Expression]) extends Expression wit copy(inputExpressions = newChildren) } - override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT) -} - -case class RS_SetSRID(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback with ExpectsInputTypes with SerdeAware { - override def nullable: Boolean = true - - override def eval(input: InternalRow): Any = { - Option(evalWithoutSerialization(input)).map(Serde.serialize).orNull - } + override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT, ArrayType, IntegerType) override def evalWithoutSerialization(input: InternalRow): GridCoverage2D = { val raster = inputExpressions(0).toRaster(input) - val srid = inputExpressions(1).eval(input).asInstanceOf[Int] if (raster == null) { - null - } else { - Functions.setSrid(raster, srid) + return null } + val band = inputExpressions(1).eval(input).asInstanceOf[ArrayData].toDoubleArray() + val bandIndex = inputExpressions(2).eval(input).asInstanceOf[Int] + MapAlgebra.addBandFromArray(raster, band, bandIndex) } - - override def dataType: DataType = RasterUDT - - override def children: Seq[Expression] = inputExpressions - - protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { - copy(inputExpressions = newChildren) - } - - override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT, IntegerType) } -case class RS_SRID(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback with ExpectsInputTypes { +case class RS_BandAsArray(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback + with ExpectsInputTypes { + override def nullable: Boolean = true override def eval(input: InternalRow): Any = { val raster = inputExpressions(0).toRaster(input) if (raster == null) { - null - } else { - Functions.srid(raster) + return null } - } - - override def dataType: DataType = IntegerType - - override def children: Seq[Expression] = inputExpressions - - protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { - copy(inputExpressions = newChildren) - } - - override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT) -} - -case class RS_Value(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback with ExpectsInputTypes { - - override def nullable: Boolean = true - - override def dataType: DataType = DoubleType - - override def eval(input: InternalRow): Any = { - val raster = inputExpressions.head.toRaster(input) - val geom = inputExpressions(1).toGeometry(input) - val band = inputExpressions(2).eval(input).asInstanceOf[Int] - if (raster == null || geom == null) { - null - } else { - Functions.value(raster, geom, band) + val bandIndex = inputExpressions(1).eval(input).asInstanceOf[Int] + val band = MapAlgebra.bandAsArray(raster, bandIndex) + if (band == null) { + return null } + new GenericArrayData(band) } - override def children: Seq[Expression] = inputExpressions - - protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { - copy(inputExpressions = newChildren) - } - - override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT, GeometryUDT, IntegerType) -} - -case class RS_Values(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback with ExpectsInputTypes { - - override def nullable: Boolean = true - override def dataType: DataType = ArrayType(DoubleType) - override def eval(input: InternalRow): Any = { - val raster = inputExpressions(0).toRaster(input) - val serializedGeometries = inputExpressions(1).eval(input).asInstanceOf[ArrayData] - val band = inputExpressions(2).eval(input).asInstanceOf[Int] - if (raster == null || serializedGeometries == null) { - null - } else { - val geometries = (0 until serializedGeometries.numElements()).map { - i => Option(serializedGeometries.getBinary(i)).map(GeometrySerializer.deserialize).orNull - } - new GenericArrayData(Functions.values(raster, java.util.Arrays.asList(geometries:_*), band).toArray) - } - } - override def children: Seq[Expression] = inputExpressions protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { copy(inputExpressions = newChildren) } - override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT, ArrayType(GeometryUDT), IntegerType) + override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT, IntegerType) } diff --git a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/PixelFunctions.scala b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/PixelFunctions.scala new file mode 100644 index 00000000..d6001651 --- /dev/null +++ b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/PixelFunctions.scala @@ -0,0 +1,85 @@ +/* + * 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.spark.sql.sedona_sql.expressions.raster + +import org.apache.sedona.common.geometrySerde.GeometrySerializer +import org.apache.sedona.common.raster.PixelFunctions +import org.apache.spark.sql.catalyst.InternalRow +import org.apache.spark.sql.catalyst.expressions.{ExpectsInputTypes, Expression} +import org.apache.spark.sql.catalyst.expressions.codegen.CodegenFallback +import org.apache.spark.sql.catalyst.util.{ArrayData, GenericArrayData} +import org.apache.spark.sql.sedona_sql.UDT.{GeometryUDT, RasterUDT} +import org.apache.spark.sql.sedona_sql.expressions.implicits.InputExpressionEnhancer +import org.apache.spark.sql.sedona_sql.expressions.raster.implicits.RasterInputExpressionEnhancer +import org.apache.spark.sql.types.{AbstractDataType, ArrayType, DataType, DoubleType, IntegerType} + +case class RS_Value(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback with ExpectsInputTypes { + + override def nullable: Boolean = true + + override def dataType: DataType = DoubleType + + override def eval(input: InternalRow): Any = { + val raster = inputExpressions.head.toRaster(input) + val geom = inputExpressions(1).toGeometry(input) + val band = inputExpressions(2).eval(input).asInstanceOf[Int] + if (raster == null || geom == null) { + null + } else { + PixelFunctions.value(raster, geom, band) + } + } + + override def children: Seq[Expression] = inputExpressions + + protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { + copy(inputExpressions = newChildren) + } + + override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT, GeometryUDT, IntegerType) +} + +case class RS_Values(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback with ExpectsInputTypes { + + override def nullable: Boolean = true + + override def dataType: DataType = ArrayType(DoubleType) + + override def eval(input: InternalRow): Any = { + val raster = inputExpressions(0).toRaster(input) + val serializedGeometries = inputExpressions(1).eval(input).asInstanceOf[ArrayData] + val band = inputExpressions(2).eval(input).asInstanceOf[Int] + if (raster == null || serializedGeometries == null) { + null + } else { + val geometries = (0 until serializedGeometries.numElements()).map { + i => Option(serializedGeometries.getBinary(i)).map(GeometrySerializer.deserialize).orNull + } + new GenericArrayData(PixelFunctions.values(raster, java.util.Arrays.asList(geometries:_*), band).toArray) + } + } + + override def children: Seq[Expression] = inputExpressions + + protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { + copy(inputExpressions = newChildren) + } + + override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT, ArrayType(GeometryUDT), IntegerType) +} \ No newline at end of file diff --git a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterAccessors.scala b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterAccessors.scala new file mode 100644 index 00000000..1f19979e --- /dev/null +++ b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterAccessors.scala @@ -0,0 +1,121 @@ +/* + * 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.spark.sql.sedona_sql.expressions.raster + +import org.apache.sedona.common.raster.RasterAccessors +import org.apache.spark.sql.catalyst.InternalRow +import org.apache.spark.sql.catalyst.expressions.{ExpectsInputTypes, Expression} +import org.apache.spark.sql.catalyst.expressions.codegen.CodegenFallback +import org.apache.spark.sql.catalyst.util.GenericArrayData +import org.apache.spark.sql.sedona_sql.UDT.{GeometryUDT, RasterUDT} +import org.apache.spark.sql.sedona_sql.expressions.implicits.GeometryEnhancer +import org.apache.spark.sql.sedona_sql.expressions.raster.implicits.RasterInputExpressionEnhancer +import org.apache.spark.sql.types.{AbstractDataType, ArrayType, DataType, DoubleType, IntegerType} + +case class RS_Envelope(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback with ExpectsInputTypes { + override def nullable: Boolean = true + + override def eval(input: InternalRow): Any = { + val raster = inputExpressions(0).toRaster(input) + if (raster == null) { + null + } else { + RasterAccessors.envelope(raster).toGenericArrayData + } + } + + override def dataType: DataType = GeometryUDT + + override def children: Seq[Expression] = inputExpressions + + protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { + copy(inputExpressions = newChildren) + } + + override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT) +} + +case class RS_NumBands(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback with ExpectsInputTypes { + override def nullable: Boolean = true + + override def eval(input: InternalRow): Any = { + val raster = inputExpressions(0).toRaster(input) + if (raster == null) { + null + } else { + RasterAccessors.numBands(raster) + } + } + + override def dataType: DataType = IntegerType + + override def children: Seq[Expression] = inputExpressions + + protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { + copy(inputExpressions = newChildren) + } + + override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT) +} + +case class RS_SRID(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback with ExpectsInputTypes { + override def nullable: Boolean = true + + override def eval(input: InternalRow): Any = { + val raster = inputExpressions(0).toRaster(input) + if (raster == null) { + null + } else { + RasterAccessors.srid(raster) + } + } + + override def dataType: DataType = IntegerType + + override def children: Seq[Expression] = inputExpressions + + protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { + copy(inputExpressions = newChildren) + } + + override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT) +} + +case class RS_Metadata(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback with ExpectsInputTypes { + override def nullable: Boolean = true + + override def eval(input: InternalRow): Any = { + val raster = inputExpressions(0).toRaster(input) + if (raster == null) { + null + } else { + new GenericArrayData(RasterAccessors.metadata(raster)) + } + } + + override def dataType: DataType = ArrayType(DoubleType) + + override def children: Seq[Expression] = inputExpressions + + protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { + copy(inputExpressions = newChildren) + } + + override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT) +} \ No newline at end of file diff --git a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/Constructors.scala b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterConstructors.scala similarity index 63% copy from sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/Constructors.scala copy to sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterConstructors.scala index c71f9963..3f2c30c9 100644 --- a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/Constructors.scala +++ b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterConstructors.scala @@ -18,14 +18,14 @@ */ package org.apache.spark.sql.sedona_sql.expressions.raster -import org.apache.sedona.common.raster.{Constructors, Serde} +import org.apache.sedona.common.raster.{RasterConstructors, Serde} import org.apache.spark.sql.catalyst.InternalRow -import org.apache.spark.sql.catalyst.expressions.{ExpectsInputTypes, Expression} import org.apache.spark.sql.catalyst.expressions.codegen.CodegenFallback +import org.apache.spark.sql.catalyst.expressions.{ExpectsInputTypes, Expression} import org.apache.spark.sql.sedona_sql.UDT.RasterUDT import org.apache.spark.sql.sedona_sql.expressions.SerdeAware -import org.apache.spark.sql.types.{AbstractDataType, BinaryType, DataType} import org.apache.spark.sql.sedona_sql.expressions.raster.implicits._ +import org.apache.spark.sql.types._ import org.geotools.coverage.grid.GridCoverage2D @@ -53,7 +53,7 @@ case class RS_FromArcInfoAsciiGrid(inputExpressions: Seq[Expression]) extends Ex if (bytes == null) { null } else { - Constructors.fromArcInfoAsciiGrid(bytes) + RasterConstructors.fromArcInfoAsciiGrid(bytes) } } } @@ -82,7 +82,40 @@ case class RS_FromGeoTiff(inputExpressions: Seq[Expression]) extends Expression if (bytes == null) { null } else { - Constructors.fromGeoTiff(bytes) + RasterConstructors.fromGeoTiff(bytes) } } } + +case class RS_MakeEmptyRaster(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback + with ExpectsInputTypes with SerdeAware { + + override def nullable: Boolean = true + + override def eval(input: InternalRow): Any = { + Option(evalWithoutSerialization(input)).map(Serde.serialize).orNull + } + + override def dataType: DataType = RasterUDT + + override def children: Seq[Expression] = inputExpressions + + protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { + copy(inputExpressions = newChildren) + } + + override def inputTypes: Seq[AbstractDataType] = Seq(IntegerType, IntegerType, DecimalType, DecimalType, DecimalType, IntegerType) + + override def evalWithoutSerialization(input: InternalRow): GridCoverage2D = { + val widthInPixels = inputExpressions(0).eval(input).asInstanceOf[Int] + val heightInPixel = inputExpressions(1).eval(input).asInstanceOf[Int] + val upperLeftX = inputExpressions(2).eval(input).asInstanceOf[Decimal].toDouble + val upperLeftY = inputExpressions(3).eval(input).asInstanceOf[Decimal].toDouble + val pixelSizeInDegree = inputExpressions(4).eval(input).asInstanceOf[Decimal].toDouble + if (inputExpressions.length == 6) { + val numBands = inputExpressions(5).eval(input).asInstanceOf[Int] + RasterConstructors.makeEmptyRaster(widthInPixels, heightInPixel, upperLeftX, upperLeftY, pixelSizeInDegree, numBands) + } + else RasterConstructors.makeEmptyRaster(widthInPixels, heightInPixel, upperLeftX, upperLeftY, pixelSizeInDegree) + } +} \ No newline at end of file diff --git a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/Constructors.scala b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterEditors.scala similarity index 56% rename from sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/Constructors.scala rename to sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterEditors.scala index c71f9963..db63fb6c 100644 --- a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/Constructors.scala +++ b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterEditors.scala @@ -18,54 +18,32 @@ */ package org.apache.spark.sql.sedona_sql.expressions.raster -import org.apache.sedona.common.raster.{Constructors, Serde} +import org.apache.sedona.common.raster.{RasterEditors, Serde} import org.apache.spark.sql.catalyst.InternalRow import org.apache.spark.sql.catalyst.expressions.{ExpectsInputTypes, Expression} import org.apache.spark.sql.catalyst.expressions.codegen.CodegenFallback import org.apache.spark.sql.sedona_sql.UDT.RasterUDT import org.apache.spark.sql.sedona_sql.expressions.SerdeAware -import org.apache.spark.sql.types.{AbstractDataType, BinaryType, DataType} -import org.apache.spark.sql.sedona_sql.expressions.raster.implicits._ +import org.apache.spark.sql.sedona_sql.expressions.raster.implicits.RasterInputExpressionEnhancer +import org.apache.spark.sql.types.{AbstractDataType, DataType, IntegerType} import org.geotools.coverage.grid.GridCoverage2D - -case class RS_FromArcInfoAsciiGrid(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback - with ExpectsInputTypes with SerdeAware { - +case class RS_SetSRID(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback with ExpectsInputTypes with SerdeAware { override def nullable: Boolean = true override def eval(input: InternalRow): Any = { Option(evalWithoutSerialization(input)).map(Serde.serialize).orNull } - override def dataType: DataType = RasterUDT - - override def children: Seq[Expression] = inputExpressions - - protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { - copy(inputExpressions = newChildren) - } - - override def inputTypes: Seq[AbstractDataType] = Seq(BinaryType) - override def evalWithoutSerialization(input: InternalRow): GridCoverage2D = { - val bytes = inputExpressions(0).eval(input).asInstanceOf[Array[Byte]] - if (bytes == null) { + val raster = inputExpressions(0).toRaster(input) + val srid = inputExpressions(1).eval(input).asInstanceOf[Int] + if (raster == null) { null } else { - Constructors.fromArcInfoAsciiGrid(bytes) + RasterEditors.setSrid(raster, srid) } } -} - -case class RS_FromGeoTiff(inputExpressions: Seq[Expression]) extends Expression with CodegenFallback - with ExpectsInputTypes with SerdeAware { - - override def nullable: Boolean = true - - override def eval(input: InternalRow): Any = { - Option(evalWithoutSerialization(input)).map(Serde.serialize).orNull - } override def dataType: DataType = RasterUDT @@ -75,14 +53,5 @@ case class RS_FromGeoTiff(inputExpressions: Seq[Expression]) extends Expression copy(inputExpressions = newChildren) } - override def inputTypes: Seq[AbstractDataType] = Seq(BinaryType) - - override def evalWithoutSerialization(input: InternalRow): GridCoverage2D = { - val bytes = inputExpressions(0).eval(input).asInstanceOf[Array[Byte]] - if (bytes == null) { - null - } else { - Constructors.fromGeoTiff(bytes) - } - } -} + override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT, IntegerType) +} \ No newline at end of file diff --git a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/Outputs.scala b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterOutputs.scala similarity index 85% rename from sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/Outputs.scala rename to sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterOutputs.scala index 632c9c46..195ed0c0 100644 --- a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/Outputs.scala +++ b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterOutputs.scala @@ -18,10 +18,10 @@ */ package org.apache.spark.sql.sedona_sql.expressions.raster -import org.apache.sedona.common.raster.Outputs +import org.apache.sedona.common.raster.RasterOutputs import org.apache.spark.sql.catalyst.InternalRow +import org.apache.spark.sql.catalyst.expressions.Expression import org.apache.spark.sql.catalyst.expressions.codegen.CodegenFallback -import org.apache.spark.sql.catalyst.expressions.{Expression, Literal} import org.apache.spark.sql.sedona_sql.expressions.raster.implicits.RasterInputExpressionEnhancer import org.apache.spark.sql.types._ import org.apache.spark.unsafe.types.UTF8String @@ -36,10 +36,10 @@ case class RS_AsGeoTiff(inputExpressions: Seq[Expression]) extends Expression wi if (raster == null) return null // If there are more than one input expressions, the additional ones are used as parameters if (inputExpressions.length > 1) { - Outputs.asGeoTiff(raster, inputExpressions(1).eval(input).asInstanceOf[UTF8String].toString, inputExpressions(2).eval(input).toString.toFloat) + RasterOutputs.asGeoTiff(raster, inputExpressions(1).eval(input).asInstanceOf[UTF8String].toString, inputExpressions(2).eval(input).toString.toFloat) } else { - Outputs.asGeoTiff(raster, null, -1) + RasterOutputs.asGeoTiff(raster, null, -1) } } @@ -61,10 +61,10 @@ case class RS_AsArcGrid(inputExpressions: Seq[Expression]) extends Expression wi if (raster == null) return null // If there are more than one input expressions, the additional ones are used as parameters if (inputExpressions.length > 1) { - Outputs.asArcGrid(raster, inputExpressions(1).eval(input).asInstanceOf[Int]) + RasterOutputs.asArcGrid(raster, inputExpressions(1).eval(input).asInstanceOf[Int]) } else { - Outputs.asArcGrid(raster, -1) + RasterOutputs.asArcGrid(raster, -1) } } diff --git a/sql/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala b/sql/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala index bc189785..3a6f3e20 100644 --- a/sql/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala +++ b/sql/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala @@ -383,5 +383,58 @@ class rasteralgebraTest extends TestBaseScala with BeforeAndAfter with GivenWhen val binaryDf = rasterDf.selectExpr("RS_FromArcInfoAsciiGrid(arc) as raster", "RS_FromArcInfoAsciiGrid(arc2) as raster2") assertEquals(rasterDf.count(), binaryDf.count()) } + + it("Passed RS_Metadata") { + val df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster/test1.tiff") + val result = df.selectExpr("RS_Metadata(RS_FromGeoTiff(content))").first().getSeq(0) + df.selectExpr("RS_Metadata(RS_FromGeoTiff(content))").show(false) + assertEquals(10, result.length) + assertEquals(512.0, result(2), 0.001) + assertEquals(517.0, result(3), 0.001) + assertEquals(1.0, result(9), 0.001) + } + + it("Passed RS_MakeEmptyRaster") { + val widthInPixel = 10 + val heightInPixel = 10 + val upperLeftX = 0.0 + val upperLeftY = 0.0 + val cellSize = 1.0 + val numBands = 2 + sparkSession.sql("SELECT RS_MakeEmptyRaster(10, 10, 0.0, 0.0, 1.0, 2)").show() + // Test with numBands + var result = sparkSession.sql(s"SELECT RS_Metadata(RS_MakeEmptyRaster($widthInPixel, $heightInPixel, $upperLeftX, $upperLeftY, $cellSize, $numBands))").first().getSeq(0) + assertEquals(numBands, result(9), 0.001) + // Test without numBands + result = sparkSession.sql(s"SELECT RS_Metadata(RS_MakeEmptyRaster($widthInPixel, $heightInPixel, $upperLeftX, $upperLeftY, $cellSize))").first().getSeq(0) + assertEquals(1, result(9), 0.001) + } + + it("Passed RS_BandAsArray") { + val df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster/test1.tiff") + val metadata = df.selectExpr("RS_Metadata(RS_FromGeoTiff(content))").first().getSeq(0) + val width = metadata(2).asInstanceOf[Double].toInt + val height = metadata(3).asInstanceOf[Double].toInt + val result = df.selectExpr("RS_BandAsArray(RS_FromGeoTiff(content), 1)").first().getSeq(0) + assertEquals(width * height, result.length) + val resultOutOfBound = df.selectExpr("RS_BandAsArray(RS_FromGeoTiff(content), 100)").first().isNullAt(0) + assertEquals(true, resultOutOfBound) + } + + it("Passed RS_AddBandFromArray") { + var df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster/test1.tiff") + df = df.selectExpr("RS_FromGeoTiff(content) as raster", "RS_MultiplyFactor(RS_BandAsArray(RS_FromGeoTiff(content), 1), 2) as band") + val bandNewExpected:Seq[Double] = df.selectExpr("band").first().getSeq(0) + val bandNewActual:Seq[Double] = df.selectExpr("RS_BandAsArray(RS_AddBandFromArray(raster, band, 1), 1)").first().getSeq(0) + df.selectExpr("RS_BandAsArray(RS_AddBandFromArray(raster, band, 1), 1) as band").show() + for (i <- bandNewExpected.indices) { + // The band value needs to be mod 256 because the ColorModel will mod 256. + assertEquals(bandNewExpected(i)%256, bandNewActual(i), 0.001) + } + // Test with out of bound band index. It should fail. + intercept[Exception] { + df.selectExpr("RS_BandAsArray(RS_AddBandFromArray(raster, band, 100), 1)").collect() + } + } } } diff --git a/sql/common/src/test/scala/org/apache/sedona/sql/serdeAwareTest.scala b/sql/common/src/test/scala/org/apache/sedona/sql/serdeAwareTest.scala index 27f0ca54..786af274 100644 --- a/sql/common/src/test/scala/org/apache/sedona/sql/serdeAwareTest.scala +++ b/sql/common/src/test/scala/org/apache/sedona/sql/serdeAwareTest.scala @@ -20,7 +20,7 @@ package org.apache.sedona.sql import org.apache.sedona.common.geometrySerde.GeometrySerializer -import org.apache.sedona.common.raster.Constructors.fromArcInfoAsciiGrid +import org.apache.sedona.common.raster.RasterConstructors.fromArcInfoAsciiGrid import org.apache.sedona.common.raster.Serde import org.apache.spark.sql.catalyst.expressions.Literal import org.apache.spark.sql.sedona_sql.expressions.{ST_Buffer, ST_GeomFromText, ST_Point, ST_Union}
