// ------------------------------- //
// -------- Start of File -------- //
// ------------------------------- //
// ----------------------------------------------------------- // 
// C++ Source Code File
// Compiler Used: MSVC, GCC
// Produced By: Douglas.Gaer@noaa.gov
// File Creation Date: 06/30/2011
// Date Last Modified: 04/22/2013
// ----------------------------------------------------------- // 
// ------------- Program Description and Details ------------- // 
// ----------------------------------------------------------- // 
/*
This software and documentation was produced within NOAA 
and is intended for internal agency use only. Please don't 
distribute any part of this software or documentation without 
first contacting the original author or NOAA agency that 
produced this software. All third party libraries used to 
build this application are subject to the licensing agreement 
stated within the source code and any documentation supplied 
with the third party library.

Interpolation routines for re-projected AWIPS grids.

Interpolation example:

Find the 4 neighbors of each grid point from the projected grid, for example:

              *  *  *
              *  *  *
  north west  *-----* north east
              |  *  |
  south west  *-----* south east
              *  *  * 
              *  *  *

We will traverse the grid starting from the southwest LAT/LON.

This program offers 3 interpolation methods and 1 testing mode:

- NONE: Will not apply any interpolation, used for testing only. 

- NEAREST: Find the nearest neighbor of each grid point and rotate 
           each data point counterclockwise wise along the X axis 
           of the output grid.

- LINEAR: Find 4 neighbors of each grid point and take an average 
          of the grid point plus its 4 neighbors.

- BILINEAR: Find 4 neighbors of each grid point a determine a weight 
            value based on the distance to each neighbor. The distance, 
            in KM, cannot exceed the grid spacing in the X or Y direction. 
            The weight value is inversely proportional to distance divided 
            by the sum of all 4 weights resulting in a 1 or 0 value. The 
            output point is equal to the sum of each grid point value 
            multiplied by its weight.                                          

NOTE: For wind and current we can only interpolate U/V vectors.
NOTE: Speed (MAG) and direction (DIR) must be converted for U/V 
NOTE: before interploating.

*/
// ----------------------------------------------------------- // 

// STDLIB includes
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

// 3pLIBS includes
#include "dfileb.h"
#include "gxstring.h"

// Our project includes
#include "interpolate.h"

void InterpolateGrid::Reset()
{
  error_string.Clear(); 
  missing_value = 9999.0;
  itype = BILINEAR_INTERP;
  interpolated_grid = 0;
  if(interpolated_grid) delete [] interpolated_grid;
  interpolated_grid = 0;
}

int InterpolateGrid::WriteInterpolatedGrid(netCDFVariables &ncv, const gxString &bin_filename, gxString &out_filename)
{
  if(ncv.projectionType == "LATLON") {
    return 1;
  }

  if(ncv.projectionType == "MERCATOR") {
    return write_interpolated_mercator_grid(ncv, bin_filename, out_filename);
  }
  if(ncv.projectionType == "LAMBERT_CONFORMAL") {
    return write_interpolated_lambert_grid(ncv, bin_filename, out_filename);
  }
  if(ncv.projectionType == "POLAR_STEREOGRAPHIC") {
    return write_interpolated_polar_grid(ncv, bin_filename, out_filename);
  }

  error_string << clear << "ERROR - Unknown AWIPS grid or type not supported";
  return 0;
}

PointData *InterpolateGrid::create_latlon_grid(netCDFVariables &ncv)
{
  return CreateLatLonGrid(ncv, missing_value);
}

int InterpolateGrid::write_interpolated_lambert_grid(netCDFVariables &ncv, const gxString &bin_filename, gxString &out_filename)
{
  return write_interpolated_grid(ncv, bin_filename, out_filename);
}

int InterpolateGrid::write_interpolated_mercator_grid(netCDFVariables &ncv, const gxString &bin_filename, gxString &out_filename)
{
  return write_interpolated_grid(ncv, bin_filename, out_filename);
}


int InterpolateGrid::write_interpolated_polar_grid(netCDFVariables &ncv, const gxString &bin_filename, gxString &out_filename)
{
  return write_interpolated_grid(ncv, bin_filename, out_filename);
}


PointData *CreateLatLonGrid(netCDFVariables &ncv, float missing_value)
{
  unsigned nx = (unsigned)ncv.gridSize[0];
  unsigned ny = (unsigned)ncv.gridSize[1];

  const unsigned gridsize = nx * ny;
  PointData *latlon_grid = new PointData[gridsize];
  unsigned x = 0;
  unsigned y = 0;
  unsigned curr_point_number = 0;
  float lat = ncv.domain.SOUTHWESTLAT;
  float lon = ncv.domain.SOUTHWESTLON;

  // Create a LATLON grid to re-project and/or interploate projected grid
  for(x = 0; x < nx; x++) {
    for(y = 0; y < ny; y++) {
      latlon_grid[curr_point_number].lat = lat;
      latlon_grid[curr_point_number].lon = lon;
      latlon_grid[curr_point_number].x = x;
      latlon_grid[curr_point_number].y = y;
      latlon_grid[curr_point_number].grid_point_number = curr_point_number;
      latlon_grid[curr_point_number].data = missing_value;
      lat += ncv.domain.NSR;
      curr_point_number++;
    }
    lon += ncv.domain.EWR;
    lat = ncv.domain.SOUTHWESTLAT;
  }

  return latlon_grid;
}

int InterpolateGrid::write_interpolated_grid(netCDFVariables &ncv, const gxString &bin_filename, gxString &out_filename)
{
  DiskFileB bin_fp, interpolated_fp;
  error_string.Clear();
  error_string.Precision(4);

  bin_fp.df_Open(bin_filename.c_str());
  if(bin_fp.df_GetError() != DiskFileB::df_NO_ERROR) {
    error_string << "ERROR - Cannot open file " <<  bin_filename.c_str();
    return 0;
  }

  float data;
  unsigned i;
  unsigned nx = (unsigned)ncv.gridSize[0];
  unsigned ny = (unsigned)ncv.gridSize[1];
  int error_flag = 0;
  PointData *awips_grid = ncv.grid->GetPointData(); 

  interpolated_fp.df_Create(out_filename.c_str());
  if(interpolated_fp.df_GetError() != DiskFileB::df_NO_ERROR) {
    error_string << "ERROR - Cannot create file " <<  out_filename.c_str();
    bin_fp.df_Close();
    delete[] interpolated_grid;
    return 0;
  }

  unsigned x, y;
  unsigned current_point = 0;
  unsigned point_location = 0;

  while(!bin_fp.df_EOF()) {
    // Load the first hour of data from the AWIPS grid
    current_point = 0;
    for(x = 0; x < nx; x++) {
      for(y = 0; y < ny; y++) {
	bin_fp.df_Read((unsigned char *)&data, sizeof(data));
	if(bin_fp.df_GetError() != DiskFileB::df_NO_ERROR) {
	  error_string << "ERROR - Cannot read from file " 
		       <<  bin_filename.c_str();
	  error_flag = 1;
	  break;
	}
	point_location = (x * ny) + y;
	if(point_location > ((nx * ny)-1)) point_location = (nx * ny)-1;

	// Set the data points in our AWIPS grid
	awips_grid[point_location].data = data;
	if(FloatEqualTo(awips_grid[point_location].data, ncv.fillValue)) awips_grid[point_location].data = missing_value;

	// Reset the data points in the interpolated grid
	interpolated_grid[current_point].data = missing_value;
	current_point++;
      }
      if(error_flag > 0) break;
    }

    if(error_flag > 0) break;

    int rv;
    if(ncv.projectionType == "MERCATOR") {
      rv = interpolate_mercator_awips_grid(ncv);
    }
    else if(ncv.projectionType == "LAMBERT_CONFORMAL") {
      rv = interpolate_lambert_awips_grid(ncv);
    }
    else if(ncv.projectionType == "POLAR_STEREOGRAPHIC") {
      rv = interpolate_polar_awips_grid(ncv);
    }
    else {
      error_string << "ERROR - Cannot interpolate projection " <<  ncv.projectionType.c_str();
      error_flag = 1;
      break;
    }

    if(!rv) {
      error_flag = 1;
      break;
    }
    
    for(i = 0; i < (nx*ny); i++) {
      interpolated_fp.df_Write(&interpolated_grid[i].data, sizeof(interpolated_grid[i].data));
      if(interpolated_fp.df_GetError() != DiskFileB::df_NO_ERROR) {
	error_string << "ERROR - Cannot write to file " <<  out_filename.c_str();
	error_flag = 1;
	break;
      }
    }

    if(error_flag > 0) break;
    interpolated_fp.df_Flush(); // Flush this hour to disk

    // This hour of data is complete, Continue to next hour or break after all hours are complete
  }
  
  bin_fp.df_Close();
  interpolated_fp.df_Close();
  if(error_flag > 0) return 0;
  
  return 1;
}

