openjump-gis / openjump

OpenJUMP, the Open Source GIS with more than one trick in its kangaroo pocket, takes the leap from svn to git. join the effort!
http://openjump.org
GNU General Public License v2.0
28 stars 14 forks source link

0.5 pixel shift in raster writing #22

Closed mukoki closed 3 years ago

mukoki commented 3 years ago

mail from Peppe

When OJ or Sextante write a new raster TIF file, they save georeferencing both in a tfw and in ModelTiepoint Tag (see for instance class org.openjump.core.rasterimage.RasterImageIO. method writeImage(File outFile, Raster raster, Envelope envelope, CellSizeXY cellSize, double noData)

Few notes about how OJ reads geolocation in a (Geo)TIF: a) it first checks ModelTiepoint Tag. It it exists, OJ reads it and ignore any sidecar tfw file b) if ModelTiepoint tag doesn't exist OJ tries to read a sidecar tfw.

Few years ago I checked the way that both OJ and Sextante were saving a TIFF and I added to both the capability to save ModelTiepoint Tag [ see classes org.openjump.core.rasterimage.RasterImageIO and es.unex.sextante.openjump.core.OpenJUMPRasterLayer(Sextante)] I left the capability to save tfw from Pirol as OpenKLEM has problems reloading auto generated raster tif files if they were without their tfw. On the other hand transition to GeoTIFF sas not complete as other GeoTIFF parameters should be implemented in GeoKeyDirectory which is quite difficult to manage. Among those parameters there is also GTRasterTypeGeoKey which defines the raster space (point versus area). If this parameter is not specified, it us used as default a RasterPixelIsPoint, which references the coordinates of the image from the center of the upper left pixel, while (possibly) OJ requires to set coordinates to the upper left corner of the upper left pixel (RasterPixelIsArea) . This problem seems more common that I suspected (see from ESRI support https://support.esri.com/en/technical-article/000003225)

Which is my solution?

Since I cannot (don't want) to work on GTRasterTypeGeoKey/GeoKeyDirectory solution I manage to define the exact position of the pixel modifying the parameter ModelTiepoint Tag on saving.

In OpenJUMP I changed the class org.openjump.core.rasterimage.RasterImageIO. method writeImage(File outFile, Raster raster, Envelope envelope, CellSizeXY cellSize, double noData)

This is new write method (I enlighted in red the modified part):

public void writeImage(File outFile, Raster raster, Envelope envelope, CellSizeXY cellSize, double noData) throws IOException {
       double east = envelope.getMinX() - (0.5 * cellSize.cellSizeX);
       double west = envelope.getMaxX() - (0.5 * cellSize.cellSizeX);
       double south = envelope.getMinY() + (0.5 * cellSize.cellSizeY);
       double north = envelope.getMaxY() + (0.5 * cellSize.cellSizeY);
       Coordinate upperLeft = new Coordinate(west, north);
      Coordinate lowerRight = new Coordinate(east, south);
      Envelope envelope2 = new Envelope(upperLeft, lowerRight);
     // Delete old .xml.aux statistics file
    final File auxXmlFile = new File(outFile.getParent(), outFile.getName() + ".aux.xml");
    if (auxXmlFile.exists() && auxXmlFile.canWrite()) {
        try {
             auxXmlFile.delete();
       } catch (final Exception ex) {
     ex.printStackTrace(System.out);
     }
   }

final SampleModel sm = raster.getSampleModel();
final ColorModel colorModel = PlanarImage.createColorModel(sm);
final BufferedImage image = new BufferedImage(colorModel, (WritableRaster) raster, false, null);

final TIFFEncodeParam param = new TIFFEncodeParam();
param.setCompression(TIFFEncodeParam.COMPRESSION_NONE);

final TIFFField[] tiffFields = new TIFFField[3];

// Cell size
tiffFields[0] = new TIFFField(GeoTiffConstants.ModelPixelScaleTag, TIFFField.TIFF_DOUBLE, 2,
new double[] { cellSize.cellSizeX, cellSize.cellSizeY });

// No data
final String noDataS = Double.toString(noData);
final byte[] bytes = noDataS.getBytes();
tiffFields[1] = new TIFFField(TiffTags.TIFFTAG_GDAL_NODATA, TIFFField.TIFF_BYTE, noDataS.length(), bytes);

// Tie point
tiffFields[2] = new TIFFField(GeoTiffConstants.ModelTiepointTag, TIFFField.TIFF_DOUBLE, 6,
new double[] { 0, 0, 0, envelope2.getMinX(), envelope2.getMaxY(), 0 });

param.setExtraFields(tiffFields);

final FileOutputStream tifOut = new FileOutputStream(outFile);
final TIFFImageEncoder encoder = (TIFFImageEncoder) TIFFCodec.createImageEncoder("tiff", tifOut, param);
encoder.encode(image);
tifOut.close();

final WorldFileHandler worldFileHandler = new WorldFileHandler(outFile.getAbsolutePath(), false);
worldFileHandler.writeWorldFile(envelope2, image.getWidth(), image.getHeight());

} 

In Sextante (OpenJUMP plugIn), in the class es.unex.sextante.openjump.core.OpenJUMPRasterLayer, method postProcess() I substituted the original method with this code:

@Override
public void postProcess() throws Exception {
    final RasterImageIO rasterImageIO = new RasterImageIO();
    rasterImageIO.writeImage(new File(m_sFilename), m_Raster, m_Layer.getWholeImageEnvelope(),
    rasterImageIO.new CellSizeXY(getLayerCellSize(), getLayerCellSize()), getNoDataValue());
}

See that the new method simply calls RasterImageIO.writeImage

OpenKLEM has not been modified

I tested the first and it is working fine. This is a screenshot of a query on two different output raster (compared to the original DEM) using both Sextante and OpenKLEM, the first two rows define the upper left coordinates of the upper left corner and they are just the same

jratike80 commented 3 years ago

Probably the implementation is right but I would like to clarify the description.

In the original GeoTIFF specification http://geotiff.maptools.org/spec/geotiff2.5.html#2.5.2 the default raster space is PixelIsArea

"PixelIsArea" Raster Space The "PixelIsArea" raster grid space R, which is the default...

In the new OGC GeoTIFF standard http://docs.opengeospatial.org/is/19-008r4/19-008r4.html#_requirements_class_gtrastertypegeokey it is recommended to use the key

The use of this geokey is highly recommended for accurate georeferencing of raster.

and it is also recommended that it is explicitly set to either PixelIsArea or PixelIsPoint

http://www.opengis.net/spec/GeoTIFF/1.1/req/GTRasterTypeGeoKey.value The GTRasterTypeGeoKey value SHALL be:

  • 0 to indicate that the Raster type is undefined or unknown; or
  • 1 to indicate that the Raster type is PixelIsArea; or
  • 2 to indicate that the Raster type is PixelIsPoint; or
  • 32767 to indicate that the Raster type is user-defined. Recommendation: the use of 0 (undefined) or 32767 (user-defined) is not recommended

If I understand right, this solution does not write the GTRasterTypeGeoKey value and then the TIFF readers will interpret that the image is of PixelIsArea. If there is also a half pixel shift between TIFF TiePoint and the coordinates in the TFW the result should be right.

However, I am not sure how OpenJUMP behaves when it opens a GeoTIFF that is PixelIsPoint. If it does not write GTRasterTypeGeoKey tag, does is still read it?