ucar.nc2.geotiff.GeotiffWriter Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of cdm Show documentation
The NetCDF-Java Library is a Java interface to NetCDF files, as well as to many other types of scientific data formats.
The newest version!
/*
 * Copyright 1998-2009 University Corporation for Atmospheric Research/Unidata
 *
 * Portions of this software were developed by the Unidata Program at the
 * University Corporation for Atmospheric Research.
 *
 * Access and use of this software shall impose the following obligations
 * and understandings on the user. The user is granted the right, without
 * any fee or cost, to use, copy, modify, alter, enhance and distribute
 * this software, and any derivative works thereof, and its supporting
 * documentation for any purpose whatsoever, provided that this entire
 * notice appears in all copies of the software, derivative works and
 * supporting documentation.  Further, UCAR requests that the user credit
 * UCAR/Unidata in any publications that result from the use of this
 * software or in any product that includes this software. The names UCAR
 * and/or Unidata, however, may not be used in any advertising or publicity
 * to endorse or promote any products or commercial entity unless specific
 * written permission is obtained from UCAR/Unidata. The user also
 * understands that UCAR/Unidata is not obligated to provide the user with
 * any support, consulting, training or assistance of any kind with regard
 * to the use, operation and performance of this software nor to provide
 * the user with any updates, revisions, new versions or "bug fixes."
 *
 * THIS SOFTWARE IS PROVIDED BY UCAR/UNIDATA "AS IS" AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL UCAR/UNIDATA BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING
 * FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
 * NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
 * WITH THE ACCESS, USE OR PERFORMANCE OF THIS SOFTWARE.
 */

package ucar.nc2.geotiff;

import ucar.ma2.*;
import ucar.nc2.dataset.CoordinateAxis1D;
import ucar.nc2.dt.GridCoordSystem;
import ucar.nc2.dt.GridDataset;
import ucar.nc2.dt.GridDatatype;
import ucar.unidata.geoloc.*;
import ucar.unidata.geoloc.projection.*;
import ucar.unidata.geoloc.projection.proj4.AlbersEqualAreaEllipse;

import java.io.IOException;


/**
 * Write GeoTIFF files.
 *
 * @author caron, yuan
 */
public class GeotiffWriter {

  protected GeoTiff geotiff;
  protected short pageNumber = 1;

  /**
   * Constructor
   *
   * @param fileOut name of output file.
   */
  public GeotiffWriter(String fileOut) {
    geotiff = new GeoTiff(fileOut);
  }

  /**
   * Write Grid data to the geotiff file.
   *
   * @param dataset   grid in contained in this dataset
   * @param grid      data is in this grid
   * @param data      2D array in YX order
   * @param greyScale if true, write greyScale image, else dataSample.
   * @throws IOException on i/o error
   */
  public void writeGrid(GridDataset dataset, GridDatatype grid, Array data, boolean greyScale) throws IOException {
    GridCoordSystem gcs = grid.getCoordinateSystem();

    if (!gcs.isRegularSpatial()) {
      throw new IllegalArgumentException("Must have 1D x and y axes for " + grid.getFullName());
    }

    CoordinateAxis1D xaxis = (CoordinateAxis1D) gcs.getXHorizAxis();
    CoordinateAxis1D yaxis = (CoordinateAxis1D) gcs.getYHorizAxis();

    //latlon coord does not need to be scaled
    double scaler = (gcs.isLatLon()) ? 1.0 : 1000.0;

    // data must go from top to bottom
    double xStart = xaxis.getCoordEdge(0) * scaler;
    double yStart = yaxis.getCoordEdge(0) * scaler;
    double xInc = xaxis.getIncrement() * scaler;
    double yInc = Math.abs(yaxis.getIncrement()) * scaler;

    if (yaxis.getCoordValue(0) < yaxis.getCoordValue(1)) {
      data = data.flip(0);
      yStart = yaxis.getCoordEdge((int) yaxis.getSize()) * scaler;
    }

    /*  remove - i think unneeded, monotonic lon handled in CoordinateAxis1D. JC 3/18/2013
     if (gcs.isLatLon()) {
      Array lon = xaxis.read();
      data = geoShiftDataAtLon(data, lon);
      xStart = geoShiftGetXstart(lon, xInc);
      //xStart = -180.0;
    }  */

    if (!xaxis.isRegular() || !yaxis.isRegular()) {
      throw new IllegalArgumentException("Must be evenly spaced grid = " + grid.getFullName());
    }

    if (pageNumber > 1) {
      geotiff.initTags();
    }

    // write it out
    writeGrid(grid, data, greyScale, xStart, yStart, xInc, yInc, pageNumber);
    pageNumber++;
  }


