// $Id: GribHorizCoordSys.java,v 1.11 2004/12/14 15:41:02 caron Exp $
/*
 * Copyright 1997-2004 Unidata Program Center/University Corporation for
 * Atmospheric Research, P.O. Box 3000, Boulder, CO 80307,
 * support@unidata.ucar.edu.
 *
 * This library is free software; you can redistribute it and/or modify it
 * under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2.1 of the License, or (at
 * your option) any later version.
 *
 * This library is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
 * General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this library; if not, write to the Free Software Foundation,
 * Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 */
package ucar.nc2.iosp.grib;

import ucar.grib.*;

import ucar.ma2.Array;
import ucar.nc2.*;
import ucar.nc2.dataset.AxisType;

import ucar.unidata.geoloc.*;
import ucar.unidata.geoloc.projection.LambertConformal;
import ucar.unidata.geoloc.projection.Stereographic;
import ucar.unidata.util.StringUtil;

import java.util.ArrayList;
import java.util.HashMap;

/**
 * A horizontal coordinate system created from a Grib2GridDefinitionSection.
 * @author caron
 *
 * <p> Note on "false_easting" and "fale_northing" projection parameters:
 * <ul><li>false_easting(northing) = The value added to all x (y) values in the rectangular coordinates for a map projection.
 * This value frequently is assigned to eliminate negative numbers.
 * Expressed in the unit of measure identified in Planar Coordinate Units.
 * <li>We dont currently use, assuming that the x and y are just fine as negetive numbers.
 */
public class GribHorizCoordSys {
  private TableLookup lookup;
  private Index.GdsRecord gdsIndex;
  private Group g;

  private String grid_name, shape_name, id;
  private boolean isLatLon = true;
  HashMap varHash = new HashMap( 200); // GribVariables that have this GribHorizCoordSys
  HashMap productHash = new HashMap(100); // List of GribVariable, sorted by product desc
  HashMap vcsHash = new HashMap(30); // GribVertCoordSys

  private double startx, starty;
  private ProjectionImpl proj;
  private ArrayList attributes = new ArrayList();

  GribHorizCoordSys (Index.GdsRecord gdsIndex, TableLookup lookup, Group g) {
    this.gdsIndex = gdsIndex;
    this.lookup = lookup;
    this.g = g;

    this.grid_name = NetcdfFile.createValidNetcdfObjectName( lookup.getGridName(gdsIndex));
    this.shape_name = lookup.getShapeName(gdsIndex);
    this.g = g;
    isLatLon = lookup.isLatLon(gdsIndex);
    grid_name = StringUtil.replace( grid_name, ' ', "_")+"_projection";
    id = (g == null) ? grid_name : g.getName();
  }

  String getID() { return id; } // unique within the file
  String getGridName() { return grid_name; } // used in CF-1 attributes
  Group getGroup() { return g; }
  boolean isLatLon() { return isLatLon; }
  int getNx() { return gdsIndex.nx; }
  int getNy() { return gdsIndex.ny; }

  void addDimensionsToNetcdfFile(NetcdfFile ncfile) {

    if (isLatLon) {
      ncfile.addDimension(g, new Dimension("lat", gdsIndex.ny, true));
      ncfile.addDimension(g, new Dimension("lon", gdsIndex.nx, true));
    } else {
      ncfile.addDimension(g, new Dimension("y", gdsIndex.ny, true));
      ncfile.addDimension(g, new Dimension("x", gdsIndex.nx, true));
    }
  }

  void addToNetcdfFile(NetcdfFile ncfile) {

    if (isLatLon) {
      addCoordAxis(ncfile, "lat", gdsIndex.ny, gdsIndex.La1, gdsIndex.dy, "degrees_north", "latitude coordinate", "latitude", AxisType.Lat);
      addCoordAxis(ncfile, "lon", gdsIndex.nx, gdsIndex.Lo1, gdsIndex.dx, "degrees_east", "longitude coordinate", "longitude", AxisType.Lon);
      add2DCoordSystem(ncfile, "latLonCoordSys", "time lat lon");
    } else {
      computeProjection(ncfile);
      double[] yData = addCoordAxis(ncfile, "y", gdsIndex.ny, starty, gdsIndex.dy / 1000.0, "km",
          "y coordinate of projection", "projection_y_coordinate", AxisType.GeoY);
      double[] xData = addCoordAxis(ncfile, "x", gdsIndex.nx, startx, gdsIndex.dx / 1000.0, "km",
          "x coordinate of projection", "projection_x_coordinate", AxisType.GeoX);
      addLatLon2D( ncfile, xData, yData);
      add2DCoordSystem(ncfile, "projectionCoordSys", "time y x");
    }
  }