int InterpolateGrid::interpolate_polar_awips_grid(netCDFVariables &ncv)
{
  unsigned nx = (unsigned)ncv.gridSize[0];
  unsigned ny = (unsigned)ncv.gridSize[1];
  unsigned x, y;
  unsigned point_location = 0;
  int nw_x, nw_y, ne_x, ne_y;
  int sw_x, sw_y, se_x, se_y;
  float nw_dist, sw_dist, ne_dist, se_dist;
  float d0, d1, d2, d3, d4;
  float sum_of_weights = 0;
  float nw_weight = 0;
  float sw_weight = 0;
  float ne_weight = 0;
  float se_weight = 0;
  float interpolated_point, point_average;
  unsigned neighboring_point_location = 0;
  unsigned neighboring_point_dataloc = 0;
  PointData *awips_grid = ncv.grid->GetPointData(); 

  PointData north_west_neighbor; // (x-1/y+1)
  PointData north_east_neighbor; // (x+1/y+1)
  PointData south_west_neighbor; // (x-1/y-1)
  PointData south_east_neighbor; // (x+1/y-1)
  
  float lat_origin = ncv.domain.SOUTHWESTLAT;
  float lon_origin = ncv.domain.SOUTHWESTLON;

  unsigned curr_point = 0;

  for(x = 0; x < nx; x++) {
    for(y = 0; y < ny; y++) {
      point_location = (x * ny) + y;
      if(point_location > ((nx * ny)-1)) point_location = (nx * ny)-1;

      // Calculate the data point value to overlay onto the re-mapped point location
      // This will rotate each data point counter clock-wise along the X axis
      unsigned dataloc = (y * nx) + x;
      if(dataloc > (nx*ny)-1) continue;

      // Calculate the i and j grid point values for each AWIPS point
      double j = (awips_grid[point_location].lon - lon_origin) / ncv.domain.EWR;
      double i = (awips_grid[point_location].lat - lat_origin) / ncv.domain.NSR;

      int int_i, int_j;
      if(x > ny) {
	Float2Int(i, int_i, 10000, 0);
	Float2Int(j, int_j, 10, 0);
      }
      else {
	Float2Int(i, int_i, 10000, 1);
	Float2Int(j, int_j, 10, 1);
      }
      if(int_i < 0) int_i -= int_i*2;
      if(int_j < 0) int_j -= int_j*2;

      // This is our re-mapped point location
      unsigned loc = (int_i * nx) + int_j;
      if(loc > (nx*ny)-1) continue;

      // Calculate the data point of our 4 neighbors
      nw_x = x-1;
      nw_y = y+1;
      ne_x = x+1;
      ne_y = y+1;
      sw_x = x-1;
      sw_y = y-1;
      se_x = x+1;
      se_y = y-1;

      nw_dist = sw_dist = ne_dist = se_dist = 0;
      nw_weight = 0; sw_weight = 0; ne_weight = 0; se_weight = 0;
      point_average = 0;
      interpolated_point = 0;
      sum_of_weights = 0;
      d0 = d1 = d2 = d3 = d4 = awips_grid[dataloc].data;

      // Set our grid spacing
      float max_grid_spacing = ncv.domain.awips_resolution_x_km + ncv.domain.awips_resolution_y_km;
      if(FloatEqualToOrLessThan(max_grid_spacing, 0)) {
	// Grid spacing info is missing so estimate the resolution
	max_grid_spacing = (ncv.domain.awips_resolution_x_deg * 100) + (ncv.domain.awips_resolution_y_deg * 100); 
      }

      if((nw_x > 0 && nw_y > 0) && (nw_x < (int)nx && nw_y < (int)ny)) {
	neighboring_point_location = (nw_x * ny) + nw_y;
	neighboring_point_dataloc = (nw_y * nx) + nw_x;
	if((neighboring_point_location < (nx * ny)) && (neighboring_point_dataloc < (nx*ny))) {
	  north_west_neighbor = awips_grid[neighboring_point_location];
	  if(FloatNotEqualTo(north_west_neighbor.data, missing_value)) {
	    if(FloatNotEqualTo(north_west_neighbor.data, ncv.fillValue)) {
	      d1 = awips_grid[neighboring_point_dataloc].data;
	      nw_dist = LatLonDistance(north_west_neighbor, awips_grid[point_location]); 
	      if((FloatGreaterThan(nw_dist, max_grid_spacing)) || 
		 (FloatEqualTo(d1, missing_value)) || (FloatEqualTo(d1, ncv.fillValue))) {
		nw_dist = 0;
		d1 = d0; 
	      }
	    }
	  }
	}
      }
      if((ne_x > 0 && ne_y > 0) && (ne_x < (int)nx && ne_y < (int)ny)) {
	neighboring_point_location = (ne_x * ny) + ne_y;
	neighboring_point_dataloc = (ne_y * nx) + ne_x;
	if((neighboring_point_location < (nx * ny)) && (neighboring_point_dataloc < (nx*ny))) {
	  north_east_neighbor = awips_grid[neighboring_point_location];
	  if(FloatNotEqualTo(north_east_neighbor.data, missing_value)) {
	    if(FloatNotEqualTo(north_east_neighbor.data, ncv.fillValue)) {
	      d2 = awips_grid[neighboring_point_dataloc].data;
	      ne_dist = LatLonDistance(north_east_neighbor, awips_grid[point_location]);  
	      if((FloatGreaterThan(ne_dist, max_grid_spacing)) ||
		 (FloatEqualTo(d2, missing_value)) || (FloatEqualTo(d2, ncv.fillValue))) {
		ne_dist = 0;
		d2 = d0;
	      }
	    }
	  }
	}
      }
      if((sw_x > 0 && sw_y > 0) && (sw_x < (int)nx && sw_y < (int)ny)) {
	neighboring_point_location = (sw_x * ny) + sw_y;
	neighboring_point_dataloc = (sw_y * nx) + sw_x;
	if((neighboring_point_location < (nx * ny)) && (neighboring_point_dataloc < (nx*ny))) {
	  south_west_neighbor = awips_grid[neighboring_point_location];
	  if(FloatNotEqualTo(south_west_neighbor.data, missing_value)) {
	    if(FloatNotEqualTo(south_west_neighbor.data, ncv.fillValue)) {
	      d3 = awips_grid[neighboring_point_dataloc].data;
	      sw_dist = LatLonDistance(south_west_neighbor, awips_grid[point_location]);  
	      if((FloatGreaterThan(sw_dist, max_grid_spacing)) ||
		 (FloatEqualTo(d3, missing_value)) || (FloatEqualTo(d3, ncv.fillValue))) {
		sw_dist = 0;
		d3 = d0;
	      }
	    }
	  }
	}
      }
      if((se_x > 0 && se_y > 0) && (se_x < (int)nx && se_y < (int)ny)) {
	neighboring_point_location = (se_x * ny) + se_y;
	neighboring_point_dataloc = (se_y * nx) + se_x;
	if((neighboring_point_location < (nx * ny)) && (neighboring_point_dataloc < (nx*ny))) {
	  south_east_neighbor = awips_grid[neighboring_point_location];
	  if(FloatNotEqualTo(south_east_neighbor.data, missing_value)) {
	    if(FloatNotEqualTo(south_east_neighbor.data, ncv.fillValue)) {
	      d4 = awips_grid[neighboring_point_dataloc].data;
	      se_dist = LatLonDistance(south_east_neighbor, awips_grid[point_location]);
	      if((FloatGreaterThan(se_dist, max_grid_spacing)) ||
		 (FloatEqualTo(d4, missing_value)) || (FloatEqualTo(d4, ncv.fillValue))) {
		se_dist = 0;
		d4 = d0;
	      }
	    }
	  }
	}
      }

      // Calculate weight values for bilinear interpolation
      sum_of_weights = nw_weight = sw_weight = ne_weight = se_weight = 0;
      
      if(FloatGreaterThan(nw_dist, 0)) sum_of_weights += 1/nw_dist;
      if(FloatGreaterThan(sw_dist, 0)) sum_of_weights += 1/sw_dist;
      if(FloatGreaterThan(ne_dist, 0)) sum_of_weights += 1/ne_dist;
      if(FloatGreaterThan(se_dist, 0)) sum_of_weights += 1/se_dist;
      
      if(FloatGreaterThan(nw_dist, 0) && FloatGreaterThan(sum_of_weights, 0)) {
	nw_weight = 1/nw_dist / sum_of_weights; 
      }
      if(FloatGreaterThan(sw_dist, 0) && FloatGreaterThan(sum_of_weights, 0)) {
	sw_weight = 1/sw_dist / sum_of_weights; 
      }
      if(FloatGreaterThan(ne_dist, 0) && FloatGreaterThan(sum_of_weights, 0)) {
	ne_weight = 1/ne_dist / sum_of_weights; 
      }
      if(FloatGreaterThan(se_dist, 0) && FloatGreaterThan(sum_of_weights, 0)) {
	se_weight = 1/se_dist / sum_of_weights; 
      }
      
      int points_to_average = 5;
      if((FloatEqualTo(d0, missing_value)) || (FloatEqualTo(d0, ncv.fillValue))) {
	points_to_average--;
	d0 = 0;
      }
      if((FloatEqualTo(d1, missing_value)) || (FloatEqualTo(d1, ncv.fillValue))) {
	points_to_average--;
	d1 = 0;
      }
      if((FloatEqualTo(d2, missing_value)) || (FloatEqualTo(d2, ncv.fillValue))) {
	points_to_average--;
	d2 = 0;
      }
      if((FloatEqualTo(d3, missing_value)) || (FloatEqualTo(d3, ncv.fillValue))) {
	points_to_average--;
	d3 = 0;
      }
      if((FloatEqualTo(d4, missing_value)) || (FloatEqualTo(d4, ncv.fillValue))) {
	points_to_average--;
	d4 = 0;
      }

      // Linear interpolation 
      if(points_to_average > 0) {
	point_average = (d0 + d1 + d2 + d3 + d4) / points_to_average;
      }
      else {
	point_average = missing_value;
      }

      // Bilinear interpolation
      interpolated_point = (nw_weight * d1) + (sw_weight * d2) + (ne_weight * d3) + (se_weight * d4);
      if(FloatEqualTo(interpolated_point, 0) && FloatNotEqualTo(awips_grid[dataloc].data, 0)) {
	interpolated_point = awips_grid[dataloc].data;
      }

      // Remove any fill values
      if(FloatEqualTo(interpolated_point, ncv.fillValue)) interpolated_point = missing_value;
      if(FloatEqualTo(point_average, ncv.fillValue)) point_average = missing_value;

      switch(itype) {
	case NO_INTERP:
	  if(curr_point < (nx*ny)) interpolated_grid[curr_point].data = awips_grid[curr_point].data;
	  break;
	case NEAREST_INTERP:
	  if(FloatEqualTo(interpolated_grid[loc].data, missing_value)) { 
	    interpolated_grid[loc].data = awips_grid[dataloc].data;
	  }
	  if(ncv.projectionType == "POLAR_STEREOGRAPHIC") {
	    loc++;
	    if(loc > (nx*ny)-1) break;
	    if(y > (ny/2)) {
	      if(FloatEqualTo(interpolated_grid[loc].data, missing_value)) { 
		interpolated_grid[loc].data = awips_grid[dataloc].data;
	      }
	    }
	  }
	  break;
	case LINEAR_INTERP:
	  if(FloatEqualTo(interpolated_grid[loc].data, missing_value)) {
	    interpolated_grid[loc].data = point_average;
	  }
	  loc++;
	  if(loc > (nx*ny)-1) break;
	  if(y > (ny/2)) { 
	    if(FloatEqualTo(interpolated_grid[loc].data, missing_value)) {
	      interpolated_grid[loc].data = point_average; 
	    }
	  }
      	  break;
	case BILINEAR_INTERP:
	  if(FloatEqualTo(interpolated_grid[loc].data, missing_value)) {
	    interpolated_grid[loc].data = interpolated_point;
	  }
	  loc++;
	  if(loc > (nx*ny)-1) break;
	  if(y > (ny/2)) {
	    if(FloatEqualTo(interpolated_grid[loc].data, missing_value)) { 
	      interpolated_grid[loc].data = interpolated_point; 
	    }
	  }
      	  break;
	case REMAPTEST_INTERP:
	  // Special testing mode used for viewing remapped grid points
	  if(FloatEqualTo(interpolated_grid[loc].data, missing_value)) { 
	    interpolated_grid[loc].data = 5;
	  }
	  loc++;
	  if(loc > (nx*ny)-1) break;
	  if(y > (ny/2)) {
	    if(FloatEqualTo(interpolated_grid[loc].data, missing_value)) { 
	      interpolated_grid[loc].data = 5;
	    }
	  }
	  break;
	default: // Default to nearest neighbor interploation
	  if(FloatEqualTo(interpolated_grid[loc].data, missing_value)) { 
	    interpolated_grid[loc].data = awips_grid[dataloc].data;
	  }
	  loc++;
	  if(loc > (nx*ny)-1) break;
	  if(y > (ny/2)) {
	    interpolated_grid[loc].data = awips_grid[dataloc].data;
	  }
	  break;
      }
      curr_point++;
    }
  }

  return 1;
}