  /**
   * Write Grid data to the geotiff file.
   * Grid currently must:
   * 
   *  have a 1D X and Y coordinate axes.
   * 
 be lat/lon or Lambert Conformal Projection
   * 
 be equally spaced
   * 
   *
   * @param grid        original grid
   * @param data        2D array in YX order
   * @param greyScale   if true, write greyScale image, else dataSample.
   * @param xStart
   * @param yStart
   * @param xInc
   * @param yInc
   * @param imageNumber
   * @throws IOException              on i/o error
   * @throws IllegalArgumentException if above assumptions not valid
   */
  void writeGrid(GridDatatype grid, Array data, boolean greyScale, double xStart, double yStart, double xInc,
                        double yInc, int imageNumber) throws IOException {

    int nextStart = 0;
    GridCoordSystem gcs = grid.getCoordinateSystem();

    // get rid of this when all projections are implemented
    if (!gcs.isLatLon()
            && !(gcs.getProjection() instanceof LambertConformal)
            && !(gcs.getProjection() instanceof Stereographic)
            && !(gcs.getProjection() instanceof Mercator)
            //  && !(gcs.getProjection() instanceof TransverseMercator)
            && !(gcs.getProjection() instanceof AlbersEqualAreaEllipse)
            && !(gcs.getProjection() instanceof AlbersEqualArea)) {
      throw new IllegalArgumentException("Unsupported projection = " + gcs.getProjection().getClass().getName());
    }

    // write the data first
    if (greyScale) {
      ArrayByte result = replaceMissingValuesAndScale(grid, data);
      nextStart = geotiff.writeData((byte[]) result.getStorage(), imageNumber);
    } else {
      ArrayFloat result = replaceMissingValues(grid, data);
      nextStart = geotiff.writeData((float[]) result.getStorage(), imageNumber);
    }

    // set the width and the height
    int elemSize = greyScale ? 1 : 4;
    int height = data.getShape()[0];         // Y
    int width = data.getShape()[1];         // X
    geotiff.addTag(new IFDEntry(Tag.ImageWidth, FieldType.SHORT).setValue(width));
    geotiff.addTag(new IFDEntry(Tag.ImageLength, FieldType.SHORT).setValue(height));

    // set the multiple images tag
    int ff = 1 << 1;
    int page = imageNumber - 1;
    geotiff.addTag(new IFDEntry(Tag.NewSubfileType, FieldType.SHORT).setValue(ff));
    geotiff.addTag(new IFDEntry(Tag.PageNumber, FieldType.SHORT).setValue(page, 2));

    // just make it all one big "row"
    geotiff.addTag(new IFDEntry(Tag.RowsPerStrip, FieldType.SHORT).setValue(1));  //height));
    // the following changes to make it viewable in ARCMAP
        /*
        int size = elemSize * height * width;  // size in bytes
        geotiff.addTag( new IFDEntry(Tag.StripByteCounts, FieldType.LONG).setValue( size));
        // data starts here, header is written at the end
        if( imageNumber == 1 )
          geotiff.addTag( new IFDEntry(Tag.StripOffsets, FieldType.LONG).setValue( 8));
        else
          geotiff.addTag( new IFDEntry(Tag.StripOffsets, FieldType.LONG).setValue(nextStart));
        */
    int[] soffset = new int[height];
    int[] sbytecount = new int[height];
    if (imageNumber == 1) {
      soffset[0] = 8;
    } else {
      soffset[0] = nextStart;
    }
    sbytecount[0] = width * elemSize;
    for (int i = 1; i < height; i++) {
      soffset[i] = soffset[i - 1] + width * elemSize;
      sbytecount[i] = width * elemSize;
    }
    geotiff.addTag(new IFDEntry(Tag.StripByteCounts, FieldType.LONG, width).setValue(sbytecount));
    geotiff.addTag(new IFDEntry(Tag.StripOffsets, FieldType.LONG, width).setValue(soffset));
    // standard tags
    geotiff.addTag(new IFDEntry(Tag.Orientation, FieldType.SHORT).setValue(1));
    geotiff.addTag(new IFDEntry(Tag.Compression, FieldType.SHORT).setValue(1));  // no compression
    geotiff.addTag(new IFDEntry(Tag.Software, FieldType.ASCII).setValue("nc2geotiff"));
    geotiff.addTag(new IFDEntry(Tag.PhotometricInterpretation, FieldType.SHORT).setValue(1));  // black is zero : not used?
    geotiff.addTag(new IFDEntry(Tag.PlanarConfiguration, FieldType.SHORT).setValue(1));

    if (greyScale) {
      // standard tags for Greyscale images ( see TIFF spec, section 4)
      geotiff.addTag(new IFDEntry(Tag.BitsPerSample, FieldType.SHORT).setValue(8));  // 8 bits per sample
      geotiff.addTag(new IFDEntry(Tag.SamplesPerPixel, FieldType.SHORT).setValue(1));

      geotiff.addTag(new IFDEntry(Tag.XResolution, FieldType.RATIONAL).setValue(1, 1));
      geotiff.addTag(new IFDEntry(Tag.YResolution, FieldType.RATIONAL).setValue(1, 1));
      geotiff.addTag(new IFDEntry(Tag.ResolutionUnit, FieldType.SHORT).setValue(1));

    } else {
      // standard tags for SampleFormat ( see TIFF spec, section 19)
      geotiff.addTag(new IFDEntry(Tag.BitsPerSample, FieldType.SHORT).setValue(32));  // 32 bits per sample
      geotiff.addTag(new IFDEntry(Tag.SampleFormat, FieldType.SHORT).setValue(3));  // Sample Format
      geotiff.addTag(new IFDEntry(Tag.SamplesPerPixel, FieldType.SHORT).setValue(1));
      MAMath.MinMax dataMinMax = grid.getMinMaxSkipMissingData(data);
      float min = (float) (dataMinMax.min);
      float max = (float) (dataMinMax.max);
      geotiff.addTag(new IFDEntry(Tag.SMinSampleValue, FieldType.FLOAT).setValue(min));
      geotiff.addTag(new IFDEntry(Tag.SMaxSampleValue, FieldType.FLOAT).setValue(max));
      geotiff.addTag(new IFDEntry(Tag.GDALNoData, FieldType.FLOAT).setValue(min - 1.f));
    }

        /*
              geotiff.addTag( new IFDEntry(Tag.Geo_ModelPixelScale, FieldType.DOUBLE).setValue(
                new double[] {5.0, 2.5, 0.0} ));
              geotiff.addTag( new IFDEntry(Tag.Geo_ModelTiepoint, FieldType.DOUBLE).setValue(
                new double[] {0.0, 0.0, 0.0, -180.0, 90.0, 0.0 } ));
              //  new double[] {0.0, 0.0, 0.0, 183.0, 90.0, 0.0} ));
              IFDEntry ifd = new IFDEntry(Tag.Geo_KeyDirectory, FieldType.SHORT).setValue(
                new int[] {1, 1, 0, 4, 1024, 0, 1, 2, 1025, 0, 1, 1, 2048, 0, 1, 4326, 2054, 0, 1, 9102} );
              geotiff.addTag( ifd);
        */

    // set the transformation from projection to pixel, add tie point tag
    geotiff.setTransform(xStart, yStart, xInc, yInc);

    if (gcs.isLatLon()) {
      addLatLonTags();
    } else if (gcs.getProjection() instanceof LambertConformal) {
      addLambertConformalTags((LambertConformal) gcs.getProjection(), xStart, yStart);
    } else if (gcs.getProjection() instanceof Stereographic) {
      addPolarStereographicTags((Stereographic) gcs.getProjection(), xStart, yStart);
    } else if (gcs.getProjection() instanceof Mercator) {
      addMercatorTags((Mercator) gcs.getProjection());
    } else if (gcs.getProjection() instanceof TransverseMercator) {
      addTransverseMercatorTags((TransverseMercator) gcs.getProjection());
    } else if (gcs.getProjection() instanceof AlbersEqualArea) {
      addAlbersEqualAreaTags((AlbersEqualArea) gcs.getProjection());
    } else if (gcs.getProjection() instanceof AlbersEqualAreaEllipse) {
      addAlbersEqualAreaEllipseTags((AlbersEqualAreaEllipse) gcs.getProjection());
    }

    geotiff.writeMetadata(imageNumber);
    //geotiff.close();

  }

