/****************************************************************************** FILE: interpolate.c PURPOSE: This file includes routines to handle interpolation for the surface temperature atmospheric parameter calculations. ******************************************************************************/ #include #include #include "interpolate.h" #include "utilities.h" /****************************************************************************** METHOD: interpolate_parameters PURPOSE: Interpolate atmospheric parameters to the location of the current pixel. ******************************************************************************/ void interpolate_parameters ( MODEL_POINT *model_points, /* I: results from model runs */ const GRID_POINTS *points, /* I: coordinate points */ int *cell_vertices, /* I: current cell vertices */ double interpolate_height, /* I: current landsat pixel height */ double interpolate_easting, /* I: interpolate to easting */ double interpolate_northing, /* I: interpolate to northing */ double *parameters /* O: interpolated pixel atmospheric parameters */ ) { int parameter; int elevation; int below; int above; int vertex; double at_heights[NUM_CELL_POINTS][AHP_NUM_PARAMETERS]; double below_parameters[AHP_NUM_PARAMETERS]; double above_parameters[AHP_NUM_PARAMETERS]; double w[NUM_CELL_POINTS]; double total = 0.0; double inv_total; double interp_ratio; /* interpolation ratio */ /* Interpolate three parameters to the height at each of the four closest points. */ for (vertex = 0; vertex < NUM_CELL_POINTS; vertex++) { MODEL_POINT point = model_points[cell_vertices[vertex]]; double *at_height = at_heights[vertex]; /* Find the heights between which the pixel height sits. The heights are in increasing order. */ for (elevation = 0; elevation < point.count; elevation++) { if (point.elevations[elevation].elevation >= interpolate_height) break; } if (elevation < point.count) { above = elevation; if (above != 0) below = above - 1; else /* All heights are above the pixel height, so use the first value. */ below = 0; } else { /* All heights are below the pixel height, so use the final value. */ above = point.count - 1; below = above; } if (above == below) { /* Use the below parameters since the same */ at_height[AHP_TRANSMISSION] = point.elevations[below].transmission; at_height[AHP_UPWELLED_RADIANCE] = point.elevations[below].upwelled_radiance; at_height[AHP_DOWNWELLED_RADIANCE] = point.elevations[below].downwelled_radiance; } else { /* Interpolate between the heights for each parameter */ interp_ratio = (interpolate_height - point.elevations[above].elevation) / (point.elevations[above].elevation - point.elevations[below].elevation); below_parameters[AHP_TRANSMISSION] = point.elevations[below].transmission; below_parameters[AHP_UPWELLED_RADIANCE] = point.elevations[below].upwelled_radiance; below_parameters[AHP_DOWNWELLED_RADIANCE] = point.elevations[below].downwelled_radiance; above_parameters[AHP_TRANSMISSION] = point.elevations[above].transmission; above_parameters[AHP_UPWELLED_RADIANCE] = point.elevations[above].upwelled_radiance; above_parameters[AHP_DOWNWELLED_RADIANCE] = point.elevations[above].downwelled_radiance; for (parameter = 0; parameter < AHP_NUM_PARAMETERS; parameter++) { at_height[parameter] = interp_ratio * (above_parameters[parameter] - below_parameters[parameter]) + above_parameters[parameter]; } } } /* vertex loop */ /* Interpolate parameters at appropriate height to location of the current pixel. */ /* Shepard's method */ for (vertex = 0; vertex < NUM_CELL_POINTS; vertex++) { /* It is possible a vertex point did not have a model run. * If so, its values should not be used; set its weight to 0 */ if (!model_points[cell_vertices[vertex]].ran_model) w[vertex] = 0.0; else { double delta_x = points->points[cell_vertices[vertex]].map_x - interpolate_easting; double delta_y = points->points[cell_vertices[vertex]].map_y - interpolate_northing; w[vertex] = 1.0 / sqrt(delta_x*delta_x + delta_y*delta_y); total += w[vertex]; } } /* Normalize the weights for each vertex. */ inv_total = 1/total; for (vertex = 0; vertex < NUM_CELL_POINTS; vertex++) { w[vertex] *= inv_total; } /* For each parameter apply each vertex's weighted value */ for (parameter = 0; parameter < AHP_NUM_PARAMETERS; parameter++) { parameters[parameter] = 0.0; for (vertex = 0; vertex < NUM_CELL_POINTS; vertex++) { parameters[parameter] += w[vertex] * at_heights[vertex][parameter]; } } } /* Below are the interpolation routines for smooth interpolation. This is intended to remove the checkerboard pattern visible when the above interpolation routines are used. The smooth interpolation is used with RTTOV. It could also be used with MODTRAN but we currently are not going to the effort to do that. */ /****************************************************************************** METHOD dec_binsearch PURPOSE: Binary search on descending array RETURN: Index that was found ******************************************************************************/ int dec_binsearch ( double *d, /* I: Array to search in */ double v, /* I: Element to search for */ int imin, /* I: Minimum of range to search */ int imax /* I: Maximum of range to search */ ) { if (imin >= imax) return imin; else { int imid = imin + (imax - imin)/2; if (imid == imin) imid = imax; if (d[imid] < v) return dec_binsearch(d, v, imin, imid - 1); else if (d[imid] > v) return dec_binsearch(d, v, imid, imax); else return imid; } } /****************************************************************************** METHOD inc_binsearch PURPOSE: Binary search on increasing array RETURN: Index that was found ******************************************************************************/ int inc_binsearch ( double *d, /* I: Array to search in */ double v, /* I: Element to search for */ int imin, /* I: Minimum of range to search */ int imax /* I: Maximum of range to search */ ) { if (imin >= imax) return imin; else { int imid = imin + (imax - imin)/2; if (imid == imin) imid = imax; if (d[imid] > v) return inc_binsearch(d, v, imin, imid - 1); else if (d[imid] < v) return inc_binsearch(d, v, imid, imax); else return imid; } } /****************************************************************************** METHOD findIx PURPOSE: Find index of a value in an ordered (ascending or descending) array RETURN: Index that was found ******************************************************************************/ int findIx ( double qtarg, /* I: Element to search for */ double *samp, /* I: Array to search in */ int nsamps /* I: Size of search array */ ) { int lowix = 0; if (samp[nsamps-1] < samp[0]) { /* Samples are decreasing */ if (qtarg >= samp[1]) { /* Use the first sample point */ lowix = 0; } else if (qtarg <= samp[nsamps-2]) { /* Use the next-to-the-last sample point */ lowix = nsamps - 2; } else { /* Search */ lowix = dec_binsearch(samp, qtarg, 0, nsamps); } } else { /* Samples are increasing */ if (qtarg <= samp[1]) { lowix = 0; } else if (qtarg >= samp[nsamps-2]) { lowix = nsamps - 2; } else { lowix = inc_binsearch(samp, qtarg, 0, nsamps); } } return lowix; } /* LERP from "Graphics Gems IV", Academic Press, 1994 * linear interpolation from l (when a=0) to h (when a=1) * (equal to (a*h)+((1-a)*l) */ #define LERP(a,l,h) ((l)+(((h)-(l))*(a))) /****************************************************************************** METHOD dval2 PURPOSE: Two dimensional lookup (e.g.: image, reanalysis grid) RETURN: Value from the source 2D array ******************************************************************************/ double dval2 ( const double *d, /* I: Look up value in this flattened 2-D array */ int x, /* I: X axis location */ int y, /* I: Y axis location */ int ysize /* I: Size of Y axis */ ) { return d[x * ysize + y]; } /****************************************************************************** METHOD: setup_smooth_interpolate PURPOSE: Set up structures for convenient use during smooth interpolation ******************************************************************************/ void setup_smooth_interpolate ( int num_points, /* I: Number of reanalysis points */ int nsampsx, /* I: Number of reanalysis points in X direction */ int nsampsy, /* I: Number of reanalysis points in Y direction */ const MODEL_POINTS *model_results, /* I: results from model runs */ bool antimeridian_crossing, /* I: 0=doesn't cross, 1=does cross */ double *upwelled_radiance, /* O: Upwelled radiance */ double *downwelled_radiance, /* O: Downwelled radiance */ double *transmission, /* O: Transmission */ double *xsamp, /* O: X linear coordinate vector */ double *ysamp /* O: Y linear coordinate vector */ ) { double *X = NULL; /* Grid point latitudes */ double *Y = NULL; /* Grid point longitudes */ int ipt; /* Reanalysis point index */ int parameter_index; /* Reanalysis index offset for antimeridian */ const char *FUNC_NAME = "setup_smooth_interpolate"; /* Allocate memory for lat/lon structures */ X = malloc(num_points * sizeof(double)); Y = malloc(num_points * sizeof(double)); if (X == NULL || Y == NULL) ERROR_MESSAGE ("Allocating interpolation structure memory", FUNC_NAME); /* Extract lat/lon and atmospheric parameter values from model. The parameter values should be 0 where the model wasn't run. */ for (ipt = 0; ipt < num_points; ipt++) { X[ipt] = model_results->points[ipt].lat; Y[ipt] = model_results->points[ipt].lon; /* In regular scenarios the ipt index is right, but for antimeridian crossings, the longitudes need to be put in order. */ if (antimeridian_crossing) { parameter_index = nsampsy * model_results->points[ipt].row + model_results->points[ipt].col; } else { parameter_index = ipt; } if (model_results->points[ipt].ran_model) { upwelled_radiance[parameter_index] = model_results->points[ipt].elevations[0].upwelled_radiance; downwelled_radiance[parameter_index] = model_results->points[ipt].elevations[0].downwelled_radiance; transmission[parameter_index] = model_results->points[ipt].elevations[0].transmission; } else { upwelled_radiance[parameter_index] = 0; downwelled_radiance[parameter_index] = 0; transmission[parameter_index] = 0; } } /* Create linear coordinate vectors */ for (ipt = 0; ipt < nsampsx; ipt++) { xsamp[ipt] = dval2(X, ipt, 0, nsampsy); } for (ipt = 0; ipt < nsampsy; ipt++) { if (antimeridian_crossing) { ysamp[model_results->points[ipt].col] = model_results->points[ipt].lon; } else { ysamp[ipt] = dval2(Y, 0, ipt, nsampsy); } } } /****************************************************************************** METHOD: smooth_interpolate_parameters PURPOSE: Interpolate atmospheric parameters to the location of the current pixel. This is based on JPL's MODIS Surface Temperature interpolation for atmospheric bands. It is intended to remove the checkerboard pattern visible in current atmospheric parameter bands. RETURN: Interpolated value for the atmospheric parameter from surrounding reanalysis points for the pixel ******************************************************************************/ double smooth_interpolate_parameters ( double pixel_lat, /* I: Latitude of the current pixel */ double pixel_lon, /* I: Longitude of the current pixel */ const double *parameter, /* I: Upwelled radiance, downwelled radiance, or atmospheric transmittance */ const double *xsamp, /* I: Reanalysis latitudes along x axis */ const double *ysamp, /* I: Reanalysis longitudes along y axis*/ int ix, /* I: X index into reanalysis grid */ int iy, /* I: Y index into reanalysis grid */ int nsampsy /* I: Number of reanalysis points along y axis */ ) { double fx; /* X axis slope */ double fy; /* Y axis slope */ double d00; /* Parameter value for first point around the pixel */ double d10; /* Parameter value for second point around the pixel */ double d01; /* Parameter value for third point around the pixel */ double d11; /* Parameter value for fourth point around the pixel */ double dy0; /* Intermediate result for the pixel */ double dy1; /* Intermediate result for the pixel */ fx = (pixel_lat - xsamp[ix]) / (xsamp[ix + 1] - xsamp[ix]); fy = (pixel_lon - ysamp[iy]) / (ysamp[iy + 1] - ysamp[iy]); d00 = dval2(parameter, ix, iy, nsampsy); d10 = dval2(parameter, ix, iy + 1, nsampsy); d01 = dval2(parameter, ix + 1, iy, nsampsy); d11 = dval2(parameter, ix + 1, iy + 1, nsampsy); dy0 = LERP(fy, d00, d10); dy1 = LERP(fy, d01, d11); return(LERP(fx, dy0, dy1)); }