  private double[] addCoordAxis( NetcdfFile ncfile, String name, int n, double start, double incr, String units,
       String desc, String standard_name, AxisType axis) {

    // ncfile.addDimension(g, new Dimension(name, n, true));

    Variable v = new Variable( ncfile, g, null, name);
    v.setDataType( DataType.DOUBLE);
    v.setDimensions(name);

    // create the data
    double[] data = new double[n];
    for (int i = 0; i < n; i++) {
      data[i] = start + incr*i;
    }
    Array dataArray = Array.factory( DataType.DOUBLE.getClassType(), new int [] {n}, data);
    v.setCachedData(dataArray, false);

    v.addAttribute( new Attribute("units", units));
    v.addAttribute( new Attribute("long_name", desc));
    v.addAttribute( new Attribute("standard_name", standard_name));
    v.addAttribute( new Attribute("grid_spacing", incr+" "+units));
    v.addAttribute( new Attribute("_CoordinateAxisType", axis.toString()));

    ncfile.addVariable( g, v);
    return data;
  }

  private void addLatLon2D( NetcdfFile ncfile, double[] xData, double[] yData) {

    Variable latVar = new Variable( ncfile, g, null, "lat");
    latVar.setDataType( DataType.DOUBLE);
    latVar.setDimensions("y x");
    latVar.addAttribute( new Attribute("units", "degrees_north"));
    latVar.addAttribute( new Attribute("long_name", "latitude coordinate"));
    latVar.addAttribute( new Attribute("standard_name", "latitude"));
    latVar.addAttribute( new Attribute("_CoordinateAxisType", AxisType.Lat.toString()));

    Variable lonVar = new Variable( ncfile, g, null, "lon");
    lonVar.setDataType( DataType.DOUBLE);
    lonVar.setDimensions("y x");
    lonVar.addAttribute( new Attribute("units", "degrees_east"));
    lonVar.addAttribute( new Attribute("long_name", "longitude coordinate"));
    lonVar.addAttribute( new Attribute("standard_name", "longitude"));
    lonVar.addAttribute( new Attribute("_CoordinateAxisType", AxisType.Lon.toString()));

    int nx = xData.length;
    int ny = yData.length;

    // create the data
    ProjectionPointImpl projPoint = new ProjectionPointImpl();
    LatLonPointImpl latlonPoint = new LatLonPointImpl();
    double[] latData = new double[nx * ny];
    double[] lonData = new double[nx * ny];
    for (int i = 0; i < ny; i++) {
      for (int j = 0; j < nx; j++) {
        projPoint.setLocation(xData[j], yData[i] );
        proj.projToLatLon( projPoint, latlonPoint);
        latData[i*nx + j] = latlonPoint.getLatitude();
        lonData[i*nx + j] = latlonPoint.getLongitude();
      }
    }
    Array latDataArray = Array.factory( DataType.DOUBLE.getClassType(), new int [] {ny, nx}, latData);
    latVar.setCachedData(latDataArray, false);

    Array lonDataArray = Array.factory( DataType.DOUBLE.getClassType(), new int [] {ny, nx}, lonData);
    lonVar.setCachedData(lonDataArray, false);

    ncfile.addVariable( g, latVar);
    ncfile.addVariable( g, lonVar);
  }

  private void computeProjection(NetcdfFile ncfile) {
     switch ( lookup.getProjectionType(gdsIndex)) {
       case TableLookup.PolarStereographic:
         computePS();
         break;
       case TableLookup.LambertConformal:
         computeLC();
         break;
       default:
         throw new UnsupportedOperationException("unknown projection = "+gdsIndex.grid_type);
     }

    Variable v = new Variable( ncfile, g, null, grid_name);
    v.setDataType(DataType.CHAR);
    v.setDimensions(new ArrayList()); // scalar
    char[] data = new char[] {'d'};
    Array dataArray = Array.factory( DataType.CHAR.getClassType(), new int[0] , data);
    v.setCachedData(dataArray, false);

    for (int i = 0; i < attributes.size(); i++) {
      Attribute att = (Attribute) attributes.get(i);
      v.addAttribute( att);
    }

    v.addAttribute( new Attribute("GRIB_earth_shape", shape_name));
    v.addAttribute( new Attribute("GRIB_earth_shape_code", new Integer(gdsIndex.grid_shape_code)));
    if (gdsIndex.grid_shape_code == 1) {
      v.addAttribute( new Attribute("GRIB_spherical_earth_radius_meters", new Double(gdsIndex.radius_spherical_earth)));
    }
    ncfile.addVariable( g, v);
  }