  public void close() throws IOException {
    geotiff.close();
  }

  /**
   * Replace missing values with dataMinMax.min - 1.0; return a floating point data array.
   *
   * @param grid GridDatatype
   * @param data input data array
   * @return floating point data array with missing values replaced.
   */
  private ArrayFloat replaceMissingValues(GridDatatype grid, Array data) {
    MAMath.MinMax dataMinMax = grid.getMinMaxSkipMissingData(data);
    float minValue = (float) (dataMinMax.min - 1.0);

    ArrayFloat floatArray = (ArrayFloat) Array.factory(float.class,
            data.getShape());
    IndexIterator dataIter = data.getIndexIterator();
    IndexIterator floatIter = floatArray.getIndexIterator();
    while (dataIter.hasNext()) {
      float v = dataIter.getFloatNext();
      if (grid.isMissingData((double) v)) {
        v = minValue;
      }
      floatIter.setFloatNext(v);
    }

    return floatArray;
  }

  /**
   * Replace missing values with 0; scale other values between 1 and 255, return a byte data array.
   *
   * @param grid GridDatatype
   * @param data input data array
   * @return byte data array with missing values replaced and data scaled from 1- 255.
   */
  private ArrayByte replaceMissingValuesAndScale(GridDatatype grid,
                                                 Array data) {
    MAMath.MinMax dataMinMax = grid.getMinMaxSkipMissingData(data);
    double scale = 254.0 / (dataMinMax.max - dataMinMax.min);

    ArrayByte byteArray = (ArrayByte) Array.factory(byte.class,
            data.getShape());
    IndexIterator dataIter = data.getIndexIterator();
    IndexIterator resultIter = byteArray.getIndexIterator();

    byte bv;
    while (dataIter.hasNext()) {
      double v = dataIter.getDoubleNext();
      if (grid.isMissingData(v)) {
        bv = 0;
      } else {
        int iv = (int) ((v - dataMinMax.min) * scale + 1);
        bv = (byte) (iv & 0xff);
      }
      resultIter.setByteNext(bv);
    }

    return byteArray;
  }