int InterpolateGrid::interpolate_lambert_awips_grid(netCDFVariables &ncv)
{
  unsigned nx = (unsigned)ncv.gridSize[0];
  unsigned ny = (unsigned)ncv.gridSize[1];
  unsigned x, y;
  unsigned point_location = 0;
  int nw_x, nw_y, ne_x, ne_y;
  int sw_x, sw_y, se_x, se_y;
  float nw_dist, sw_dist, ne_dist, se_dist;
  float d0, d1, d2, d3, d4;
  float sum_of_weights = 0;
  float nw_weight = 0;
  float sw_weight = 0;
  float ne_weight = 0;
  float se_weight = 0;
  float interpolated_point, point_average;
  unsigned neighboring_point_location = 0;
  unsigned neighboring_point_dataloc = 0;
  PointData *awips_grid = ncv.grid->GetPointData(); 

  // Set our grid spacing
  float max_grid_spacing = ncv.domain.awips_resolution_x_km + ncv.domain.awips_resolution_y_km;
  if(FloatEqualToOrLessThan(max_grid_spacing, 0)) {
    // Grid spacing info is missing so estimate the resolution
    max_grid_spacing = (ncv.domain.awips_resolution_x_deg * 100) + (ncv.domain.awips_resolution_y_deg * 100); 
  }

  PointData north_west_neighbor; // (x-1/y+1)
  PointData north_east_neighbor; // (x+1/y+1)
  PointData south_west_neighbor; // (x-1/y-1)
  PointData south_east_neighbor; // (x+1/y-1)
  
  float lat_origin = ncv.domain.SOUTHWESTLAT;
  float lon_origin = ncv.domain.SOUTHWESTLON;

  unsigned curr_point = 0;

  for(x = 0; x < nx; x++) {
    for(y = 0; y < ny; y++) {
      point_location = (x * ny) + y;
      if(point_location > ((nx * ny)-1)) point_location = (nx * ny)-1;

      // Calculate the data point value to overlay onto the re-mapped point location
      // This will rotate each data point counter clock-wise along the X axis
      unsigned dataloc = (y * nx) + x;
      if(dataloc > (nx*ny)-1) continue;

      // Calculate the i and j grid point values for each AWIPS point
      double j = (awips_grid[point_location].lon - lon_origin) / ncv.domain.EWR;
      double i = (awips_grid[point_location].lat - lat_origin) / ncv.domain.NSR;

      int int_i, int_j;
      Float2Int(i, int_i, 10, 1);
      Float2Int(j, int_j, 10, 1);
      if(int_i < 0) int_i -= int_i*2;
      if(int_j < 0) int_j -= int_j*2;

      // This is our re-mapped point location for point 1
      unsigned loc1 = (int_i * nx) + int_j;
      if(loc1 > (nx*ny)-1) continue;

      Float2Int(i, int_i, 10, 0);
      Float2Int(j, int_j, 10, 0);

      // This is our re-mapped point location for point 2
      unsigned loc2 = (int_i * nx) + int_j;
      if(loc2 > (nx*ny)-1) continue;

      // NOTE: The order for remap Lambert is write LOC1 for missings
      // NOTE: and write LOC2 for data location.

      // Calculate the data point of our 4 neighbors
      nw_x = x-1;
      nw_y = y+1;
      ne_x = x+1;
      ne_y = y+1;
      sw_x = x-1;
      sw_y = y-1;
      se_x = x+1;
      se_y = y-1;

      nw_dist = sw_dist = ne_dist = se_dist = 0;
      nw_weight = 0; sw_weight = 0; ne_weight = 0; se_weight = 0;
      point_average = 0;
      interpolated_point = 0;
      sum_of_weights = 0;
      d0 = d1 = d2 = d3 = d4 = awips_grid[dataloc].data;

      if((nw_x > 0 && nw_y > 0) && (nw_x < (int)nx && nw_y < (int)ny)) {
	neighboring_point_location = (nw_x * ny) + nw_y;
	neighboring_point_dataloc = (nw_y * nx) + nw_x;
	if((neighboring_point_location < (nx * ny)) && (neighboring_point_dataloc < (nx*ny))) {
	  north_west_neighbor = awips_grid[neighboring_point_location];
	  if(FloatNotEqualTo(north_west_neighbor.data, missing_value)) {
	    if(FloatNotEqualTo(north_west_neighbor.data, ncv.fillValue)) {
	      d1 = awips_grid[neighboring_point_dataloc].data;
	      nw_dist = LatLonDistance(north_west_neighbor, awips_grid[point_location]); 
	      if((FloatGreaterThan(nw_dist, max_grid_spacing)) || 
		 (FloatEqualTo(d1, missing_value)) || (FloatEqualTo(d1, ncv.fillValue))) {
		nw_dist = 0;
		d1 = d0; 
	      }
	    }
	  }
	}
      }
      if((ne_x > 0 && ne_y > 0) && (ne_x < (int)nx && ne_y < (int)ny)) {
	neighboring_point_location = (ne_x * ny) + ne_y;
	neighboring_point_dataloc = (ne_y * nx) + ne_x;
	if((neighboring_point_location < (nx * ny)) && (neighboring_point_dataloc < (nx*ny))) {
	  north_east_neighbor = awips_grid[neighboring_point_location];
	  if(FloatNotEqualTo(north_east_neighbor.data, missing_value)) {
	    if(FloatNotEqualTo(north_east_neighbor.data, ncv.fillValue)) {
	      d2 = awips_grid[neighboring_point_dataloc].data;
	      ne_dist = LatLonDistance(north_east_neighbor, awips_grid[point_location]);  
	      if((FloatGreaterThan(ne_dist, max_grid_spacing)) ||
		 (FloatEqualTo(d2, missing_value)) || (FloatEqualTo(d2, ncv.fillValue))) {
		ne_dist = 0;
		d2 = d0;
	      }
	    }
	  }
	}
      }
      if((sw_x > 0 && sw_y > 0) && (sw_x < (int)nx && sw_y < (int)ny)) {
	neighboring_point_location = (sw_x * ny) + sw_y;
	neighboring_point_dataloc = (sw_y * nx) + sw_x;
	if((neighboring_point_location < (nx * ny)) && (neighboring_point_dataloc < (nx*ny))) {
	  south_west_neighbor = awips_grid[neighboring_point_location];
	  if(FloatNotEqualTo(south_west_neighbor.data, missing_value)) {
	    if(FloatNotEqualTo(south_west_neighbor.data, ncv.fillValue)) {
	      d3 = awips_grid[neighboring_point_dataloc].data;
	      sw_dist = LatLonDistance(south_west_neighbor, awips_grid[point_location]);  
	      if((FloatGreaterThan(sw_dist, max_grid_spacing)) ||
		 (FloatEqualTo(d3, missing_value)) || (FloatEqualTo(d3, ncv.fillValue))) {
		sw_dist = 0;
		d3 = d0;
	      }
	    }
	  }
	}
      }
      if((se_x > 0 && se_y > 0) && (se_x < (int)nx && se_y < (int)ny)) {
	neighboring_point_location = (se_x * ny) + se_y;
	neighboring_point_dataloc = (se_y * nx) + se_x;
	if((neighboring_point_location < (nx * ny)) && (neighboring_point_dataloc < (nx*ny))) {
	  south_east_neighbor = awips_grid[neighboring_point_location];
	  if(FloatNotEqualTo(south_east_neighbor.data, missing_value)) {
	    if(FloatNotEqualTo(south_east_neighbor.data, ncv.fillValue)) {
	      d4 = awips_grid[neighboring_point_dataloc].data;
	      se_dist = LatLonDistance(south_east_neighbor, awips_grid[point_location]);
	      if((FloatGreaterThan(se_dist, max_grid_spacing)) ||
		 (FloatEqualTo(d4, missing_value)) || (FloatEqualTo(d4, ncv.fillValue))) {
		se_dist = 0;
		d4 = d0;
	      }
	    }
	  }
	}
      }

      // Calculate weight values for bilinear interpolation
      sum_of_weights = nw_weight = sw_weight = ne_weight = se_weight = 0;
      
      if(FloatGreaterThan(nw_dist, 0)) sum_of_weights += 1/nw_dist;
      if(FloatGreaterThan(sw_dist, 0)) sum_of_weights += 1/sw_dist;
      if(FloatGreaterThan(ne_dist, 0)) sum_of_weights += 1/ne_dist;
      if(FloatGreaterThan(se_dist, 0)) sum_of_weights += 1/se_dist;
      
      if(FloatGreaterThan(nw_dist, 0) && FloatGreaterThan(sum_of_weights, 0)) {
	nw_weight = 1/nw_dist / sum_of_weights; 
      }
      if(FloatGreaterThan(sw_dist, 0) && FloatGreaterThan(sum_of_weights, 0)) {
	sw_weight = 1/sw_dist / sum_of_weights; 
      }
      if(FloatGreaterThan(ne_dist, 0) && FloatGreaterThan(sum_of_weights, 0)) {
	ne_weight = 1/ne_dist / sum_of_weights; 
      }
      if(FloatGreaterThan(se_dist, 0) && FloatGreaterThan(sum_of_weights, 0)) {
	se_weight = 1/se_dist / sum_of_weights; 
      }
      
      int points_to_average = 5;
      if((FloatEqualTo(d0, missing_value)) || (FloatEqualTo(d0, ncv.fillValue))) {
	points_to_average--;
	d0 = 0;
      }
      if((FloatEqualTo(d1, missing_value)) || (FloatEqualTo(d1, ncv.fillValue))) {
	points_to_average--;
	d1 = 0;
      }
      if((FloatEqualTo(d2, missing_value)) || (FloatEqualTo(d2, ncv.fillValue))) {
	points_to_average--;
	d2 = 0;
      }
      if((FloatEqualTo(d3, missing_value)) || (FloatEqualTo(d3, ncv.fillValue))) {
	points_to_average--;
	d3 = 0;
      }
      if((FloatEqualTo(d4, missing_value)) || (FloatEqualTo(d4, ncv.fillValue))) {
	points_to_average--;
	d4 = 0;
      }

      // Linear interpolation 
      if(points_to_average > 0) {
	point_average = (d0 + d1 + d2 + d3 + d4) / points_to_average;
      }
      else {
	point_average = missing_value;
      }

      // Bilinear interpolation
      interpolated_point = (nw_weight * d1) + (sw_weight * d2) + (ne_weight * d3) + (se_weight * d4);
      if(FloatEqualTo(interpolated_point, 0) && FloatNotEqualTo(awips_grid[dataloc].data, 0)) {
	interpolated_point = awips_grid[dataloc].data;
      }

      // Remove any fill values
      if(FloatEqualTo(interpolated_point, ncv.fillValue)) interpolated_point = missing_value;
      if(FloatEqualTo(point_average, ncv.fillValue)) point_average = missing_value;

      switch(itype) {
	case NO_INTERP:
	  if(curr_point < (nx*ny)) interpolated_grid[curr_point].data = awips_grid[curr_point].data;
	  break;
	case NEAREST_INTERP:
	  if(FloatEqualTo(interpolated_grid[loc1].data, missing_value)) { 
	    interpolated_grid[loc1].data = awips_grid[dataloc].data;
	  }
	  if(FloatEqualTo(interpolated_grid[loc2].data, missing_value)) {  
	    interpolated_grid[loc2].data = awips_grid[dataloc].data;
	  }
	  break;
	case LINEAR_INTERP:
	  if(FloatEqualTo(interpolated_grid[loc1].data, missing_value)) { 
	    interpolated_grid[loc1].data = point_average;
	  }
	  if(FloatEqualTo(interpolated_grid[loc2].data, missing_value)) { 
	    interpolated_grid[loc2].data = point_average;
	  }
	  break;
	case BILINEAR_INTERP:
	  if(FloatEqualTo(interpolated_grid[loc1].data, missing_value)) { 
	    interpolated_grid[loc1].data = interpolated_point;
	  }
	  if(FloatEqualTo(interpolated_grid[loc2].data, missing_value)) { 
	    interpolated_grid[loc2].data = interpolated_point;
	  }
	  break;
	case REMAPTEST_INTERP:
	  // Special testing mode used for viewing remapped grid points
	  if(FloatEqualTo(interpolated_grid[loc1].data, missing_value)) { 
	    interpolated_grid[loc1].data = 5.0;
	  }
	  if(FloatEqualTo(interpolated_grid[loc2].data, missing_value)) { 
	    interpolated_grid[loc2].data = 5.0;
	  }
	  break;
	default: // Default to nearest neighbor interploation
	  if(FloatEqualTo(interpolated_grid[loc1].data, missing_value)) { 
	    interpolated_grid[loc1].data = awips_grid[dataloc].data;
	  }
	  if(FloatEqualTo(interpolated_grid[loc2].data, missing_value)) { 
	    interpolated_grid[loc2].data = awips_grid[dataloc].data;
	  }
	  break;
      }
      curr_point++;
    }
  }

  return 1;
}

