Skip to content

Commit

Permalink
[SEDONA-338] Fix RS_MakeEmptyRaster and make affine transformation …
Browse files Browse the repository at this point in the history
…parameters compliant with the GDAL/PostGIS convention (#940)
  • Loading branch information
Kontinuation authored Aug 2, 2023
1 parent c340a97 commit 8d5c76d
Show file tree
Hide file tree
Showing 9 changed files with 234 additions and 129 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -18,25 +18,15 @@
*/
package org.apache.sedona.common.raster;

import com.sun.media.imageioimpl.common.BogusColorSpace;
import org.geotools.coverage.CoverageFactoryFinder;
import org.apache.sedona.common.utils.RasterUtils;
import org.geotools.coverage.GridSampleDimension;
import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.GridCoverageFactory;

import javax.media.jai.RasterFactory;

import java.awt.Point;
import java.awt.Transparency;
import java.awt.color.ColorSpace;
import java.awt.image.BufferedImage;
import java.awt.image.ColorModel;
import java.awt.image.ComponentColorModel;
import java.awt.image.DataBuffer;
import java.awt.image.Raster;
import java.awt.image.RenderedImage;
import java.awt.image.WritableRaster;
import java.util.Arrays;

public class MapAlgebra
{
Expand Down Expand Up @@ -134,7 +124,7 @@ private static GridCoverage2D copyRasterAndAppendBand(GridCoverage2D gridCoverag
System.arraycopy(originalSampleDimensions, 0, sampleDimensions, 0, originalSampleDimensions.length);
sampleDimensions[numBand - 1] = new GridSampleDimension("band" + numBand);
// Construct a GridCoverage2D with the copied image.
return createCompatibleGridCoverage2D(gridCoverage2D, wr, sampleDimensions);
return RasterUtils.create(wr, gridCoverage2D.getGridGeometry(), sampleDimensions);
}

private static GridCoverage2D copyRasterAndReplaceBand(GridCoverage2D gridCoverage2D, int bandIndex, double[] bandValues) {
Expand All @@ -155,20 +145,6 @@ private static GridCoverage2D copyRasterAndReplaceBand(GridCoverage2D gridCovera
}
}
// Create a new GridCoverage2D with the copied image
return createCompatibleGridCoverage2D(gridCoverage2D, wr, gridCoverage2D.getSampleDimensions());
}

private static GridCoverage2D createCompatibleGridCoverage2D(GridCoverage2D gridCoverage2D, WritableRaster wr, GridSampleDimension[] bands) {
int rasterDataType = wr.getDataBuffer().getDataType();
int numBand = wr.getNumBands();
final ColorSpace cs = new BogusColorSpace(numBand);
final int[] nBits = new int[numBand];
Arrays.fill(nBits, DataBuffer.getDataTypeSize(rasterDataType));
ColorModel colorModel =
new ComponentColorModel(cs, nBits, false, true, Transparency.OPAQUE, rasterDataType);
final RenderedImage image = new BufferedImage(colorModel, wr, false, null);
GridCoverageFactory gridCoverageFactory = CoverageFactoryFinder.getGridCoverageFactory(null);
return gridCoverageFactory.create(gridCoverage2D.getName(), image, gridCoverage2D.getCoordinateReferenceSystem(),
gridCoverage2D.getGridGeometry().getGridToCRS(), bands, null, null);
return RasterUtils.create(wr, gridCoverage2D.getGridGeometry(), gridCoverage2D.getSampleDimensions());
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -18,9 +18,9 @@
*/
package org.apache.sedona.common.raster;

import org.apache.sedona.common.utils.RasterUtils;
import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.GridGeometry2D;
import org.geotools.coverage.processing.operation.Affine;
import org.geotools.coverage.grid.GridEnvelope2D;
import org.geotools.geometry.Envelope2D;
import org.geotools.referencing.CRS;
import org.geotools.referencing.crs.DefaultEngineeringCRS;
Expand All @@ -31,9 +31,7 @@
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
Expand Down Expand Up @@ -64,30 +62,21 @@ public static int getHeight(GridCoverage2D raster) {
}

public static double getUpperLeftX(GridCoverage2D raster) {
Envelope2D envelope2D = raster.getEnvelope2D();
return envelope2D.getMinX();
AffineTransform2D affine = RasterUtils.getGDALAffineTransform(raster);
return affine.getTranslateX();
}

public static double getUpperLeftY(GridCoverage2D raster) {
Envelope2D envelope2D = raster.getEnvelope2D();
return envelope2D.getMaxY();
AffineTransform2D affine = RasterUtils.getGDALAffineTransform(raster);
return affine.getTranslateY();
}

public static double getScaleX(GridCoverage2D raster) {
return getAffineTransform(raster).getScaleX();
return RasterUtils.getGDALAffineTransform(raster).getScaleX();
}

public static double getScaleY(GridCoverage2D raster) {
return getAffineTransform(raster).getScaleY();
}

private static AffineTransform2D getAffineTransform(GridCoverage2D raster) throws UnsupportedOperationException {
GridGeometry2D gridGeometry2D = raster.getGridGeometry();
MathTransform crsTransform = gridGeometry2D.getGridToCRS2D();
if (!(crsTransform instanceof AffineTransform2D)) {
throw new UnsupportedOperationException("Only AffineTransform2D is supported");
}
return (AffineTransform2D) crsTransform;
return RasterUtils.getGDALAffineTransform(raster).getScaleY();
}

public static Geometry envelope(GridCoverage2D raster) throws FactoryException {
Expand All @@ -100,40 +89,38 @@ public static Geometry envelope(GridCoverage2D raster) throws FactoryException {

/**
* Returns the metadata of a raster as an array of doubles.
* @param raster
* @param raster the raster
* @return double[] with the following values:
* 0: minX: upper left x
* 1: maxY: upper left y
* 0: upperLeftX: upper left x
* 1: upperLeftY: 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
* 6: skewX: skew on x axis
* 7: skewY: 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 Geo-reference metadata
GridEnvelope2D gridRange = raster.getGridGeometry().getGridRange2D();
AffineTransform2D affine = RasterUtils.getGDALAffineTransform(raster);

// 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");
}
// Get the affine parameters
double upperLeftX = affine.getTranslateX();
double upperLeftY = affine.getTranslateY();
double scaleX = affine.getScaleX();
double scaleY = affine.getScaleY();
double skewX = affine.getShearX();
double skewY = affine.getShearY();
return new double[] {
upperLeftX, upperLeftY,
gridRange.getWidth(), gridRange.getHeight(),
scaleX, scaleY, skewX, skewY,
srid(raster), raster.getNumSampleDimensions()};
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -14,19 +14,18 @@
package org.apache.sedona.common.raster;

import org.apache.sedona.common.raster.inputstream.ByteArrayImageInputStream;
import org.geotools.coverage.CoverageFactoryFinder;
import org.apache.sedona.common.utils.RasterUtils;
import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.GridCoverageFactory;
import org.geotools.coverage.grid.GridEnvelope2D;
import org.geotools.coverage.grid.GridGeometry2D;
import org.geotools.gce.arcgrid.ArcGridReader;
import org.geotools.gce.geotiff.GeoTiffReader;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.referencing.CRS;
import org.geotools.referencing.crs.DefaultEngineeringCRS;
import org.geotools.referencing.operation.transform.AffineTransform2D;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.datum.PixelInCell;
import org.opengis.referencing.operation.MathTransform;

import javax.media.jai.RasterFactory;
Expand All @@ -47,26 +46,28 @@ public static GridCoverage2D fromGeoTiff(byte[] bytes) throws IOException {
}

/**
* 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
* 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, minY (lowerLeftY) = upperLeftY - height * pixelSize
* 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
* SRID is default to 0 which means the default CRS (Cartesian 2D)
* <ul>
* <li>The raster is defined by the width and height
* <li>The upper left corner is defined by the upperLeftX and upperLeftY
* <li>The scale is defined by pixelSize. The scaleX is equal to pixelSize and scaleY is equal to -pixelSize
* <li>skewX and skewY are zero, which means no shear or rotation.
* <li>SRID is default to 0 which means the default CRS (Generic 2D)
* </ul>
* @param numBand the number of bands
* @param widthInPixel
* @param heightInPixel
* @param widthInPixel the width of the raster, in pixel
* @param heightInPixel the height of the raster, in pixel
* @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 upperLeftY the upper left corner of the raster. Note that: the minY of the envelope is equal to the upperLeftY - height * pixelSize
* @param pixelSize the size of the pixel in the unit of the CRS
* @return
* @return the new empty raster
*/
public static GridCoverage2D makeEmptyRaster(int numBand, int widthInPixel, int heightInPixel, double upperLeftX, double upperLeftY, double pixelSize)
throws FactoryException
{
return makeEmptyRaster(numBand, widthInPixel, heightInPixel, upperLeftX, upperLeftY, pixelSize, pixelSize, 0, 0, 0);
return makeEmptyRaster(numBand, widthInPixel, heightInPixel, upperLeftX, upperLeftY, pixelSize, -pixelSize, 0, 0, 0);
}

/**
Expand All @@ -75,13 +76,13 @@ public static GridCoverage2D makeEmptyRaster(int numBand, int widthInPixel, int
* @param widthInPixel the width of the raster, in pixel
* @param heightInPixel the height of the raster, in pixel
* @param upperLeftX the upper left corner of the raster, in the CRS unit. Note that: the minX of the envelope is equal to the upperLeftX
* @param upperLeftY the upper left corner of the raster, in the CRS unit. Note that: the minY of the envelope is equal to the upperLeftY - height * scaleY
* @param upperLeftY the upper left corner of the raster, in the CRS unit. Note that: the minY of the envelope is equal to the upperLeftY + height * scaleY
* @param scaleX the scale of the raster (pixel size on X), in the CRS unit
* @param scaleY the scale of the raster (pixel size on Y), in the CRS unit
* @param skewX the skew of the raster on X, in the CRS unit
* @param skewY the skew of the raster on Y, in the CRS unit
* @param srid the srid of the CRS. 0 means the default CRS (Cartesian 2D)
* @return
* @return the new empty raster
* @throws FactoryException
*/
public static GridCoverage2D makeEmptyRaster(int numBand, int widthInPixel, int heightInPixel, double upperLeftX, double upperLeftY, double scaleX, double scaleY, double skewX, double skewY, int srid)
Expand All @@ -93,13 +94,14 @@ public static GridCoverage2D makeEmptyRaster(int numBand, int widthInPixel, int
} else {
crs = CRS.decode("EPSG:" + srid);
}

// Create a new empty raster
WritableRaster raster = RasterFactory.createBandedRaster(DataBuffer.TYPE_DOUBLE, widthInPixel, heightInPixel, numBand, null);
MathTransform transform = new AffineTransform2D(scaleX, skewY, skewX, -scaleY, upperLeftX + scaleX / 2, upperLeftY - scaleY / 2);
GridGeometry2D gridGeometry = new GridGeometry2D(new GridEnvelope2D(0, 0, widthInPixel, heightInPixel), transform, crs);
ReferencedEnvelope referencedEnvelope = new ReferencedEnvelope(gridGeometry.getEnvelope2D());
// Create a new coverage
GridCoverageFactory gridCoverageFactory = CoverageFactoryFinder.getGridCoverageFactory(null);
return gridCoverageFactory.create("genericCoverage", raster, referencedEnvelope);
MathTransform transform = new AffineTransform2D(scaleX, skewY, skewX, scaleY, upperLeftX, upperLeftY);
GridGeometry2D gridGeometry = new GridGeometry2D(
new GridEnvelope2D(0, 0, widthInPixel, heightInPixel),
PixelInCell.CELL_CORNER,
transform, crs, null);
return RasterUtils.create(raster, gridGeometry, null);
}
}
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@
/*
* 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.utils;

import com.sun.media.imageioimpl.common.BogusColorSpace;
import org.geotools.coverage.CoverageFactoryFinder;
import org.geotools.coverage.GridSampleDimension;
import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.GridCoverageFactory;
import org.geotools.coverage.grid.GridGeometry2D;
import org.geotools.referencing.operation.transform.AffineTransform2D;
import org.opengis.metadata.spatial.PixelOrientation;
import org.opengis.referencing.operation.MathTransform;

import java.awt.Transparency;
import java.awt.color.ColorSpace;
import java.awt.image.BufferedImage;
import java.awt.image.ColorModel;
import java.awt.image.ComponentColorModel;
import java.awt.image.DataBuffer;
import java.awt.image.RenderedImage;
import java.awt.image.WritableRaster;
import java.util.Arrays;

/**
* Utility functions for working with GridCoverage2D objects.
*/
public class RasterUtils {
private RasterUtils() {}

private static final GridCoverageFactory gridCoverageFactory = CoverageFactoryFinder.getGridCoverageFactory(null);

/**
* Create a new empty raster from the given WritableRaster object.
* @param raster The raster object to be wrapped as an image.
* @param gridGeometry The grid geometry of the raster.
* @param bands The bands of the raster.
* @return A new GridCoverage2D object.
*/
public static GridCoverage2D create(WritableRaster raster, GridGeometry2D gridGeometry, GridSampleDimension[] bands) {
int numBand = raster.getNumBands();
int rasterDataType = raster.getDataBuffer().getDataType();

// Construct a color model for the rendered image. This color model should be able to be serialized and
// deserialized. The color model object automatically constructed by grid coverage factory may not be
// serializable, please refer to https://issues.apache.org/jira/browse/SEDONA-319 for more details.
final ColorSpace cs = new BogusColorSpace(numBand);
final int[] nBits = new int[numBand];
Arrays.fill(nBits, DataBuffer.getDataTypeSize(rasterDataType));
ColorModel colorModel =
new ComponentColorModel(cs, nBits, false, true, Transparency.OPAQUE, rasterDataType);

final RenderedImage image = new BufferedImage(colorModel, raster, false, null);
return gridCoverageFactory.create("genericCoverage", image, gridGeometry, bands, null, null);
}

/**
* Get a GDAL-compliant affine transform from the given raster, where the grid coordinate indicates the upper left
* corner of the pixel. PostGIS also follows GDAL convention.
* @param raster The raster to get the affine transform from.
* @return The affine transform.
*/
public static AffineTransform2D getGDALAffineTransform(GridCoverage2D raster) {
return getAffineTransform(raster, PixelOrientation.UPPER_LEFT);
}

public static AffineTransform2D getAffineTransform(GridCoverage2D raster, PixelOrientation orientation) throws UnsupportedOperationException {
GridGeometry2D gridGeometry2D = raster.getGridGeometry();
MathTransform crsTransform = gridGeometry2D.getGridToCRS2D(orientation);
if (!(crsTransform instanceof AffineTransform2D)) {
throw new UnsupportedOperationException("Only AffineTransform2D is supported");
}
return (AffineTransform2D) crsTransform;
}
}
Loading

0 comments on commit 8d5c76d

Please sign in to comment.