  private void addLatLonTags1() {
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GTModelTypeGeoKey,
            GeoKey.TagValue.ModelType_Geographic));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeogGeodeticDatumGeoKey,
            GeoKey.TagValue.GeogGeodeticDatum6267));
  }

  private void addLatLonTags() {
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GTModelTypeGeoKey, GeoKey.TagValue.ModelType_Geographic));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GTRasterTypeGeoKey, GeoKey.TagValue.RasterType_Area));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeographicTypeGeoKey, GeoKey.TagValue.GeographicType_WGS_84));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeogPrimeMeridianGeoKey, GeoKey.TagValue.GeogPrimeMeridian_GREENWICH));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeogAngularUnitsGeoKey, GeoKey.TagValue.GeogAngularUnits_DEGREE));
  }


  private void addPolarStereographicTags(Stereographic proj, double FalseEasting, double FalseNorthing) {
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GTModelTypeGeoKey, GeoKey.TagValue.ModelType_Projected));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GTRasterTypeGeoKey, GeoKey.TagValue.RasterType_Area));

    // define the "geographic Coordinate System"
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeographicTypeGeoKey, GeoKey.TagValue.GeographicType_WGS_84));
    //geotiff.addGeoKey( new GeoKey( GeoKey.Tag.GeogPrimeMeridianGeoKey, GeoKey.TagValue.GeogPrimeMeridian_GREENWICH));
    //geotiff.addGeoKey( new GeoKey( GeoKey.Tag.GeogAngularUnitsGeoKey, GeoKey.TagValue.GeogAngularUnits_DEGREE));

    // define the "coordinate transformation"
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjectedCSTypeGeoKey, GeoKey.TagValue.ProjectedCSType_UserDefined));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.PCSCitationGeoKey, "Snyder"));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjectionGeoKey, GeoKey.TagValue.ProjectedCSType_UserDefined));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjLinearUnitsGeoKey, GeoKey.TagValue.ProjLinearUnits_METER));
    //geotiff.addGeoKey( new GeoKey( GeoKey.Tag.ProjLinearUnitsSizeGeoKey, 1.0)); // units of km

    // the specifics for Polar Stereographic
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjCoordTransGeoKey, GeoKey.TagValue.ProjCoordTrans_Stereographic));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjCenterLongGeoKey, 0.0));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjNatOriginLatGeoKey, 90.0));
    //geotiff.addGeoKey( new GeoKey( GeoKey.Tag.ProjNatOriginLongGeoKey, proj.getTangentLon()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjScaleAtNatOriginGeoKey, 1.0));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjFalseEastingGeoKey, 0.0));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjFalseNorthingGeoKey, 0.0));
  }

  private void addLambertConformalTags(LambertConformal proj, double FalseEasting, double FalseNorthing) {
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GTModelTypeGeoKey, GeoKey.TagValue.ModelType_Projected));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GTRasterTypeGeoKey, GeoKey.TagValue.RasterType_Area));

    // define the "geographic Coordinate System"
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeographicTypeGeoKey, GeoKey.TagValue.GeographicType_WGS_84));
    //geotiff.addGeoKey( new GeoKey( GeoKey.Tag.GeogPrimeMeridianGeoKey, GeoKey.TagValue.GeogPrimeMeridian_GREENWICH));
    //geotiff.addGeoKey( new GeoKey( GeoKey.Tag.GeogAngularUnitsGeoKey, GeoKey.TagValue.GeogAngularUnits_DEGREE));

    // define the "coordinate transformation"
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjectedCSTypeGeoKey, GeoKey.TagValue.ProjectedCSType_UserDefined));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.PCSCitationGeoKey, "Snyder"));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjectionGeoKey, GeoKey.TagValue.ProjectedCSType_UserDefined));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjLinearUnitsGeoKey, GeoKey.TagValue.ProjLinearUnits_METER));
    //geotiff.addGeoKey( new GeoKey( GeoKey.Tag.ProjLinearUnitsSizeGeoKey, 1.0)); // units of km

    // the specifics for lambert conformal
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjCoordTransGeoKey, GeoKey.TagValue.ProjCoordTrans_LambertConfConic_2SP));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjStdParallel1GeoKey, proj.getParallelOne()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjStdParallel2GeoKey, proj.getParallelTwo()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjCenterLongGeoKey, proj.getOriginLon()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjNatOriginLatGeoKey, proj.getOriginLat()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjNatOriginLongGeoKey, proj.getOriginLon()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjScaleAtNatOriginGeoKey, 1.0));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjFalseEastingGeoKey, 0.0));  // LOOK why not FalseNorthing ??
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjFalseNorthingGeoKey, 0.0));
  }

  private void addMercatorTags(Mercator proj) {
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GTModelTypeGeoKey, GeoKey.TagValue.ModelType_Projected));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GTRasterTypeGeoKey, GeoKey.TagValue.RasterType_Area));

    // define the "geographic Coordinate System"
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeographicTypeGeoKey, GeoKey.TagValue.GeographicType_WGS_84));
    // geotiff.addGeoKey( new GeoKey( GeoKey.Tag.GeogLinearUnitsGeoKey, GeoKey.TagValue.ProjLinearUnits_METER));
    // geotiff.addGeoKey( new GeoKey( GeoKey.Tag.GeogAngularUnitsGeoKey, GeoKey.TagValue.GeogAngularUnits_DEGREE));

    // define the "coordinate transformation"
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjectedCSTypeGeoKey, GeoKey.TagValue.ProjectedCSType_UserDefined));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.PCSCitationGeoKey, "Mercator"));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjectionGeoKey, GeoKey.TagValue.ProjectedCSType_UserDefined));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjLinearUnitsGeoKey, GeoKey.TagValue.ProjLinearUnits_METER));
    //geotiff.addGeoKey( new GeoKey( GeoKey.Tag.ProjLinearUnitsSizeGeoKey, 1.0)); // units of km

    // the specifics for mercator
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjCoordTransGeoKey, GeoKey.TagValue.ProjCoordTrans_Mercator));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjNatOriginLongGeoKey, proj.getOriginLon()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjNatOriginLatGeoKey, proj.getParallel()));
    //    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjStdParallel1GeoKey, proj.getParallel()));
    //   geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjScaleAtNatOriginGeoKey, 1));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjFalseEastingGeoKey, 0.0));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjFalseNorthingGeoKey, 0.0));
  }

  private void addTransverseMercatorTags(TransverseMercator proj) {
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GTModelTypeGeoKey, GeoKey.TagValue.ModelType_Projected));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GTRasterTypeGeoKey, GeoKey.TagValue.RasterType_Area));

    // define the "geographic Coordinate System"
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeographicTypeGeoKey, GeoKey.TagValue.GeographicType_WGS_84));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeogLinearUnitsGeoKey, GeoKey.TagValue.ProjLinearUnits_METER));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeogAngularUnitsGeoKey, GeoKey.TagValue.GeogAngularUnits_DEGREE));

    // define the "coordinate transformation"
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjectedCSTypeGeoKey, GeoKey.TagValue.ProjectedCSType_UserDefined));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.PCSCitationGeoKey, "Transvers Mercator"));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjectionGeoKey, GeoKey.TagValue.ProjectedCSType_UserDefined));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjLinearUnitsGeoKey, GeoKey.TagValue.ProjLinearUnits_METER));
    //geotiff.addGeoKey( new GeoKey( GeoKey.Tag.ProjLinearUnitsSizeGeoKey, 1.0)); // units of km

    // the specifics for mercator
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjCoordTransGeoKey, GeoKey.TagValue.ProjCoordTrans_TransverseMercator));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjNatOriginLatGeoKey, proj.getOriginLat()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjNatOriginLongGeoKey, proj.getTangentLon()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjScaleAtNatOriginGeoKey, proj.getScale()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjScaleAtNatOriginGeoKey, 1.0));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjFalseEastingGeoKey, 0.0));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjFalseNorthingGeoKey, 0.0));
  }

  private void addAlbersEqualAreaEllipseTags(AlbersEqualAreaEllipse proj) {
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GTModelTypeGeoKey, GeoKey.TagValue.ModelType_Projected));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GTRasterTypeGeoKey, GeoKey.TagValue.RasterType_Area));

    // define the "geographic Coordinate System"
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeographicTypeGeoKey, GeoKey.TagValue.GeographicType_WGS_84));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeogLinearUnitsGeoKey, GeoKey.TagValue.ProjLinearUnits_METER));

    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeogSemiMajorAxisGeoKey, proj.getEarth().getMajor()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeogSemiMinorAxisGeoKey, proj.getEarth().getMinor()));

    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeogAngularUnitsGeoKey, GeoKey.TagValue.GeogAngularUnits_DEGREE));

    // define the "coordinate transformation"
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjectedCSTypeGeoKey, GeoKey.TagValue.ProjectedCSType_UserDefined));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.PCSCitationGeoKey, "Albers Conial Equal Area"));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjectionGeoKey, GeoKey.TagValue.ProjectedCSType_UserDefined));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjLinearUnitsGeoKey, GeoKey.TagValue.ProjLinearUnits_METER));
    //geotiff.addGeoKey( new GeoKey( GeoKey.Tag.ProjLinearUnitsSizeGeoKey, 1.0)); // units of km

    // the specifics for mercator
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjCoordTransGeoKey, GeoKey.TagValue.ProjCoordTrans_AlbersEqualAreaEllipse));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjNatOriginLatGeoKey, proj.getOriginLat()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjNatOriginLongGeoKey, proj.getOriginLon()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjStdParallel1GeoKey, proj.getParallelOne()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjStdParallel2GeoKey, proj.getParallelTwo()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjFalseEastingGeoKey, 0.0));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjFalseNorthingGeoKey, 0.0));
  }

  private void addAlbersEqualAreaTags(AlbersEqualArea proj) {
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GTModelTypeGeoKey, GeoKey.TagValue.ModelType_Projected));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GTRasterTypeGeoKey, GeoKey.TagValue.RasterType_Area));

    // define the "geographic Coordinate System"
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeographicTypeGeoKey, GeoKey.TagValue.GeographicType_WGS_84));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeogLinearUnitsGeoKey, GeoKey.TagValue.ProjLinearUnits_METER));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.GeogAngularUnitsGeoKey, GeoKey.TagValue.GeogAngularUnits_DEGREE));

    // define the "coordinate transformation"
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjectedCSTypeGeoKey, GeoKey.TagValue.ProjectedCSType_UserDefined));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.PCSCitationGeoKey, "Albers Conial Equal Area"));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjectionGeoKey, GeoKey.TagValue.ProjectedCSType_UserDefined));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjLinearUnitsGeoKey, GeoKey.TagValue.ProjLinearUnits_METER));
    //geotiff.addGeoKey( new GeoKey( GeoKey.Tag.ProjLinearUnitsSizeGeoKey, 1.0)); // units of km

    // the specifics for mercator
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjCoordTransGeoKey, GeoKey.TagValue.ProjCoordTrans_AlbersConicalEqualArea));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjNatOriginLatGeoKey, proj.getOriginLat()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjNatOriginLongGeoKey, proj.getOriginLon()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjStdParallel1GeoKey, proj.getParallelOne()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjStdParallel2GeoKey, proj.getParallelTwo()));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjFalseEastingGeoKey, 0.0));
    geotiff.addGeoKey(new GeoKey(GeoKey.Tag.ProjFalseNorthingGeoKey, 0.0));
  }

  private void dump(Array data, int col) {
    int[] shape = data.getShape();
    Index ima = data.getIndex();

    for (int j = 0; j < shape[0]; j++) {
      float dd = data.getFloat(ima.set(j, col));
      System.out.println(j + " value= " + dd);
    }
  }

  // WTF ?? is this the seam crossing ??
  private double geoShiftGetXstart(Array lon, double inc) {
    Index ilon = lon.getIndex();
    int[] lonShape = lon.getShape();
    IndexIterator lonIter = lon.getIndexIterator();
    double xlon = 0.0;

    LatLonPoint p0 = new LatLonPointImpl(0, lon.getFloat(ilon.set(0)));
    LatLonPoint pN = new LatLonPointImpl(0, lon.getFloat(ilon.set(lonShape[0] - 1)));

    xlon = p0.getLongitude();
    while (lonIter.hasNext()) {
      float l = lonIter.getFloatNext();
      LatLonPoint pn = new LatLonPointImpl(0, l);
      if (pn.getLongitude() < xlon) {
        xlon = pn.getLongitude();
      }
    }

    if (p0.getLongitude() == pN.getLongitude()) {
      xlon = xlon - inc;
    }

    return xlon;
  }

  Array geoShiftDataAtLon(Array data, Array lon) {
    int count = 0;
    int[] shape = data.getShape();
    Index ima = data.getIndex();
    Index ilon = lon.getIndex();
    int[] lonShape = lon.getShape();
    ArrayFloat adata = new ArrayFloat(new int[]{shape[0], shape[1]});
    Index imaa = adata.getIndex();
    IndexIterator lonIter = lon.getIndexIterator();

    LatLonPoint p0 = new LatLonPointImpl(0, lon.getFloat(ilon.set(lonShape[0] - 1)));
    LatLonPoint pN = new LatLonPointImpl(0, lon.getFloat(ilon.set(0)));

    while (lonIter.hasNext()) {
      float l = lonIter.getFloatNext();
      if (l > 180.0) {
        count++;
      }
    }

    //checking if the 0 point and the N point are the same point
    int spoint = 0;
    if (p0.getLongitude() == pN.getLongitude()) {
      spoint = shape[1] - count - 1;
    } else {
      spoint = shape[1] - count;
    }

    if ((count > 0) && (shape[1] > count)) {
      for (int j = 1; j < shape[1]; j++) {
        int jj = 0;

        if (j >= count) {
          jj = j - count;
        } else {
          jj = j + spoint;
        }

        for (int i = 0; i < shape[0]; i++) {
          float dd = data.getFloat(ima.set(i, jj));
          adata.setFloat(imaa.set(i, j), dd);
        }
      }

      if (p0.getLongitude() == pN.getLongitude()) {
        for (int i = 0; i < shape[0]; i++) {
          float dd = adata.getFloat(imaa.set(i, shape[1] - 1));
          adata.setFloat(imaa.set(i, 0), dd);
        }
      }
      return adata;

    } else {
      return data;
    }
  }


}