int InterpolateGrid::interpolate_mercator_awips_grid(netCDFVariables &ncv)
{
  unsigned nx = (unsigned)ncv.gridSize[0];
  unsigned ny = (unsigned)ncv.gridSize[1];
  unsigned x, y;
  unsigned point_location = 0;
  int nw_x, nw_y, ne_x, ne_y;
  int sw_x, sw_y, se_x, se_y;
  float nw_dist, sw_dist, ne_dist, se_dist;
  float d0, d1, d2, d3, d4;
  float sum_of_weights = 0;
  float nw_weight = 0;
  float sw_weight = 0;
  float ne_weight = 0;
  float se_weight = 0;
  float interpolated_point, point_average;
  unsigned neighboring_point_location = 0;
  unsigned neighboring_point_dataloc = 0;
  PointData *awips_grid = ncv.grid->GetPointData(); 

  PointData north_west_neighbor; // (x-1/y+1)
  PointData north_east_neighbor; // (x+1/y+1)
  PointData south_west_neighbor; // (x-1/y-1)
  PointData south_east_neighbor; // (x+1/y-1)
  
  float lat_origin = ncv.domain.SOUTHWESTLAT;
  float lon_origin = ncv.domain.SOUTHWESTLON;

  unsigned curr_point = 0;

  for(x = 0; x < nx; x++) {
    for(y = 0; y < ny; y++) {
      point_location = (x * ny) + y;
      if(point_location > ((nx * ny)-1)) point_location = (nx * ny)-1;

      // Calculate the data point value to overlay onto the re-mapped point location
      // This will rotate each data point counter clock-wise along the X axis
      unsigned dataloc = (y * nx) + x;
      if(dataloc > (nx*ny)-1) continue;

      // Calculate the i and j grid point values for each AWIPS point
      double j = (awips_grid[point_location].lon - lon_origin) / ncv.domain.EWR;
      double i = (awips_grid[point_location].lat - lat_origin) / ncv.domain.NSR;

      int epsilon_rep = 0;
      int ceil_value = 0;

      int epsilon_rep_i_loc1 = 0;
      int epsilon_rep_j_loc1 = 0;

      int epsilon_rep_i_loc2 = 0;
      int epsilon_rep_j_loc2 = 10;

      int epsilon_rep_i_loc3 = 0;
      int epsilon_rep_j_loc3 = 100;

      int epsilon_rep_i_loc4 = 0;
      int epsilon_rep_j_loc4 = 1000;

      int epsilon_rep_i_loc5 = 0;
      int epsilon_rep_j_loc5 = 10000;

      int epsilon_rep_i_loc6 = 10;
      int epsilon_rep_j_loc6 = 10;

      int epsilon_rep_i_loc7 = 100;
      int epsilon_rep_j_loc7 = 10;

      int epsilon_rep_i_loc8 = 1000;
      int epsilon_rep_j_loc8 = 10;

      int epsilon_rep_i_loc9 = 10000;
      int epsilon_rep_j_loc9 = 100;

      int epsilon_rep_i_loc10 = 100000;
      int epsilon_rep_j_loc10 = 1000;

      int int_i, int_j;
      Float2Int(i, int_i, epsilon_rep_i_loc1, ceil_value);
      Float2Int(j, int_j, epsilon_rep_j_loc1, ceil_value);
      if(int_i < 0) int_i -= int_i*2;
      if(int_j < 0) int_j -= int_j*2;
      // This is our re-mapped point location for point 1
      unsigned loc1 = (int_i * nx) + int_j;
      if(loc1 > (nx*ny)-1) continue;

      Float2Int(i, int_i, epsilon_rep_i_loc2, ceil_value);
      Float2Int(j, int_j, epsilon_rep_j_loc2, ceil_value);
      // This is our re-mapped point location for point 2
      unsigned loc2 = (int_i * nx) + int_j;
      if(loc2 > (nx*ny)-1) continue;

      Float2Int(i, int_i, epsilon_rep_i_loc3, ceil_value);
      Float2Int(j, int_j, epsilon_rep_j_loc3, ceil_value);
      // This is our re-mapped point location for point 3
      unsigned loc3 = (int_i * nx) + int_j;
      if(loc3 > (nx*ny)-1) continue;

      Float2Int(i, int_i, epsilon_rep_i_loc4, ceil_value);
      Float2Int(j, int_j, epsilon_rep_j_loc4, ceil_value);
      // This is our re-mapped point location for point 4
      unsigned loc4 = (int_i * nx) + int_j;
      if(loc4 > (nx*ny)-1) continue;

      Float2Int(i, int_i, epsilon_rep_i_loc5, ceil_value);
      Float2Int(j, int_j, epsilon_rep_j_loc5, ceil_value);
      // This is our re-mapped point location for point 5
      unsigned loc5 = (int_i * nx) + int_j;
      if(loc5 > (nx*ny)-1) continue;

      Float2Int(i, int_i, epsilon_rep_i_loc6, ceil_value);
      Float2Int(j, int_j, epsilon_rep_j_loc6, ceil_value);
      // This is our re-mapped point location for point 6
      unsigned loc6 = (int_i * nx) + int_j;
      if(loc6 > (nx*ny)-1) continue;
    
      Float2Int(i, int_i, epsilon_rep_i_loc7, ceil_value);
      Float2Int(j, int_j, epsilon_rep_j_loc7, ceil_value);
      // This is our re-mapped point location for point 7
      unsigned loc7 = (int_i * nx) + int_j;
      if(loc7 > (nx*ny)-1) continue;

      Float2Int(i, int_i, epsilon_rep_i_loc8, ceil_value);
      Float2Int(j, int_j, epsilon_rep_j_loc8, ceil_value);
      // This is our re-mapped point location for point 8
      unsigned loc8 = (int_i * nx) + int_j;
      if(loc8 > (nx*ny)-1) continue;

      Float2Int(i, int_i, epsilon_rep_i_loc9, ceil_value);
      Float2Int(j, int_j, epsilon_rep_j_loc9, ceil_value);
      // This is our re-mapped point location for point 9
      unsigned loc9 = (int_i * nx) + int_j;
      if(loc9 > (nx*ny)-1) continue;

      Float2Int(i, int_i, epsilon_rep_i_loc10, ceil_value);
      Float2Int(j, int_j, epsilon_rep_j_loc10, ceil_value);
      // This is our re-mapped point location for point 10
      unsigned loc10 = (int_i * nx) + int_j;
      if(loc10 > (nx*ny)-1) continue;

      // Calculate the data point of our 4 neighbors
      nw_x = x-1;
      nw_y = y+1;
      ne_x = x+1;
      ne_y = y+1;
      sw_x = x-1;
      sw_y = y-1;
      se_x = x+1;
      se_y = y-1;

      nw_dist = sw_dist = ne_dist = se_dist = 0;
      nw_weight = 0; sw_weight = 0; ne_weight = 0; se_weight = 0;
      point_average = 0;
      interpolated_point = 0;
      sum_of_weights = 0;
      d0 = d1 = d2 = d3 = d4 = awips_grid[dataloc].data;

      // Set our grid spacing
      float max_grid_spacing = ncv.domain.awips_resolution_x_km + ncv.domain.awips_resolution_y_km;
      if(FloatEqualToOrLessThan(max_grid_spacing, 0, epsilon_rep)) {
	// Grid spacing info is missing so estimate the resolution
	max_grid_spacing = (ncv.domain.awips_resolution_x_deg * 100) + (ncv.domain.awips_resolution_y_deg * 100); 
      }

      if((nw_x > 0 && nw_y > 0) && (nw_x < (int)nx && nw_y < (int)ny)) {
	neighboring_point_location = (nw_x * ny) + nw_y;
	neighboring_point_dataloc = (nw_y * nx) + nw_x;
	if((neighboring_point_location < (nx * ny)) && (neighboring_point_dataloc < (nx*ny))) {
	  north_west_neighbor = awips_grid[neighboring_point_location];
	  if(FloatNotEqualTo(north_west_neighbor.data, missing_value, epsilon_rep)) {
	    if(FloatNotEqualTo(north_west_neighbor.data, ncv.fillValue, epsilon_rep)) {
	      d1 = awips_grid[neighboring_point_dataloc].data;
	      nw_dist = LatLonDistance(north_west_neighbor, awips_grid[point_location]); 
	      if((FloatGreaterThan(nw_dist, max_grid_spacing, epsilon_rep)) || 
		 (FloatEqualTo(d1, missing_value, epsilon_rep)) || (FloatEqualTo(d1, ncv.fillValue, epsilon_rep))) {
		nw_dist = 0;
		d1 = d0; 
	      }
	    }
	  }
	}
      }
      if((ne_x > 0 && ne_y > 0) && (ne_x < (int)nx && ne_y < (int)ny)) {
	neighboring_point_location = (ne_x * ny) + ne_y;
	neighboring_point_dataloc = (ne_y * nx) + ne_x;
	if((neighboring_point_location < (nx * ny)) && (neighboring_point_dataloc < (nx*ny))) {
	  north_east_neighbor = awips_grid[neighboring_point_location];
	  if(FloatNotEqualTo(north_east_neighbor.data, missing_value, epsilon_rep)) {
	    if(FloatNotEqualTo(north_east_neighbor.data, ncv.fillValue, epsilon_rep)) {
	      d2 = awips_grid[neighboring_point_dataloc].data;
	      ne_dist = LatLonDistance(north_east_neighbor, awips_grid[point_location]);  
	      if((FloatGreaterThan(ne_dist, max_grid_spacing, epsilon_rep)) ||
		 (FloatEqualTo(d2, missing_value, epsilon_rep)) || (FloatEqualTo(d2, ncv.fillValue, epsilon_rep))) {
		ne_dist = 0;
		d2 = d0;
	      }
	    }
	  }
	}
      }
      if((sw_x > 0 && sw_y > 0) && (sw_x < (int)nx && sw_y < (int)ny)) {
	neighboring_point_location = (sw_x * ny) + sw_y;
	neighboring_point_dataloc = (sw_y * nx) + sw_x;
	if((neighboring_point_location < (nx * ny)) && (neighboring_point_dataloc < (nx*ny))) {
	  south_west_neighbor = awips_grid[neighboring_point_location];
	  if(FloatNotEqualTo(south_west_neighbor.data, missing_value, epsilon_rep)) {
	    if(FloatNotEqualTo(south_west_neighbor.data, ncv.fillValue, epsilon_rep)) {
	      d3 = awips_grid[neighboring_point_dataloc].data;
	      sw_dist = LatLonDistance(south_west_neighbor, awips_grid[point_location]);  
	      if((FloatGreaterThan(sw_dist, max_grid_spacing, epsilon_rep)) ||
		 (FloatEqualTo(d3, missing_value, epsilon_rep)) || (FloatEqualTo(d3, ncv.fillValue, epsilon_rep))) {
		sw_dist = 0;
		d3 = d0;
	      }
	    }
	  }
	}
      }
      if((se_x > 0 && se_y > 0) && (se_x < (int)nx && se_y < (int)ny)) {
	neighboring_point_location = (se_x * ny) + se_y;
	neighboring_point_dataloc = (se_y * nx) + se_x;
	if((neighboring_point_location < (nx * ny)) && (neighboring_point_dataloc < (nx*ny))) {
	  south_east_neighbor = awips_grid[neighboring_point_location];
	  if(FloatNotEqualTo(south_east_neighbor.data, missing_value, epsilon_rep)) {
	    if(FloatNotEqualTo(south_east_neighbor.data, ncv.fillValue, epsilon_rep)) {
	      d4 = awips_grid[neighboring_point_dataloc].data;
	      se_dist = LatLonDistance(south_east_neighbor, awips_grid[point_location]);
	      if((FloatGreaterThan(se_dist, max_grid_spacing, epsilon_rep)) ||
		 (FloatEqualTo(d4, missing_value, epsilon_rep)) || (FloatEqualTo(d4, ncv.fillValue, epsilon_rep))) {
		se_dist = 0;
		d4 = d0;
	      }
	    }
	  }
	}
      }

      // Calculate weight values for bilinear interpolation
      sum_of_weights = nw_weight = sw_weight = ne_weight = se_weight = 0;
      
      if(FloatGreaterThan(nw_dist, 0, epsilon_rep)) sum_of_weights += 1/nw_dist;
      if(FloatGreaterThan(sw_dist, 0, epsilon_rep)) sum_of_weights += 1/sw_dist;
      if(FloatGreaterThan(ne_dist, 0, epsilon_rep)) sum_of_weights += 1/ne_dist;
      if(FloatGreaterThan(se_dist, 0, epsilon_rep)) sum_of_weights += 1/se_dist;
      
      if(FloatGreaterThan(nw_dist, 0, epsilon_rep) && FloatGreaterThan(sum_of_weights, 0, epsilon_rep)) {
	nw_weight = 1/nw_dist / sum_of_weights; 
      }
      if(FloatGreaterThan(sw_dist, 0, epsilon_rep) && FloatGreaterThan(sum_of_weights, 0, epsilon_rep)) {
	sw_weight = 1/sw_dist / sum_of_weights; 
      }
      if(FloatGreaterThan(ne_dist, 0, epsilon_rep) && FloatGreaterThan(sum_of_weights, 0, epsilon_rep)) {
	ne_weight = 1/ne_dist / sum_of_weights; 
      }
      if(FloatGreaterThan(se_dist, 0, epsilon_rep) && FloatGreaterThan(sum_of_weights, 0, epsilon_rep)) {
	se_weight = 1/se_dist / sum_of_weights; 
      }
      
      int points_to_average = 5;
      if((FloatEqualTo(d0, missing_value, epsilon_rep)) || (FloatEqualTo(d0, ncv.fillValue, epsilon_rep))) {
	points_to_average--;
	d0 = 0;
      }
      if((FloatEqualTo(d1, missing_value, epsilon_rep)) || (FloatEqualTo(d1, ncv.fillValue, epsilon_rep))) {
	points_to_average--;
	d1 = 0;
      }
      if((FloatEqualTo(d2, missing_value, epsilon_rep)) || (FloatEqualTo(d2, ncv.fillValue, epsilon_rep))) {
	points_to_average--;
	d2 = 0;
      }
      if((FloatEqualTo(d3, missing_value, epsilon_rep)) || (FloatEqualTo(d3, ncv.fillValue, epsilon_rep))) {
	points_to_average--;
	d3 = 0;
      }
      if((FloatEqualTo(d4, missing_value, epsilon_rep)) || (FloatEqualTo(d4, ncv.fillValue, epsilon_rep))) {
	points_to_average--;
	d4 = 0;
      }

      // Linear interpolation 
      if(points_to_average > 0) {
	point_average = (d0 + d1 + d2 + d3 + d4) / points_to_average;
      }
      else {
	point_average = missing_value;
      }

      // Bilinear interpolation
      interpolated_point = (nw_weight * d1) + (sw_weight * d2) + (ne_weight * d3) + (se_weight * d4);
      if(FloatEqualTo(interpolated_point, 0, epsilon_rep) && 
	 FloatNotEqualTo(awips_grid[dataloc].data, 0, epsilon_rep)) {
	interpolated_point = awips_grid[dataloc].data;
      }

      // Remove any fill values
      if(FloatEqualTo(interpolated_point, ncv.fillValue, epsilon_rep)) interpolated_point = missing_value;
      if(FloatEqualTo(point_average, ncv.fillValue, epsilon_rep)) point_average = missing_value;

      switch(itype) {
	case NO_INTERP:
	  if(curr_point < (nx*ny)) interpolated_grid[curr_point].data = awips_grid[curr_point].data;
	  break;
	case NEAREST_INTERP:
	  if(FloatEqualTo(interpolated_grid[loc1].data, missing_value)) 
	    interpolated_grid[loc1].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc2].data, missing_value)) 
	    interpolated_grid[loc2].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc3].data, missing_value)) 
	    interpolated_grid[loc3].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc4].data, missing_value)) 
	    interpolated_grid[loc4].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc5].data, missing_value)) 
	    interpolated_grid[loc5].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc6].data, missing_value)) 
	    interpolated_grid[loc6].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc7].data, missing_value)) 
	    interpolated_grid[loc7].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc8].data, missing_value)) 
	    interpolated_grid[loc8].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc9].data, missing_value)) 
	    interpolated_grid[loc9].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc10].data, missing_value)) 
	    interpolated_grid[loc10].data = awips_grid[dataloc].data;
	  break;

	case LINEAR_INTERP:
	  if(FloatEqualTo(interpolated_grid[loc1].data, missing_value)) 
	    interpolated_grid[loc1].data = point_average;
	  if(FloatEqualTo(interpolated_grid[loc2].data, missing_value)) 
	    interpolated_grid[loc2].data = point_average;
	  if(FloatEqualTo(interpolated_grid[loc3].data, missing_value)) 
	    interpolated_grid[loc3].data = point_average;
	  if(FloatEqualTo(interpolated_grid[loc4].data, missing_value)) 
	    interpolated_grid[loc4].data = point_average;
	  if(FloatEqualTo(interpolated_grid[loc5].data, missing_value)) 
	    interpolated_grid[loc5].data = point_average;
	  if(FloatEqualTo(interpolated_grid[loc6].data, missing_value)) 
	    interpolated_grid[loc6].data = point_average;
	  if(FloatEqualTo(interpolated_grid[loc7].data, missing_value)) 
	    interpolated_grid[loc7].data = point_average;
	  if(FloatEqualTo(interpolated_grid[loc8].data, missing_value)) 
	    interpolated_grid[loc8].data = point_average;
	  if(FloatEqualTo(interpolated_grid[loc9].data, missing_value)) 
	    interpolated_grid[loc9].data = point_average;
	  if(FloatEqualTo(interpolated_grid[loc10].data, missing_value)) 
	    interpolated_grid[loc10].data = point_average;
	  break;

	case BILINEAR_INTERP:
	  if(FloatEqualTo(interpolated_grid[loc1].data, missing_value)) 
	    interpolated_grid[loc1].data = interpolated_point;
	  if(FloatEqualTo(interpolated_grid[loc2].data, missing_value)) 
	    interpolated_grid[loc2].data = interpolated_point;
	  if(FloatEqualTo(interpolated_grid[loc3].data, missing_value)) 
	    interpolated_grid[loc3].data = interpolated_point;
	  if(FloatEqualTo(interpolated_grid[loc4].data, missing_value)) 
	    interpolated_grid[loc4].data = interpolated_point;
	  if(FloatEqualTo(interpolated_grid[loc5].data, missing_value)) 
	    interpolated_grid[loc5].data = interpolated_point;
	  if(FloatEqualTo(interpolated_grid[loc6].data, missing_value)) 
	    interpolated_grid[loc6].data = interpolated_point;
	  if(FloatEqualTo(interpolated_grid[loc7].data, missing_value)) 
	    interpolated_grid[loc7].data = interpolated_point;
	  if(FloatEqualTo(interpolated_grid[loc8].data, missing_value)) 
	    interpolated_grid[loc8].data = interpolated_point;
	  if(FloatEqualTo(interpolated_grid[loc9].data, missing_value)) 
	    interpolated_grid[loc9].data = interpolated_point;
	  if(FloatEqualTo(interpolated_grid[loc10].data, missing_value)) 
	    interpolated_grid[loc10].data = interpolated_point;
	  break;

	case REMAPTEST_INTERP:
	  // Special testing mode used for viewing remapped grid points
	  if(FloatEqualTo(interpolated_grid[loc1].data, missing_value)) interpolated_grid[loc1].data = 5;
	  if(FloatEqualTo(interpolated_grid[loc2].data, missing_value)) interpolated_grid[loc2].data = 5;
	  if(FloatEqualTo(interpolated_grid[loc3].data, missing_value)) interpolated_grid[loc3].data = 5;
	  if(FloatEqualTo(interpolated_grid[loc4].data, missing_value)) interpolated_grid[loc4].data = 5;
	  if(FloatEqualTo(interpolated_grid[loc5].data, missing_value)) interpolated_grid[loc5].data = 5;
	  if(FloatEqualTo(interpolated_grid[loc6].data, missing_value)) interpolated_grid[loc6].data = 5;
	  if(FloatEqualTo(interpolated_grid[loc7].data, missing_value)) interpolated_grid[loc7].data = 5;
	  if(FloatEqualTo(interpolated_grid[loc8].data, missing_value)) interpolated_grid[loc8].data = 5;
	  if(FloatEqualTo(interpolated_grid[loc9].data, missing_value)) interpolated_grid[loc9].data = 5;
	  if(FloatEqualTo(interpolated_grid[loc10].data, missing_value)) interpolated_grid[loc10].data = 5;
	  break;

	default: // Default to nearest neighbor interploation
	  if(FloatEqualTo(interpolated_grid[loc1].data, missing_value)) 
	    interpolated_grid[loc1].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc2].data, missing_value)) 
	    interpolated_grid[loc2].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc3].data, missing_value)) 
	    interpolated_grid[loc3].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc4].data, missing_value)) 
	    interpolated_grid[loc4].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc5].data, missing_value)) 
	    interpolated_grid[loc5].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc6].data, missing_value)) 
	    interpolated_grid[loc6].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc7].data, missing_value)) 
	    interpolated_grid[loc7].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc8].data, missing_value)) 
	    interpolated_grid[loc8].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc9].data, missing_value)) 
	    interpolated_grid[loc9].data = awips_grid[dataloc].data;
	  if(FloatEqualTo(interpolated_grid[loc10].data, missing_value)) 
	    interpolated_grid[loc10].data = awips_grid[dataloc].data;
	  break;
      }
      curr_point++;
    }
  }

  return 1;
}
// ----------------------------------------------------------- //
// ------------------------------- //
// --------- End of File --------- //
// ------------------------------- //