/***************************************************************************** DESCRIPTION: Generate data for the surface temperature uncertainty band. *****************************************************************************/ #include #include #include #include #include #include "const.h" #include "utilities.h" #include "input.h" #include "interpolate.h" /* Calculate radiance image given temperature image using the inverse Planck Function, modified to use K1 and K2 constants for a specific instrument. */ static double convert_radiance_to_temperature(double radiance, double k1, double k2) { return k2/log(k1/radiance + 1); } /* Calculate the uncertainty term, which is part of the surface temperature uncertainty estimation. */ static double get_parm_uncertainty ( double val, AT_HEIGHT_PARAMETERS parm, double coeffs[NUM_QA_COEFFS], /* quadratic coefficients */ double real_data_end ) { double uncertainty; /* computed uncertainty */ /* Set a boundary value (lower for transmission and upper for upwelled and downwelled radiances), beyond which a quadratic fit will be made, and below which last value of the fit line is extended as a constant. The lower bound was set as simply the smallest transmission value from the validation set that was used. Also set the polynomial (quadratic) coefficients calculated from MODTRAN simulations using MERRA. */ switch(parm) { case AHP_TRANSMISSION: if (val < real_data_end) val = real_data_end; break; case AHP_UPWELLED_RADIANCE: case AHP_DOWNWELLED_RADIANCE: if (val > real_data_end) val = real_data_end; break; default: RETURN_ERROR("Unsupported uncertainty type.", "get_parm_uncertainty", EXIT_FAILURE); } /* Calculate uncertainty values. */ uncertainty = coeffs[0] * val * val + coeffs[1] * val + coeffs[2]; return uncertainty; } /* Calculate the cross correlation term, which is part of the surface temperature uncertainty estimation. */ static double get_cross_correlation(double dlt_dtau, double dlt_dlu, double dlt_dld, double s_tau, double s_lu, double s_ld) { /* Correlation coefficients from MODTRAN simulations using MERRA. */ double corr_tau_lu = -0.9899; double corr_tau_ld = -0.9857; double corr_lu_ld = 0.9965; /* Calculate cross correlation terms. */ double corr_term[3] = {dlt_dtau*s_tau, dlt_dlu*s_lu, dlt_dld*s_ld}; double part_tau_lu = 2*corr_tau_lu*corr_term[0]*corr_term[1]; double part_tau_ld = 2*corr_tau_ld*corr_term[0]*corr_term[2]; double part_lu_ld = 2*corr_lu_ld*corr_term[1]*corr_term[2]; /* Calculate cross correlation. */ double cross_correlation = part_tau_lu + part_tau_ld + part_lu_ld; return cross_correlation; } /* Display help/usage information. */ static void usage() { printf("Surface Temperature - calculate_unc\n"); printf("\n"); printf("Generate the surface temperature uncertainty (QA) band.\n"); printf("\n"); printf("The arguments are listed below (in order).\nAll are required " "unless identified as optional:\n"); printf(" satellite ID (string)\n"); printf(" number of lines in band (integer)\n"); printf(" number of samples in band (integer)\n"); printf(" observed radiance filename (string)\n"); printf(" transmission filename (string)\n"); printf(" upwelled filename (string)\n"); printf(" downwelled filename (string)\n"); printf(" emissivity filename (string)\n"); printf(" emissivity standard deviation filename (string)\n"); printf(" thermal input filename (string)\n"); printf(" thermal input 2 (ETM) filename (string)\n"); printf(" uncertainty (output) filename (string)\n"); printf(" K1 thermal constant (float)\n"); printf(" K2 thermal constant (float)\n"); printf(" transmittance coefficient 1 (float)\n"); printf(" transmittance coefficient 2 (float)\n"); printf(" transmittance coefficient 3 (float)\n"); printf(" transmittance lower bound (float)\n"); printf(" upwelled radiance coefficient 1 (float)\n"); printf(" upwelled radiance coefficient 2 (float)\n"); printf(" upwelled radiance coefficient 3 (float)\n"); printf(" upwelled radiance upper bound (float)\n"); printf(" downwelled radiance coefficient 1 (float)\n"); printf(" downwelled radiance coefficient 2 (float)\n"); printf(" downwelled radiance coefficient 3 (float)\n"); printf(" downwelled radiance upper bound (float)\n"); printf(" no data fill value (integer)\n"); printf(" integer scaling conversion factor (float)\n"); printf(" number of threads to use (integer)(optional)\n"); printf("\n"); printf ("Note: This application must run from the directory" " where the input data\n is located.\n\n"); } int main(int argc, char *argv[]) { const char *FUNC_NAME = "CALCULATE_UNC"; const char *satellite;/* satellite string */ const char *uncertainty_filename; /* output filename for temperature uncertainty band */ FILE *uncertainty_fp; /* output temperature uncertainty file */ int num_args = 29; /* number of input arguments */ char msg[PATH_MAX]; /* status message */ /* data buffers */ float *observed_radiance, /* observed radiance band */ *tau, /* atmospheric transmittance radiance layer */ *upwelled, /* upwelled radiance layer */ *downwelled, /* downwelled radiance layer */ *emis, /* emissivity layer */ *emis_stdev; /* emissivity standard deviation layer */ uint8_t *thermal_1 = NULL;/* level 1 thermal band (TM/ETM) */ uint8_t *thermal_2 = NULL;/* level 1 thermal band 2 (ETM) */ uint16_t *thermal_1_oli = NULL; /* level 1 thermal band (OLI) */ int16_t *st_uncertainty; /* surface temperature uncertainty (QA) band */ int16_t fill; /* fill value */ double k1, k2; /* thermal constants */ double tau_coeffs[NUM_QA_COEFFS]; /* Transmittance coefficients */ double tau_lower_bound; /* Transmittance boundary */ double upwelled_coeffs[NUM_QA_COEFFS]; /* Upwelled radiance coefficients */ double upwelled_upper_bound; /* Upwelled radiance boundary */ double downwelled_coeffs[NUM_QA_COEFFS]; /* Downwelled radiance coefficients */ double downwelled_upper_bound; /* Downwelled radiance boundary */ double landsat_uncertainty; /* satellite thermal band uncertainty (radiance) */ double emis_regfit, /* RMSE of the linear regression fit of the spectral emissivity adjustment procedure */ emis_regfit_sqr; /* squared value used in pixel calculations */ const double eret13 = 0.0164, /* components of Temperature Emissivity Separation uncertainty */ eret14 = 0.0174; const double eret_sqr = (eret13*eret13 + eret14*eret14)/2; /* total algorithmic uncertainty of the Temperature Emissivity Separation algorithm used to produce the ASTER GED emissivities described in Hulley et al. 2012 (Actually, this is the squared value, since we only need the squared value.) */ double mult_factor; /* multiplicative factor for converting float to int */ int i; /* buffer index */ int coeff_num; /* Nuumber of the coefficient being processed */ int nlines, nsamps; /* number of lines and samples in bands */ int pixel_count; /* number of pixels in bands */ int saturation; /* Value indicating input is saturated */ #ifdef _OPENMP int nthreads; /* number of threads to use if using OpenMP */ LOG_MESSAGE("Application was compiled using multi-threading via " "-fopenmp", FUNC_NAME); #else LOG_MESSAGE("Application was not compiled using multi-threading", FUNC_NAME); #endif /* Read the command-line arguments. The last argument (number of threads) is optional. See the usage() function for a description of the arguments. */ if (argc < num_args || argc > num_args + 1) { usage(); RETURN_ERROR("Wrong number of arguments.\n" "This app is for internal use only and should not be " "run from the command line.", FUNC_NAME, EXIT_FAILURE); } i = 1; satellite = argv[i++]; if (sscanf(argv[i++], "%d", &nlines) != 1) RETURN_ERROR("Unable to retrieve number of lines from command line.", FUNC_NAME, EXIT_FAILURE); if (sscanf(argv[i++], "%d", &nsamps) != 1) RETURN_ERROR("Unable to retrieve number of samples from command line.", FUNC_NAME, EXIT_FAILURE); pixel_count = nlines*nsamps; observed_radiance = read_unc_input_float(argv[i++], pixel_count); tau = read_unc_input_float(argv[i++], pixel_count); upwelled = read_unc_input_float(argv[i++], pixel_count); downwelled = read_unc_input_float(argv[i++], pixel_count); emis = read_unc_input_float(argv[i++], pixel_count); emis_stdev = read_unc_input_float(argv[i++], pixel_count); if ((strcmp(satellite, "LANDSAT_4") == 0) || (strcmp(satellite, "LANDSAT_5") == 0) || (strcmp(satellite, "LANDSAT_7") == 0)) { thermal_1 = read_unc_input_uint8(argv[i++], pixel_count); } if ((strcmp(satellite, "LANDSAT_8") == 0) || (strcmp(satellite, "LANDSAT_9") == 0)) { thermal_1_oli = read_unc_input_uint16(argv[i++], pixel_count); } if (strcmp(satellite, "LANDSAT_7") == 0) { thermal_2 = read_unc_input_uint8(argv[i++], pixel_count); } else { /* Skip this band since it only exists for Landsat 7 */ thermal_2 = NULL; i++; } uncertainty_filename = argv[i++]; st_uncertainty = malloc(pixel_count*sizeof(int16_t)); if (observed_radiance == NULL || tau == NULL || upwelled == NULL || downwelled == NULL || emis == NULL || emis_stdev == NULL || ((strcmp(satellite, "LANDSAT_4") == 0) && (thermal_1 == NULL)) || ((strcmp(satellite, "LANDSAT_5") == 0) && (thermal_1 == NULL)) || ((strcmp(satellite, "LANDSAT_7") == 0) && (thermal_1 == NULL)) || ((strcmp(satellite, "LANDSAT_8") == 0) && (thermal_1_oli == NULL)) || ((strcmp(satellite, "LANDSAT_9") == 0) && (thermal_1_oli == NULL)) || ((strcmp(satellite, "LANDSAT_7") == 0) && (thermal_2 == NULL)) || st_uncertainty == NULL) { free(observed_radiance); free(tau); free(upwelled); free(downwelled); free(emis); free(emis_stdev); if ((strcmp(satellite, "LANDSAT_4") == 0) || (strcmp(satellite, "LANDSAT_5") == 0) || (strcmp(satellite, "LANDSAT_7") == 0)) { free(thermal_1); } if ((strcmp(satellite, "LANDSAT_8") == 0) || (strcmp(satellite, "LANDSAT_9") == 0)) { free(thermal_1_oli); } if (strcmp(satellite, "LANDSAT_7") == 0) { free(thermal_2); } free(st_uncertainty); RETURN_ERROR("Unable to read all input files", FUNC_NAME, EXIT_FAILURE); } if (sscanf(argv[i++], "%lf", &k1) != 1) RETURN_ERROR("Unable to retrieve thermal constant K1 from " "command line.", FUNC_NAME, EXIT_FAILURE); if (sscanf(argv[i++], "%lf", &k2) != 1) RETURN_ERROR("Unable to retrieve thermal constant K2 from " "command line.", FUNC_NAME, EXIT_FAILURE); for (coeff_num = 0; coeff_num < NUM_QA_COEFFS; coeff_num++) { if (sscanf(argv[i++], "%lf", &tau_coeffs[coeff_num]) != 1) RETURN_ERROR("Unable to retrieve trasmittance coefficient from " "command line.", FUNC_NAME, EXIT_FAILURE); } if (sscanf(argv[i++], "%lf", &tau_lower_bound) != 1) RETURN_ERROR("Unable to retrieve transmittance lower bound from " "command line.", FUNC_NAME, EXIT_FAILURE); for (coeff_num = 0; coeff_num < NUM_QA_COEFFS; coeff_num++) { if (sscanf(argv[i++], "%lf", &upwelled_coeffs[coeff_num]) != 1) RETURN_ERROR("Unable to retrieve upwelled radiance coefficient " "from command line.", FUNC_NAME, EXIT_FAILURE); } if (sscanf(argv[i++], "%lf", &upwelled_upper_bound) != 1) RETURN_ERROR("Unable to retrieve upwelled radiance upper bound from " "command line.", FUNC_NAME, EXIT_FAILURE); for (coeff_num = 0; coeff_num < NUM_QA_COEFFS; coeff_num++) { if (sscanf(argv[i++], "%lf", &downwelled_coeffs[coeff_num]) != 1) RETURN_ERROR("Unable to retrieve downwelled radiance coefficient " "from command line.", FUNC_NAME, EXIT_FAILURE); } if (sscanf(argv[i++], "%lf", &downwelled_upper_bound) != 1) RETURN_ERROR("Unable to retrieve downwelled radiance upper bound from " "command line.", FUNC_NAME, EXIT_FAILURE); if (sscanf(argv[i++], "%hd", &fill) != 1) RETURN_ERROR("Unable to retrieve \"no data\" fill value from " "command line.", FUNC_NAME, EXIT_FAILURE); if (sscanf(argv[i++], "%lf", &mult_factor) != 1) RETURN_ERROR("Unable to retrieve conversion scaling factor from " "command line.", FUNC_NAME, EXIT_FAILURE); #ifdef _OPENMP if (argc == num_args) nthreads = 1; else if (sscanf(argv[i], "%d", &nthreads) != 1) RETURN_ERROR("Unable to retrieve the number of threads from " "command line.", FUNC_NAME, EXIT_FAILURE); #endif /* Look up satellite thermal band uncertainty in radiance based on satellite. These are based on the following uncertainty values in K: Landsat 4: 0.8 Landsat 5: 0.8 Landsat 7: 0.4 Landsat 8: 0.3 Landsat 9: 0.1 Get the RMSE of the linear regression fit of the spectral emissivity adjustment procedure (which is the emis_data calculation in estimate_landsat_emissivity). */ if (strcmp(satellite, "LANDSAT_4") == 0) { landsat_uncertainty = 0.08959421; emis_regfit = 0.00085135; } else if (strcmp(satellite, "LANDSAT_5") == 0) { landsat_uncertainty = 0.08959421; emis_regfit = 0.0013; } else if (strcmp(satellite, "LANDSAT_7") == 0) { landsat_uncertainty = 0.04471454; emis_regfit = 0.0011; } else if (strcmp(satellite, "LANDSAT_8") == 0) { landsat_uncertainty = 0.03577072; emis_regfit = 0.00093909; } else if (strcmp(satellite, "LANDSAT_9") == 0) { landsat_uncertainty = 0.0144; /* From RIT using a 300K LUT */ emis_regfit = 0.00084125; } else { free(observed_radiance); free(tau); free(upwelled); free(downwelled); free(emis); free(emis_stdev); if ((strcmp(satellite, "LANDSAT_4") == 0) || (strcmp(satellite, "LANDSAT_5") == 0) || (strcmp(satellite, "LANDSAT_7") == 0)) { free(thermal_1); } if ((strcmp(satellite, "LANDSAT_8") == 0) || (strcmp(satellite, "LANDSAT_9") == 0)) { free(thermal_1_oli); } if (strcmp(satellite, "LANDSAT_7") == 0) { free(thermal_2); } free(st_uncertainty); RETURN_ERROR("Unsupported satellite specified.", FUNC_NAME, EXIT_FAILURE); } emis_regfit_sqr = emis_regfit * emis_regfit; /* square value for later use*/ /* Define the input saturation value for the satellite */ if ((strcmp(satellite, "LANDSAT_4") == 0) || (strcmp(satellite, "LANDSAT_5") == 0) || (strcmp(satellite, "LANDSAT_7") == 0)) { saturation = ST_UNC_L1_THERMAL_SATURATION; } if ((strcmp(satellite, "LANDSAT_8") == 0) || (strcmp(satellite, "LANDSAT_9") == 0)) { saturation = ST_UNC_L1_THERMAL_SATURATION_OLI; } #ifdef _OPENMP #pragma omp parallel for private (i) \ num_threads(nthreads) #endif for (i = 0; i < pixel_count; i++) { double le_image; /* radiance image value */ double dlt_dtau, dlt_dlu, dlt_dld, dlt_dlobs, dlt_demis; /* partial derivatives for Jacobian of radiance */ double s_tau, s_lu, s_ld; /* transmission, upwelled, and downwelled radiance uncertainty */ double s_a; /* atmospheric uncertainty */ double s_i; /* instrument uncertainty */ double s_e; /* emissivity uncertainty */ double s_p; /* cross correlation */ double st_uncertainty_radiance; /* uncertainty in radiance units */ double delta_radiance; /* radiance difference from nominal */ double radiance_minus_delta, /* delta offsets from nominal radiance */ radiance_plus_delta; double temp_uncertainty_high, /* high and low temperature values */ temp_uncertainty_low; double temp; /* temporary (intermediate) value */ /* Skip fill values. */ if (observed_radiance[i] == fill || emis[i] == fill || emis_stdev[i] == fill) { st_uncertainty[i] = fill; continue; } /* Skip values that are saturated in the input. */ if ((strcmp(satellite, "LANDSAT_4") == 0) || (strcmp(satellite, "LANDSAT_5") == 0)) { if (thermal_1[i] == saturation) { st_uncertainty[i] = ST_NO_DATA_VALUE; continue; } } else if ((strcmp(satellite, "LANDSAT_8") == 0) || (strcmp(satellite, "LANDSAT_9") == 0)) { if (thermal_1_oli[i] == saturation) { st_uncertainty[i] = ST_NO_DATA_VALUE; continue; } } else /* satellite is "LANDSAT_7" */ { if ((thermal_1[i] == saturation) && (thermal_2[i] == saturation)) { st_uncertainty[i] = ST_NO_DATA_VALUE; continue; } } /* Compute the radiance that will be used as nominal value around which the uncertainty value will operate. */ temp = emis[i]*tau[i]; le_image = (observed_radiance[i] - upwelled[i] - downwelled[i]*(tau[i] - temp))/temp; if (le_image > 30) le_image = 0; /* Calculate the partial derviatives of the radiance. */ dlt_dtau = (upwelled[i] - observed_radiance[i])/(temp*tau[i]); dlt_dlu = -1/temp; dlt_dld = (emis[i] - 1)/emis[i]; dlt_dlobs = -dlt_dlu; dlt_demis = ((upwelled[i] - observed_radiance[i] + downwelled[i] * tau[i]) / (temp * emis[i])); /* Calculate transmission, upwelled, and downwelled radiance uncertainty. */ s_tau = get_parm_uncertainty(tau[i], AHP_TRANSMISSION, tau_coeffs, tau_lower_bound); s_lu = get_parm_uncertainty(upwelled[i], AHP_UPWELLED_RADIANCE, upwelled_coeffs, upwelled_upper_bound); s_ld = get_parm_uncertainty(downwelled[i], AHP_DOWNWELLED_RADIANCE, downwelled_coeffs, downwelled_upper_bound); /* Calculate the atmospheric uncertainty. */ s_a = (dlt_dtau*s_tau)*(dlt_dtau*s_tau) + (dlt_dlu*s_lu)*(dlt_dlu*s_lu) + (dlt_dld*s_ld)*(dlt_dld*s_ld); /* Calculate instrument uncertainty. */ s_i = (dlt_dlobs*landsat_uncertainty)*(dlt_dlobs*landsat_uncertainty); /* Calculate the emissivity uncertainty. */ s_e = dlt_demis * sqrt(emis_stdev[i] * emis_stdev[i] + emis_regfit_sqr + eret_sqr); s_e = s_e * s_e; /* Calculate the cross correlation. */ s_p = get_cross_correlation(dlt_dtau, dlt_dlu, dlt_dld, s_tau, s_lu, s_ld); st_uncertainty_radiance = sqrt(s_a + s_i + s_e + s_p); /* Again, the st_uncertainty value is a residual about some nominal radiance value, so we must convert to residual temperature about some nominal reference radiance, which we get from the original radiance image. */ delta_radiance = 0.5*st_uncertainty_radiance; radiance_minus_delta = le_image - delta_radiance; radiance_plus_delta = le_image + delta_radiance; /* Convert from radiance to temperature (K). */ temp_uncertainty_high = convert_radiance_to_temperature(radiance_plus_delta, k1, k2); temp_uncertainty_low = convert_radiance_to_temperature(radiance_minus_delta, k1, k2); if (isnan(temp_uncertainty_high) || isnan(temp_uncertainty_low)) { /* Set ST uncertainty to invalid pixel since it can't be calculated. */ st_uncertainty[i] = ST_UNC_INVALID_PIXEL; } else { /* The total uncertainty is the difference between the high and low values, scaled by the multiplication factor. */ st_uncertainty[i] = roundf((temp_uncertainty_high - temp_uncertainty_low) * mult_factor); } if (st_uncertainty[i] > 100*mult_factor) st_uncertainty[i] = ST_UNC_INVALID_PIXEL; } /* pixel loop */ /* Free memory. */ free(observed_radiance); free(tau); free(upwelled); free(downwelled); free(emis); free(emis_stdev); if ((strcmp(satellite, "LANDSAT_4") == 0) || (strcmp(satellite, "LANDSAT_5") == 0) || (strcmp(satellite, "LANDSAT_7") == 0)) { free(thermal_1); } if ((strcmp(satellite, "LANDSAT_8") == 0) || (strcmp(satellite, "LANDSAT_9") == 0)) { free(thermal_1_oli); } if (strcmp(satellite, "LANDSAT_7") == 0) { free(thermal_2); } /* Write uncertainty data to disk. */ uncertainty_fp = fopen(uncertainty_filename, "wb"); if (uncertainty_fp == NULL) { free(st_uncertainty); sprintf(msg, "Opening intermediate file: %s", uncertainty_filename); RETURN_ERROR(msg, FUNC_NAME, EXIT_FAILURE); } if (fwrite(st_uncertainty, sizeof(int16_t), pixel_count, uncertainty_fp) != pixel_count) { free(st_uncertainty); sprintf(msg, "Writing to %s", uncertainty_filename); RETURN_ERROR(msg, FUNC_NAME, EXIT_FAILURE); } free(st_uncertainty); if (fclose(uncertainty_fp)) { sprintf(msg, "Closing file %s", uncertainty_filename); RETURN_ERROR(msg, FUNC_NAME, EXIT_FAILURE); } return EXIT_SUCCESS; }