  private void add2DCoordSystem(NetcdfFile ncfile, String name, String dims) {
    Variable v = new Variable( ncfile, g, null, name);
    v.setDataType(DataType.CHAR);
    v.setDimensions(new ArrayList()); // scalar
    Array dataArray = Array.factory( DataType.CHAR.getClassType(), new int[0] , new char[] {'0'});
    v.setCachedData(dataArray, false);
    v.addAttribute( new Attribute("_CoordinateAxes", dims));
    if (!isLatLon())
      v.addAttribute( new Attribute("_CoordinateTransforms", getGridName()));
    ncfile.addVariable( g, v);
  }

  // lambert conformal
  private void computeLC() {
    // we have to project in order to find the origin
    proj = new LambertConformal(gdsIndex.latin1, gdsIndex.LoV, gdsIndex.latin1, gdsIndex.latin2);
    ProjectionPointImpl start = (ProjectionPointImpl) proj.latLonToProj( new LatLonPointImpl( gdsIndex.La1, gdsIndex.Lo1));
    startx = start.getX();
    starty = start.getY();

    attributes.add( new Attribute("grid_mapping_name", "lambert_conformal_conic"));
    if (gdsIndex.latin1 == gdsIndex.latin2)
      attributes.add( new Attribute("standard_parallel", new Double(gdsIndex.latin1)));
    else {
      double[] data = new double[] {gdsIndex.latin1, gdsIndex.latin2};
      attributes.add( new Attribute("standard_parallel",
          Array.factory( DataType.DOUBLE.getClassType(), new int[] {2}, data)));
    }
    attributes.add( new Attribute("longitude_of_central_meridian", new Double(gdsIndex.LoV)));
    attributes.add( new Attribute("latitude_of_projection_origin", new Double(gdsIndex.latin1)));
    //attributes.add( new Attribute("false_easting", new Double(startx)));
    //attributes.add( new Attribute("false_northing", new Double(starty)));
  }

  // polar stereographic
  private void computePS() {
    double scale = .933;

    // Why the scale factor?. accordining to GRIB docs:
    // "Grid lengths are in units of meters, at the 60 degree latitude circle nearest to the pole"
    // since the scale factor at 60 degrees = k = 2*k0/(1+sin(60))  [Snyder,Working Manual p157]
    // then to make scale = 1 at 60 degrees, k0 = (1+sin(60))/2 = .933
    proj = new Stereographic(90.0, gdsIndex.LoV, scale);

    // we have to project in order to find the origin
    ProjectionPointImpl start = (ProjectionPointImpl) proj.latLonToProj( new LatLonPointImpl( gdsIndex.La1, gdsIndex.Lo1));
    if (GribServiceProvider.debugProj) System.out.println("start at proj coord "+start);
    startx = start.getX();
    starty = start.getY();

    attributes.add( new Attribute("grid_mapping_name", "polar_stereographic"));
    attributes.add( new Attribute("scale_factor_at_projection_origin", new Double(scale)));
    //attributes.add( new Attribute("false_easting", new Double(startx)));
    //attributes.add( new Attribute("false_northing", new Double(starty)));
  }

}

/* Change History:
   $Log: GribHorizCoordSys.java,v $
   Revision 1.11  2004/12/14 15:41:02  caron
   *** empty log message ***

   Revision 1.10  2004/12/07 02:43:22  caron
   *** empty log message ***

   Revision 1.9  2004/12/07 01:29:31  caron
   redo convention parsing, use _Coordinate encoding.

   Revision 1.8  2004/12/01 05:53:42  caron
   ncml pass 2, new convention parsing

   Revision 1.7  2004/11/10 17:00:28  caron
   no message

   Revision 1.6  2004/10/23 21:36:10  caron
   no message

   Revision 1.5  2004/10/13 19:45:12  caron
   add strict NCDump

   Revision 1.4  2004/10/12 02:57:05  caron
   refactor for grib1/grib2: move common functionality up to ucar.grib
   split GribServiceProvider

   Revision 1.3  2004/10/02 20:54:41  caron
   *** empty log message ***

   Revision 1.2  2004/09/30 20:49:06  caron
   *** empty log message ***

   Revision 1.1  2004/09/30 00:33:41  caron
   *** empty log message ***

*/