[GRASS-SVN] r32986 - in grass-addons/gipe: i.dn2full.l5 i.dn2full.l7
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Aug 21 19:12:57 EDT 2008
Author: ychemin
Date: 2008-08-21 19:12:57 -0400 (Thu, 21 Aug 2008)
New Revision: 32986
Modified:
grass-addons/gipe/i.dn2full.l5/date2doy.c
grass-addons/gipe/i.dn2full.l5/dn2rad_landsat5.c
grass-addons/gipe/i.dn2full.l5/l5inread.c
grass-addons/gipe/i.dn2full.l5/main.c
grass-addons/gipe/i.dn2full.l5/rad2ref_landsat5.c
grass-addons/gipe/i.dn2full.l5/tempk_landsat5.c
grass-addons/gipe/i.dn2full.l7/date2doy.c
grass-addons/gipe/i.dn2full.l7/main.c
grass-addons/gipe/i.dn2full.l7/rad2ref_landsat7.c
grass-addons/gipe/i.dn2full.l7/tempk_landsat7.c
Log:
GRASS GIS standard indentation parsing
Modified: grass-addons/gipe/i.dn2full.l5/date2doy.c
===================================================================
--- grass-addons/gipe/i.dn2full.l5/date2doy.c 2008-08-21 23:06:30 UTC (rev 32985)
+++ grass-addons/gipe/i.dn2full.l5/date2doy.c 2008-08-21 23:12:57 UTC (rev 32986)
@@ -1,66 +1,70 @@
+
/*********************************************/
-/*This program converts day/month/year to doy*/
+/*This program converts day/month/year to doy */
+
/*********************************************/
+
/*********************************************/
int date2doy(int day, int month, int year)
{
- int leap=0;
- int day_month_tot=0;
- int doy;
+ int leap = 0;
- doy=0;
+ int day_month_tot = 0;
-/*printf("Date is %i/%i/%i\n", day, month, year);*/
+ int doy;
- if (month == 1) {
- day_month_tot = 0;
- }
- else if (month == 2) {
- day_month_tot = 31;
- }
- else if (month == 3) {
- day_month_tot = 59;
- }
- else if (month == 4) {
- day_month_tot = 90;
- }
- else if (month == 5) {
- day_month_tot = 120;
- }
- else if (month == 6) {
- day_month_tot = 151;
- }
- else if (month == 7) {
- day_month_tot = 181;
- }
- else if (month == 8) {
- day_month_tot = 212;
- }
- else if (month == 9) {
- day_month_tot = 243;
- }
- else if (month == 10) {
- day_month_tot = 273;
- }
- else if (month == 11) {
- day_month_tot = 304;
- }
- else if (month == 12) {
- day_month_tot = 334;
- }
-
- /* Leap year if dividing by 4 leads % 0.0*/
-
- if (year/4*4 == year) {
- leap = 1;
- }
-
- doy=day_month_tot+day;
- if(doy>59){
- doy=day_month_tot+day+leap;
- }
+ doy = 0;
- return(doy);
+ /*printf("Date is %i/%i/%i\n", day, month, year); */
+
+ if (month == 1) {
+ day_month_tot = 0;
+ }
+ else if (month == 2) {
+ day_month_tot = 31;
+ }
+ else if (month == 3) {
+ day_month_tot = 59;
+ }
+ else if (month == 4) {
+ day_month_tot = 90;
+ }
+ else if (month == 5) {
+ day_month_tot = 120;
+ }
+ else if (month == 6) {
+ day_month_tot = 151;
+ }
+ else if (month == 7) {
+ day_month_tot = 181;
+ }
+ else if (month == 8) {
+ day_month_tot = 212;
+ }
+ else if (month == 9) {
+ day_month_tot = 243;
+ }
+ else if (month == 10) {
+ day_month_tot = 273;
+ }
+ else if (month == 11) {
+ day_month_tot = 304;
+ }
+ else if (month == 12) {
+ day_month_tot = 334;
+ }
+
+ /* Leap year if dividing by 4 leads % 0.0 */
+
+ if (year / 4 * 4 == year) {
+ leap = 1;
+ }
+
+ doy = day_month_tot + day;
+ if (doy > 59) {
+ doy = day_month_tot + day + leap;
+ }
+
+ return (doy);
}
-
Modified: grass-addons/gipe/i.dn2full.l5/dn2rad_landsat5.c
===================================================================
--- grass-addons/gipe/i.dn2full.l5/dn2rad_landsat5.c 2008-08-21 23:06:30 UTC (rev 32985)
+++ grass-addons/gipe/i.dn2full.l5/dn2rad_landsat5.c 2008-08-21 23:12:57 UTC (rev 32986)
@@ -1,162 +1,170 @@
-#include<stdio.h>
-#include<math.h>
-#include<stdlib.h>
-
-/* Conversion of DN to Radiance for Landsat 5TM */
-
-double dn2rad_landsat5(int c_year,int c_month,int c_day,int year,int month,int day,int band,int DN)
-{
- double result, gain, bias;
- int gain_mode;
-
- if(c_year<2003){
- gain_mode=1;
- } else if (c_year==2003){
- if(c_month<5){
- gain_mode=1;
- } else if(c_month==5){
- if(c_day<5){
- gain_mode=1;
- }else{
- gain_mode=2;
- }
- } else {
- gain_mode=2;
- }
- } else if (c_year>2003&&c_year<2007){
- gain_mode=2;
- } else if (c_year==2007){
- if(c_month<5){
- gain_mode=2;
- } else if(year<1992){
- gain_mode=3;
- } else {
- gain_mode=4;
- }
- }
- if (gain_mode==1){
- if (band==1){
- gain = 0.602431;
- bias = -1.52;
- }
- if (band==2){
- gain = 1.175100;
- bias = -2.84;
- }
- if (band==3){
- gain = 0.805765;
- bias = -1.17;
- }
- if (band==4){
- gain = 0.814549;
- bias = -1.51;
- }
- if (band==5){
- gain = 0.108078;
- bias = -0.37;
- }
- if (band==6){
- gain = 0.055158;
- bias = 1.2378;
- }
- if (band==7){
- gain = 0.056980;
- bias = -0.15;
- }
- }
- if (gain_mode==2){
- if (band==1){
- gain = 0.762824;
- bias = -1.52;
- }
- if (band==2){
- gain = 1.442510;
- bias = -2.84;
- }
- if (band==3){
- gain = 1.039880;
- bias = -1.17;
- }
- if (band==4){
- gain = 0.872588;
- bias = -1.51;
- }
- if (band==5){
- gain = 0.119882;
- bias = -0.37;
- }
- if (band==6){
- gain = 0.055158;
- bias = 1.2378;
- }
- if (band==7){
- gain = 0.065294;
- bias = -0.15;
- }
- }
- if (gain_mode==3){
- if (band==1){
- gain = 0.668706;
- bias = -1.52;
- }
- if (band==2){
- gain = 1.317020;
- bias = -2.84;
- }
- if (band==3){
- gain = 1.039880;
- bias = -1.17;
- }
- if (band==4){
- gain = 0.872588;
- bias = -1.51;
- }
- if (band==5){
- gain = 0.119882;
- bias = -0.37;
- }
- if (band==6){
- gain = 0.055158;
- bias = 1.2378;
- }
- if (band==7){
- gain = 0.065294;
- bias = -0.15;
- }
- }
- if (gain_mode==4){
- if (band==1){
- gain = 0.762824;
- bias = -1.52;
- }
- if (band==2){
- gain = 1.442510;
- bias = -2.84;
- }
- if (band==3){
- gain = 1.039880;
- bias = -1.17;
- }
- if (band==4){
- gain = 0.872588;
- bias = -1.51;
- }
- if (band==5){
- gain = 0.119882;
- bias = -0.37;
- }
- if (band==6){
- gain = 0.055158;
- bias = 1.2378;
- }
- if (band==7){
- gain = 0.065294;
- bias = -0.15;
- }
- }
-
- result = gain * (double) DN + bias ;
-
- return result;
-}
-
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+
+ /* Conversion of DN to Radiance for Landsat 5TM */
+double dn2rad_landsat5(int c_year, int c_month, int c_day, int year,
+ int month, int day, int band, int DN)
+{
+ double result, gain, bias;
+
+ int gain_mode;
+
+ if (c_year < 2003) {
+ gain_mode = 1;
+ }
+ else if (c_year == 2003) {
+ if (c_month < 5) {
+ gain_mode = 1;
+ }
+ else if (c_month == 5) {
+ if (c_day < 5) {
+ gain_mode = 1;
+ }
+ else {
+ gain_mode = 2;
+ }
+ }
+ else {
+ gain_mode = 2;
+ }
+ }
+ else if (c_year > 2003 && c_year < 2007) {
+ gain_mode = 2;
+ }
+ else if (c_year == 2007) {
+ if (c_month < 5) {
+ gain_mode = 2;
+ }
+ else if (year < 1992) {
+ gain_mode = 3;
+ }
+ else {
+ gain_mode = 4;
+ }
+ }
+ if (gain_mode == 1) {
+ if (band == 1) {
+ gain = 0.602431;
+ bias = -1.52;
+ }
+ if (band == 2) {
+ gain = 1.175100;
+ bias = -2.84;
+ }
+ if (band == 3) {
+ gain = 0.805765;
+ bias = -1.17;
+ }
+ if (band == 4) {
+ gain = 0.814549;
+ bias = -1.51;
+ }
+ if (band == 5) {
+ gain = 0.108078;
+ bias = -0.37;
+ }
+ if (band == 6) {
+ gain = 0.055158;
+ bias = 1.2378;
+ }
+ if (band == 7) {
+ gain = 0.056980;
+ bias = -0.15;
+ }
+ }
+ if (gain_mode == 2) {
+ if (band == 1) {
+ gain = 0.762824;
+ bias = -1.52;
+ }
+ if (band == 2) {
+ gain = 1.442510;
+ bias = -2.84;
+ }
+ if (band == 3) {
+ gain = 1.039880;
+ bias = -1.17;
+ }
+ if (band == 4) {
+ gain = 0.872588;
+ bias = -1.51;
+ }
+ if (band == 5) {
+ gain = 0.119882;
+ bias = -0.37;
+ }
+ if (band == 6) {
+ gain = 0.055158;
+ bias = 1.2378;
+ }
+ if (band == 7) {
+ gain = 0.065294;
+ bias = -0.15;
+ }
+ }
+ if (gain_mode == 3) {
+ if (band == 1) {
+ gain = 0.668706;
+ bias = -1.52;
+ }
+ if (band == 2) {
+ gain = 1.317020;
+ bias = -2.84;
+ }
+ if (band == 3) {
+ gain = 1.039880;
+ bias = -1.17;
+ }
+ if (band == 4) {
+ gain = 0.872588;
+ bias = -1.51;
+ }
+ if (band == 5) {
+ gain = 0.119882;
+ bias = -0.37;
+ }
+ if (band == 6) {
+ gain = 0.055158;
+ bias = 1.2378;
+ }
+ if (band == 7) {
+ gain = 0.065294;
+ bias = -0.15;
+ }
+ }
+ if (gain_mode == 4) {
+ if (band == 1) {
+ gain = 0.762824;
+ bias = -1.52;
+ }
+ if (band == 2) {
+ gain = 1.442510;
+ bias = -2.84;
+ }
+ if (band == 3) {
+ gain = 1.039880;
+ bias = -1.17;
+ }
+ if (band == 4) {
+ gain = 0.872588;
+ bias = -1.51;
+ }
+ if (band == 5) {
+ gain = 0.119882;
+ bias = -0.37;
+ }
+ if (band == 6) {
+ gain = 0.055158;
+ bias = 1.2378;
+ }
+ if (band == 7) {
+ gain = 0.065294;
+ bias = -0.15;
+ }
+ }
+ result = gain * (double)DN + bias;
+ return result;
+}
+
+
Modified: grass-addons/gipe/i.dn2full.l5/l5inread.c
===================================================================
--- grass-addons/gipe/i.dn2full.l5/l5inread.c 2008-08-21 23:06:30 UTC (rev 32985)
+++ grass-addons/gipe/i.dn2full.l5/l5inread.c 2008-08-21 23:12:57 UTC (rev 32986)
@@ -4,89 +4,92 @@
#include <grass/gis.h>
#include <grass/glocale.h>
-int l5_in_read(char *metfName,int *path, int *row, double *latitude,double *longitude,double *sun_elevation, double *sun_azimuth,int *c_year, int *c_month, int *c_day,int *day, int *month, int *year,double *decimal_hour)
+int l5_in_read(char *metfName, int *path, int *row, double *latitude,
+ double *longitude, double *sun_elevation, double *sun_azimuth,
+ int *c_year, int *c_month, int *c_day, int *day, int *month,
+ int *year, double *decimal_hour)
{
- FILE *f;
- char s[1000], *ptr;
- int i=0;
- char *p;
- int hour, minute;
- float second;
+ FILE *f;
- if((f=fopen(metfName,"r"))==NULL){
- G_fatal_error (_("NLAPS report file [%s] not found, are you in the proper directory?"), metfName);
- }
+ char s[1000], *ptr;
- while (fgets(s,1000,f)!=NULL)
- {
- ptr = strstr(s, "Strip no.");
- if (ptr != NULL)
- {
- p = strtok(ptr, ":");
- p = strtok(NULL, " ");
- *path = atoi(p);
- p = strtok(NULL, " ");
- p = strtok(NULL, "Start Row no.:");
- *row = atoi(p);
- }
- ptr = strstr(s, "center lat");
- if (ptr != NULL)
- {
- p = strtok(ptr, ":");
- p = strtok(NULL, " ");
- *latitude = atof(p);
- p = strtok(NULL, " deg");
- p = strtok(NULL, " Scene center long:");
- *longitude = atof(p);
- }
- ptr = strstr(s, "Sun Elevation");
- if (ptr != NULL)
- {
- p = strtok(ptr, ":");
- p = strtok(NULL, " ");
- *sun_elevation = atof(p);
- p = strtok(NULL, " deg");
- p = strtok(NULL, "Sun Azimuth:");
- *sun_azimuth = atof(p);
- }
- ptr = strstr(s, "center date");
- if (ptr != NULL)
- {
- p = strtok(ptr, ":");
- p = strtok(NULL, " ");
- *year = atoi(p);
- p = strtok(NULL, " ");
- *month = atoi(p);
- p = strtok(NULL, " ");
- *day = atoi(p);
- p = strtok(NULL, " Scene center time:");
- hour = atoi(p);
- p = strtok(NULL, ":");
- minute = atoi(p);
- p = strtok(NULL, ":");
- second = atof(p);
- }
- ptr = strstr(s, "Completion date");
- if (ptr != NULL)
- {
- p = strtok(ptr, ":");
- p = strtok(NULL, " ");
- *c_year = atoi(p);
- p = strtok(NULL, " ");
- *c_month = atoi(p);
- p = strtok(NULL, " ");
- *c_day = atoi(p);
- }
+ int i = 0;
+
+ char *p;
+
+ int hour, minute;
+
+ float second;
+
+ if ((f = fopen(metfName, "r")) == NULL) {
+ G_fatal_error(_("NLAPS report file [%s] not found, are you in the proper directory?"),
+ metfName);
+ }
+
+ while (fgets(s, 1000, f) != NULL) {
+ ptr = strstr(s, "Strip no.");
+ if (ptr != NULL) {
+ p = strtok(ptr, ":");
+ p = strtok(NULL, " ");
+ *path = atoi(p);
+ p = strtok(NULL, " ");
+ p = strtok(NULL, "Start Row no.:");
+ *row = atoi(p);
}
+ ptr = strstr(s, "center lat");
+ if (ptr != NULL) {
+ p = strtok(ptr, ":");
+ p = strtok(NULL, " ");
+ *latitude = atof(p);
+ p = strtok(NULL, " deg");
+ p = strtok(NULL, " Scene center long:");
+ *longitude = atof(p);
+ }
+ ptr = strstr(s, "Sun Elevation");
+ if (ptr != NULL) {
+ p = strtok(ptr, ":");
+ p = strtok(NULL, " ");
+ *sun_elevation = atof(p);
+ p = strtok(NULL, " deg");
+ p = strtok(NULL, "Sun Azimuth:");
+ *sun_azimuth = atof(p);
+ }
+ ptr = strstr(s, "center date");
+ if (ptr != NULL) {
+ p = strtok(ptr, ":");
+ p = strtok(NULL, " ");
+ *year = atoi(p);
+ p = strtok(NULL, " ");
+ *month = atoi(p);
+ p = strtok(NULL, " ");
+ *day = atoi(p);
+ p = strtok(NULL, " Scene center time:");
+ hour = atoi(p);
+ p = strtok(NULL, ":");
+ minute = atoi(p);
+ p = strtok(NULL, ":");
+ second = atof(p);
+ }
+ ptr = strstr(s, "Completion date");
+ if (ptr != NULL) {
+ p = strtok(ptr, ":");
+ p = strtok(NULL, " ");
+ *c_year = atoi(p);
+ p = strtok(NULL, " ");
+ *c_month = atoi(p);
+ p = strtok(NULL, " ");
+ *c_day = atoi(p);
+ }
+ }
-/* printf("year = %i\n", year);
- printf("month = %i\n", month);
- printf("day = %i\n", day);
- printf("sun azimuth = %f\n", sun_azimuth);
- printf("sun elevation = %f\n", sun_elevation);
- printf("latitude = %f\n", latitude);
- printf("longitude = %f\n", longitude);
-*/
- (void)fclose(f);
- return;
+ /* printf("year = %i\n", year);
+ printf("month = %i\n", month);
+ printf("day = %i\n", day);
+ printf("sun azimuth = %f\n", sun_azimuth);
+ printf("sun elevation = %f\n", sun_elevation);
+ printf("latitude = %f\n", latitude);
+ printf("longitude = %f\n", longitude);
+ */
+ (void)fclose(f);
+ return;
}
Modified: grass-addons/gipe/i.dn2full.l5/main.c
===================================================================
--- grass-addons/gipe/i.dn2full.l5/main.c 2008-08-21 23:06:30 UTC (rev 32985)
+++ grass-addons/gipe/i.dn2full.l5/main.c 2008-08-21 23:12:57 UTC (rev 32986)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: i.dn2full.l5
@@ -22,441 +23,519 @@
#define MAXFILES 7
-//sun exo-atmospheric irradiance
+/*sun exo-atmospheric irradiance */
#define KEXO1 1969.0
#define KEXO2 1840.0
#define KEXO3 1551.0
#define KEXO4 1044.0
#define KEXO5 225.7
-#define KEXO6 00.0//filling only
+#define KEXO6 00.0 /*filling only */
#define KEXO7 82.07
#define PI 3.1415926
-int l5_in_read(char *metfName, int *path, int *row, double *latitude,double *longitude,double *sun_elevation, double *sun_azimuth,int *c_year, int *c_month, int *c_day, int *day, int *month, int *year,double *decimal_hour);
+int l5_in_read(char *metfName, int *path, int *row, double *latitude,
+ double *longitude, double *sun_elevation, double *sun_azimuth,
+ int *c_year, int *c_month, int *c_day, int *day, int *month,
+ int *year, double *decimal_hour);
int date2doy(int day, int month, int year);
-double dn2rad_landsat5(int c_year, int c_month, int c_day, int year, int month, int day, int band, int DN );
-double rad2ref_landsat5( double radiance, double doy,double sun_elevation, double k_exo );
-double tempk_landsat5( double l6 );
-int
-main(int argc, char *argv[])
+double dn2rad_landsat5(int c_year, int c_month, int c_day, int year,
+ int month, int day, int band, int DN);
+double rad2ref_landsat5(double radiance, double doy, double sun_elevation,
+ double k_exo);
+double tempk_landsat5(double l6);
+
+int main(int argc, char *argv[])
{
- struct Cell_head cellhd;/*region+header info*/
- char *mapset; /*mapset name*/
- int nrows, ncols;
- int row,col;
+ struct Cell_head cellhd; /*region+header info */
- struct GModule *module;
- struct Option *input, *input1, *output;
-
- struct Flag *flag1;
- struct History history; /*metadata*/
+ char *mapset; /*mapset name */
+ int nrows, ncols;
+
+ int row, col;
+
+ struct GModule *module;
+
+ struct Option *input, *input1, *output;
+
+ struct Flag *flag1;
+
+ struct History history; /*metadata */
+
/************************************/
- /* FMEO Declarations*****************/
- char history_buf[200];
- char *name; /*input raster name*/
- char *result; /*output raster name*/
- /*Prepare new names for output files*/
- char result1[80], result2[80], result3[80], result4[80];
- char result5[80], result6[80], result7[80];
-
- /*File Descriptors*/
- int infd[MAXFILES];
- int outfd[MAXFILES];
+ /* FMEO Declarations**************** */
+ char history_buf[200];
- char **names;
- char **ptr;
-
- int i=0,j=0;
-
- void *inrast[MAXFILES];
- DCELL *outrast[MAXFILES];
- RASTER_MAP_TYPE in_data_type[MAXFILES];
- RASTER_MAP_TYPE out_data_type = DCELL_TYPE; /* 0=numbers 1=text */
+ char *name; /*input raster name */
- double kexo[MAXFILES];
- /*Metfile*/
- char *metfName; /*NLAPS report file, header in text format*/
- double sun_elevation;
- double sun_azimuth;/*not useful here, only for parser()*/
- int c_day,c_month,c_year;/*NLAPS processing date*/
- int day,month,year;
- double decimal_hour;
- double latitude;
- double longitude;
- int l5path, l5row;
- /*EndofMetfile*/
- int doy;
- char b1[80],b2[80],b3[80];
- char b4[80],b5[80];
- char b6[80],b7[80];/*Load .tif names*/
- int temp;
+ char *result; /*output raster name */
+
+ /*Prepare new names for output files */
+ char result1[80], result2[80], result3[80], result4[80];
+
+ char result5[80], result6[80], result7[80];
+
+ /*File Descriptors */
+ int infd[MAXFILES];
+
+ int outfd[MAXFILES];
+
+ char **names;
+
+ char **ptr;
+
+ int i = 0, j = 0;
+
+ void *inrast[MAXFILES];
+
+ DCELL *outrast[MAXFILES];
+
+ RASTER_MAP_TYPE in_data_type[MAXFILES];
+
+ RASTER_MAP_TYPE out_data_type = DCELL_TYPE; /* 0=numbers 1=text */
+
+ double kexo[MAXFILES];
+
+ /*Metfile */
+ char *metfName; /*NLAPS report file, header in text format */
+
+ double sun_elevation;
+
+ double sun_azimuth; /*not useful here, only for parser() */
+
+ int c_day, c_month, c_year; /*NLAPS processing date */
+
+ int day, month, year;
+
+ double decimal_hour;
+
+ double latitude;
+
+ double longitude;
+
+ int l5path, l5row;
+
+ /*EndofMetfile */
+ int doy;
+
+ char b1[80], b2[80], b3[80];
+
+ char b4[80], b5[80];
+
+ char b6[80], b7[80]; /*Load .tif names */
+
+ int temp;
+
/************************************/
- G_gisinit(argv[0]);
+ G_gisinit(argv[0]);
- module = G_define_module();
- module->keywords = _("DN, reflectance, temperature, import");
- module->description =
- _("Calculates Top of Atmosphere Reflectance/Temperature from Landsat 5 DN.\n");
+ module = G_define_module();
+ module->keywords = _("DN, reflectance, temperature, import");
+ module->description =
+ _("Calculates Top of Atmosphere Reflectance/Temperature from Landsat 5 DN.\n");
- /* Define the different options */
- input = G_define_option() ;
- input->key = _("file");
- input->type = TYPE_STRING;
- input->required = YES;
- input->gisprompt = _("old_file,file,file");
- input->description= _("Landsat 5TM NLAPS processing report File (.txt)");
+ /* Define the different options */
+ input = G_define_option();
+ input->key = _("file");
+ input->type = TYPE_STRING;
+ input->required = YES;
+ input->gisprompt = _("old_file,file,file");
+ input->description = _("Landsat 5TM NLAPS processing report File (.txt)");
- input1 = G_define_standard_option(G_OPT_R_INPUTS) ;
- input1->key = _("input");
- input1->required = NO;
- input1->description= _("Names of all L5 DN layers (1,2,3,4,5,6,7)");
+ input1 = G_define_standard_option(G_OPT_R_INPUTS);
+ input1->key = _("input");
+ input1->required = NO;
+ input1->description = _("Names of all L5 DN layers (1,2,3,4,5,6,7)");
- output = G_define_standard_option(G_OPT_R_OUTPUT) ;
- output->key = _("output");
- output->description= _("Base name of the output layers (will add .x)");
+ output = G_define_standard_option(G_OPT_R_OUTPUT);
+ output->key = _("output");
+ output->description = _("Base name of the output layers (will add .x)");
+
/********************/
- if (G_parser(argc, argv))
- exit (-1);
-
- metfName = input->answer;
- if(input1->answers){
- names = input1->answers;
- ptr = input1->answers;
- }
- result = output->answer;
+ if (G_parser(argc, argv))
+ exit(-1);
+
+ metfName = input->answer;
+ if (input1->answers) {
+ names = input1->answers;
+ ptr = input1->answers;
+ }
+ result = output->answer;
+
/******************************************/
- /*Fetch parameters for DN2Rad2Ref correction*/
- l5_in_read(metfName,&l5path,&l5row,&latitude,&longitude,&sun_elevation,&sun_azimuth,&c_year,&c_month,&c_day,&day,&month,&year,&decimal_hour);
- G_message("l5path=%i",l5path);
- G_message("l5row=%i",l5row);
- G_message("latitude=%f",latitude);
- G_message("longitude=%f",longitude);
- G_message("sun_elevation=%f",sun_elevation);
- G_message("sun_azimuth=%f",sun_azimuth);
- G_message("c_year=%i",c_year);
- G_message("c_month=%i",c_month);
- G_message("c_day=%i",c_day);
- G_message("year=%i",year);
- G_message("month=%i",month);
- G_message("day=%i",day);
- G_message("decimal hour=%f",decimal_hour);
+ /*Fetch parameters for DN2Rad2Ref correction */
+ l5_in_read(metfName, &l5path, &l5row, &latitude, &longitude,
+ &sun_elevation, &sun_azimuth, &c_year, &c_month, &c_day, &day,
+ &month, &year, &decimal_hour);
+ G_message("l5path=%i", l5path);
+ G_message("l5row=%i", l5row);
+ G_message("latitude=%f", latitude);
+ G_message("longitude=%f", longitude);
+ G_message("sun_elevation=%f", sun_elevation);
+ G_message("sun_azimuth=%f", sun_azimuth);
+ G_message("c_year=%i", c_year);
+ G_message("c_month=%i", c_month);
+ G_message("c_day=%i", c_day);
+ G_message("year=%i", year);
+ G_message("month=%i", month);
+ G_message("day=%i", day);
+ G_message("decimal hour=%f", decimal_hour);
+
/********************/
- /*Prepare the input file names */
+ /*Prepare the input file names */
+
/********************/
- doy = date2doy(day,month,year);
- G_message("doy=%i",doy);
- if(input1->answers){
- int index=1;
- for (; *ptr != NULL; ptr++){
- name = *ptr;
- if(index==1)
- sprintf(b1,"%s",name);
- if(index==2)
- sprintf(b2,"%s",name);
- if(index==3)
- sprintf(b3,"%s",name);
- if(index==4)
- sprintf(b4,"%s",name);
- if(index==5)
- sprintf(b5,"%s",name);
- if(index==6)
- sprintf(b6,"%s",name);
- if(index==7)
- sprintf(b7,"%s",name);
- index++;
- }
- } else {
- if(year<2000){
- temp = year - 1900;
- } else {
- temp = year - 2000;
- }
- if (temp >=10){
- sprintf(b1, "%s%d%s%d%s%d%d%s","LT5",l5path,"0",l5row,"00",temp,doy,"50_B1.tif");
- sprintf(b2, "%s%d%s%d%s%d%d%s","LT5",l5path,"0",l5row,"00",temp,doy,"50_B2.tif");
- sprintf(b3, "%s%d%s%d%s%d%d%s","LT5",l5path,"0",l5row,"00",temp,doy,"50_B3.tif");
- sprintf(b4, "%s%d%s%d%s%d%d%s","LT5",l5path,"0",l5row,"00",temp,doy,"50_B4.tif");
- sprintf(b5, "%s%d%s%d%s%d%d%s","LT5",l5path,"0",l5row,"00",temp,doy,"50_B5.tif");
- sprintf(b6, "%s%d%s%d%s%d%d%s","LT5",l5path,"0",l5row,"00",temp,doy,"50_B6.tif");
- sprintf(b7, "%s%d%s%d%s%d%d%s","LT5",l5path,"0",l5row,"00",temp,doy,"50_B7.tif");
- } else {
- sprintf(b1, "%s%d%s%d%s%d%d%s","LT5",l5path,"0",l5row,"000",temp,doy,"50_B1.tif");
- sprintf(b2, "%s%d%s%d%s%d%d%s","LT5",l5path,"0",l5row,"000",temp,doy,"50_B2.tif");
- sprintf(b3, "%s%d%s%d%s%d%d%s","LT5",l5path,"0",l5row,"000",temp,doy,"50_B3.tif");
- sprintf(b4, "%s%d%s%d%s%d%d%s","LT5",l5path,"0",l5row,"000",temp,doy,"50_B4.tif");
- sprintf(b5, "%s%d%s%d%s%d%d%s","LT5",l5path,"0",l5row,"000",temp,doy,"50_B5.tif");
- sprintf(b6, "%s%d%s%d%s%d%d%s","LT5",l5path,"0",l5row,"000",temp,doy,"50_B6.tif");
- sprintf(b7, "%s%d%s%d%s%d%d%s","LT5",l5path,"0",l5row,"000",temp,doy,"50_B7.tif");
- }
+ doy = date2doy(day, month, year);
+ G_message("doy=%i", doy);
+ if (input1->answers) {
+ int index = 1;
+
+ for (; *ptr != NULL; ptr++) {
+ name = *ptr;
+ if (index == 1)
+ sprintf(b1, "%s", name);
+ if (index == 2)
+ sprintf(b2, "%s", name);
+ if (index == 3)
+ sprintf(b3, "%s", name);
+ if (index == 4)
+ sprintf(b4, "%s", name);
+ if (index == 5)
+ sprintf(b5, "%s", name);
+ if (index == 6)
+ sprintf(b6, "%s", name);
+ if (index == 7)
+ sprintf(b7, "%s", name);
+ index++;
}
- G_message("b1_in:%s",b1);
- G_message("b2_in:%s",b2);
- G_message("b3_in:%s",b3);
- G_message("b4_in:%s",b4);
- G_message("b5_in:%s",b5);
- G_message("b6_in:%s",b6);
- G_message("b7_in:%s",b7);
-
+ }
+ else {
+ if (year < 2000) {
+ temp = year - 1900;
+ }
+ else {
+ temp = year - 2000;
+ }
+ if (temp >= 10) {
+ sprintf(b1, "%s%d%s%d%s%d%d%s", "LT5", l5path, "0", l5row, "00",
+ temp, doy, "50_B1.tif");
+ sprintf(b2, "%s%d%s%d%s%d%d%s", "LT5", l5path, "0", l5row, "00",
+ temp, doy, "50_B2.tif");
+ sprintf(b3, "%s%d%s%d%s%d%d%s", "LT5", l5path, "0", l5row, "00",
+ temp, doy, "50_B3.tif");
+ sprintf(b4, "%s%d%s%d%s%d%d%s", "LT5", l5path, "0", l5row, "00",
+ temp, doy, "50_B4.tif");
+ sprintf(b5, "%s%d%s%d%s%d%d%s", "LT5", l5path, "0", l5row, "00",
+ temp, doy, "50_B5.tif");
+ sprintf(b6, "%s%d%s%d%s%d%d%s", "LT5", l5path, "0", l5row, "00",
+ temp, doy, "50_B6.tif");
+ sprintf(b7, "%s%d%s%d%s%d%d%s", "LT5", l5path, "0", l5row, "00",
+ temp, doy, "50_B7.tif");
+ }
+ else {
+ sprintf(b1, "%s%d%s%d%s%d%d%s", "LT5", l5path, "0", l5row, "000",
+ temp, doy, "50_B1.tif");
+ sprintf(b2, "%s%d%s%d%s%d%d%s", "LT5", l5path, "0", l5row, "000",
+ temp, doy, "50_B2.tif");
+ sprintf(b3, "%s%d%s%d%s%d%d%s", "LT5", l5path, "0", l5row, "000",
+ temp, doy, "50_B3.tif");
+ sprintf(b4, "%s%d%s%d%s%d%d%s", "LT5", l5path, "0", l5row, "000",
+ temp, doy, "50_B4.tif");
+ sprintf(b5, "%s%d%s%d%s%d%d%s", "LT5", l5path, "0", l5row, "000",
+ temp, doy, "50_B5.tif");
+ sprintf(b6, "%s%d%s%d%s%d%d%s", "LT5", l5path, "0", l5row, "000",
+ temp, doy, "50_B6.tif");
+ sprintf(b7, "%s%d%s%d%s%d%d%s", "LT5", l5path, "0", l5row, "000",
+ temp, doy, "50_B7.tif");
+ }
+ }
+ G_message("b1_in:%s", b1);
+ G_message("b2_in:%s", b2);
+ G_message("b3_in:%s", b3);
+ G_message("b4_in:%s", b4);
+ G_message("b5_in:%s", b5);
+ G_message("b6_in:%s", b6);
+ G_message("b7_in:%s", b7);
+
/********************/
- /*Prepare the ouput file names */
+ /*Prepare the ouput file names */
+
/********************/
- G_message("result=%s",result);
- snprintf(result1, 80, "%s%s",result,".1");
- G_message("%s",result1);
- snprintf(result2, 80, "%s%s",result,".2");
- G_message("%s",result2);
- snprintf(result3, 80, "%s%s",result,".3");
- G_message("%s",result3);
- snprintf(result4, 80, "%s%s",result,".4");
- G_message("%s",result4);
- snprintf(result5, 80, "%s%s",result,".5");
- G_message("%s",result5);
- snprintf(result6, 80, "%s%s",result,".6");
- G_message("%s",result6);
- snprintf(result7, 80, "%s%s",result,".7");
- G_message("%s",result7);
+ G_message("result=%s", result);
+ snprintf(result1, 80, "%s%s", result, ".1");
+ G_message("%s", result1);
+ snprintf(result2, 80, "%s%s", result, ".2");
+ G_message("%s", result2);
+ snprintf(result3, 80, "%s%s", result, ".3");
+ G_message("%s", result3);
+ snprintf(result4, 80, "%s%s", result, ".4");
+ G_message("%s", result4);
+ snprintf(result5, 80, "%s%s", result, ".5");
+ G_message("%s", result5);
+ snprintf(result6, 80, "%s%s", result, ".6");
+ G_message("%s", result6);
+ snprintf(result7, 80, "%s%s", result, ".7");
+ G_message("%s", result7);
/********************/
- /*Prepare sun exo-atm irradiance*/
+ /*Prepare sun exo-atm irradiance */
+
/********************/
-
- kexo[0]=KEXO1;
- kexo[1]=KEXO2;
- kexo[2]=KEXO3;
- kexo[3]=KEXO4;
- kexo[4]=KEXO5;
- kexo[5]=KEXO6;/*filling only*/
- kexo[6]=KEXO7;
-
+
+ kexo[0] = KEXO1;
+ kexo[1] = KEXO2;
+ kexo[2] = KEXO3;
+ kexo[3] = KEXO4;
+ kexo[4] = KEXO5;
+ kexo[5] = KEXO6; /*filling only */
+ kexo[6] = KEXO7;
+
/***************************************************/
- /*Band1*/
- /* find map in mapset */
- mapset = G_find_cell2 (b1, "");
- if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b1);
- }
- if (G_legal_filename (b1) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b1);
- }
- infd[0] = G_open_cell_old (b1, mapset);
- /* Allocate input buffer */
- in_data_type[0] = G_raster_map_type(b1, mapset);
- if( (infd[0] = G_open_cell_old(b1,mapset)) < 0){
- G_fatal_error(_("Cannot open cell file [%s]"), b1);
- }
- if( (G_get_cellhd(b1,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b1);
- }
- inrast[0] = G_allocate_raster_buf(in_data_type[0]);
+ /*Band1 */
+ /* find map in mapset */
+ mapset = G_find_cell2(b1, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), b1);
+ }
+ if (G_legal_filename(b1) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), b1);
+ }
+ infd[0] = G_open_cell_old(b1, mapset);
+ /* Allocate input buffer */
+ in_data_type[0] = G_raster_map_type(b1, mapset);
+ if ((infd[0] = G_open_cell_old(b1, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), b1);
+ }
+ if ((G_get_cellhd(b1, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), b1);
+ }
+ inrast[0] = G_allocate_raster_buf(in_data_type[0]);
+
/***************************************************/
+
/***************************************************/
- /*Band2*/
- /* find map in mapset */
- mapset = G_find_cell2 (b2, "");
- if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b2);
- }
- if (G_legal_filename (b2) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b2);
- }
- infd[1] = G_open_cell_old (b2, mapset);
- /* Allocate input buffer */
- in_data_type[1] = G_raster_map_type(b2, mapset);
- if( (infd[1] = G_open_cell_old(b2,mapset)) < 0){
- G_fatal_error(_("Cannot open cell file [%s]"), b2);
- }
- if( (G_get_cellhd(b2,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b2);
- }
- inrast[1] = G_allocate_raster_buf(in_data_type[1]);
+ /*Band2 */
+ /* find map in mapset */
+ mapset = G_find_cell2(b2, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), b2);
+ }
+ if (G_legal_filename(b2) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), b2);
+ }
+ infd[1] = G_open_cell_old(b2, mapset);
+ /* Allocate input buffer */
+ in_data_type[1] = G_raster_map_type(b2, mapset);
+ if ((infd[1] = G_open_cell_old(b2, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), b2);
+ }
+ if ((G_get_cellhd(b2, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), b2);
+ }
+ inrast[1] = G_allocate_raster_buf(in_data_type[1]);
+
/***************************************************/
+
/***************************************************/
- /*Band3*/
- /* find map in mapset */
- mapset = G_find_cell2 (b3, "");
- if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b3);
- }
- if (G_legal_filename (b3) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b3);
- }
- infd[2] = G_open_cell_old (b3, mapset);
- /* Allocate input buffer */
- in_data_type[2] = G_raster_map_type(b3, mapset);
- if( (infd[2] = G_open_cell_old(b3,mapset)) < 0){
- G_fatal_error(_("Cannot open cell file [%s]"), b3);
- }
- if( (G_get_cellhd(b3,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b3);
- }
- inrast[2] = G_allocate_raster_buf(in_data_type[2]);
+ /*Band3 */
+ /* find map in mapset */
+ mapset = G_find_cell2(b3, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), b3);
+ }
+ if (G_legal_filename(b3) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), b3);
+ }
+ infd[2] = G_open_cell_old(b3, mapset);
+ /* Allocate input buffer */
+ in_data_type[2] = G_raster_map_type(b3, mapset);
+ if ((infd[2] = G_open_cell_old(b3, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), b3);
+ }
+ if ((G_get_cellhd(b3, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), b3);
+ }
+ inrast[2] = G_allocate_raster_buf(in_data_type[2]);
+
/***************************************************/
+
/***************************************************/
- /*Band4*/
- /* find map in mapset */
- mapset = G_find_cell2 (b4, "");
- if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b4);
- }
- if (G_legal_filename (b4) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b4);
- }
- infd[3] = G_open_cell_old (b4, mapset);
- /* Allocate input buffer */
- in_data_type[3] = G_raster_map_type(b4, mapset);
- if( (infd[3] = G_open_cell_old(b4,mapset)) < 0){
- G_fatal_error(_("Cannot open cell file [%s]"), b4);
- }
- if( (G_get_cellhd(b4,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b4);
- }
- inrast[3] = G_allocate_raster_buf(in_data_type[3]);
+ /*Band4 */
+ /* find map in mapset */
+ mapset = G_find_cell2(b4, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), b4);
+ }
+ if (G_legal_filename(b4) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), b4);
+ }
+ infd[3] = G_open_cell_old(b4, mapset);
+ /* Allocate input buffer */
+ in_data_type[3] = G_raster_map_type(b4, mapset);
+ if ((infd[3] = G_open_cell_old(b4, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), b4);
+ }
+ if ((G_get_cellhd(b4, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), b4);
+ }
+ inrast[3] = G_allocate_raster_buf(in_data_type[3]);
+
/***************************************************/
+
/***************************************************/
- /*Band5*/
- /* find map in mapset */
- mapset = G_find_cell2 (b5, "");
- if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b5);
- }
- if (G_legal_filename (b5) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b5);
- }
- infd[4] = G_open_cell_old (b5, mapset);
- /* Allocate input buffer */
- in_data_type[4] = G_raster_map_type(b5, mapset);
- if( (infd[4] = G_open_cell_old(b5,mapset)) < 0){
- G_fatal_error(_("Cannot open cell file [%s]"), b5);
- }
- if( (G_get_cellhd(b5,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b5);
- }
- inrast[4] = G_allocate_raster_buf(in_data_type[4]);
+ /*Band5 */
+ /* find map in mapset */
+ mapset = G_find_cell2(b5, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), b5);
+ }
+ if (G_legal_filename(b5) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), b5);
+ }
+ infd[4] = G_open_cell_old(b5, mapset);
+ /* Allocate input buffer */
+ in_data_type[4] = G_raster_map_type(b5, mapset);
+ if ((infd[4] = G_open_cell_old(b5, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), b5);
+ }
+ if ((G_get_cellhd(b5, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), b5);
+ }
+ inrast[4] = G_allocate_raster_buf(in_data_type[4]);
+
/***************************************************/
+
/***************************************************/
- /*Band6*/
- /* find map in mapset */
- mapset = G_find_cell2 (b6, "");
- if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b6);
- }
- if (G_legal_filename (b6) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b6);
- }
- infd[5] = G_open_cell_old (b6, mapset);
- /* Allocate input buffer */
- in_data_type[5] = G_raster_map_type(b6, mapset);
- if( (infd[5] = G_open_cell_old(b6,mapset)) < 0){
- G_fatal_error(_("Cannot open cell file [%s]"), b6);
- }
- if( (G_get_cellhd(b6,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b6);
- }
- inrast[5] = G_allocate_raster_buf(in_data_type[5]);
+ /*Band6 */
+ /* find map in mapset */
+ mapset = G_find_cell2(b6, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), b6);
+ }
+ if (G_legal_filename(b6) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), b6);
+ }
+ infd[5] = G_open_cell_old(b6, mapset);
+ /* Allocate input buffer */
+ in_data_type[5] = G_raster_map_type(b6, mapset);
+ if ((infd[5] = G_open_cell_old(b6, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), b6);
+ }
+ if ((G_get_cellhd(b6, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), b6);
+ }
+ inrast[5] = G_allocate_raster_buf(in_data_type[5]);
+
/***************************************************/
+
/***************************************************/
- /*Band7*/
- /* find map in mapset */
- mapset = G_find_cell2 (b7, "");
- if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b7);
- }
- if (G_legal_filename (b7) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b7);
- }
- infd[6] = G_open_cell_old (b7, mapset);
- /* Allocate input buffer */
- in_data_type[6] = G_raster_map_type(b7, mapset);
- if( (infd[6] = G_open_cell_old(b7,mapset)) < 0){
- G_fatal_error(_("Cannot open cell file [%s]"), b7);
- }
- if( (G_get_cellhd(b7,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b7);
- }
- inrast[6] = G_allocate_raster_buf(in_data_type[6]);
+ /*Band7 */
+ /* find map in mapset */
+ mapset = G_find_cell2(b7, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), b7);
+ }
+ if (G_legal_filename(b7) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), b7);
+ }
+ infd[6] = G_open_cell_old(b7, mapset);
+ /* Allocate input buffer */
+ in_data_type[6] = G_raster_map_type(b7, mapset);
+ if ((infd[6] = G_open_cell_old(b7, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), b7);
+ }
+ if ((G_get_cellhd(b7, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), b7);
+ }
+ inrast[6] = G_allocate_raster_buf(in_data_type[6]);
+
/***************************************************/
+
/***************************************************/
- /* Allocate output buffer, use input map data_type */
- nrows = G_window_rows();
- ncols = G_window_cols();
- out_data_type=DCELL_TYPE;
- for(i=0;i<MAXFILES;i++){
- outrast[i] = G_allocate_raster_buf(out_data_type);
+ /* Allocate output buffer, use input map data_type */
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ out_data_type = DCELL_TYPE;
+ for (i = 0; i < MAXFILES; i++) {
+ outrast[i] = G_allocate_raster_buf(out_data_type);
+ }
+ if ((outfd[0] = G_open_raster_new(result1, 1)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result1);
+ if ((outfd[1] = G_open_raster_new(result2, 1)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result2);
+ if ((outfd[2] = G_open_raster_new(result3, 1)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result3);
+ if ((outfd[3] = G_open_raster_new(result4, 1)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result4);
+ if ((outfd[4] = G_open_raster_new(result5, 1)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result5);
+ if ((outfd[5] = G_open_raster_new(result6, 1)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result6);
+ if ((outfd[6] = G_open_raster_new(result7, 1)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result7);
+ /* Process pixels */
+ DCELL dout[MAXFILES];
+
+ DCELL d[MAXFILES];
+
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 2);
+ /* read input map */
+ for (i = 0; i < MAXFILES; i++) {
+ if ((G_get_raster_row(infd[i], inrast[i], row, in_data_type[i])) <
+ 0) {
+ G_fatal_error(_("Could not read from <%s>"), i);
+ }
}
- if ( (outfd[0] = G_open_raster_new (result1,1)) < 0)
- G_fatal_error (_("Could not open <%s>"),result1);
- if ( (outfd[1] = G_open_raster_new (result2,1)) < 0)
- G_fatal_error (_("Could not open <%s>"),result2);
- if ( (outfd[2] = G_open_raster_new (result3,1)) < 0)
- G_fatal_error (_("Could not open <%s>"),result3);
- if ( (outfd[3] = G_open_raster_new (result4,1)) < 0)
- G_fatal_error (_("Could not open <%s>"),result4);
- if ( (outfd[4] = G_open_raster_new (result5,1)) < 0)
- G_fatal_error (_("Could not open <%s>"),result5);
- if ( (outfd[5] = G_open_raster_new (result6,1)) < 0)
- G_fatal_error (_("Could not open <%s>"),result6);
- if ( (outfd[6] = G_open_raster_new (result7,1)) < 0)
- G_fatal_error (_("Could not open <%s>"),result7);
- /* Process pixels */
- DCELL dout[MAXFILES];
- DCELL d[MAXFILES];
- for (row = 0; row < nrows; row++)
- {
- G_percent (row, nrows, 2);
- /* read input map */
- for (i=0;i<MAXFILES;i++)
- {
- if((G_get_raster_row(infd[i],inrast[i],row,in_data_type[i])) < 0){
- G_fatal_error (_("Could not read from <%s>"),i);
- }
+ /*process the data */
+ for (col = 0; col < ncols; col++) {
+ for (i = 0; i < MAXFILES; i++) {
+ switch (in_data_type[i]) {
+ case CELL_TYPE:
+ d[i] = (double)((CELL *) inrast[i])[col];
+ break;
+ case FCELL_TYPE:
+ d[i] = (double)((FCELL *) inrast[i])[col];
+ break;
+ case DCELL_TYPE:
+ d[i] = (double)((DCELL *) inrast[i])[col];
+ break;
}
- /*process the data */
- for (col=0; col < ncols; col++)
- {
- for(i=0;i<MAXFILES;i++)
- {
- switch(in_data_type[i]){
- case CELL_TYPE:
- d[i] = (double) ((CELL *) inrast[i])[col];
- break;
- case FCELL_TYPE:
- d[i] = (double) ((FCELL *) inrast[i])[col];
- break;
- case DCELL_TYPE:
- d[i] = (double) ((DCELL *) inrast[i])[col];
- break;
- }
- dout[i]=dn2rad_landsat5(c_year,c_month,c_day,year,month,day,i+1,d[i]);
- if(i==5){/*if band 6, process brightness temperature*/
- if(dout[i]<=0.0){
- G_set_d_null_value(&dout[i],1);
- } else {
- dout[i]=tempk_landsat5(dout[i]);
- }
- } else {/*process reflectance*/
- dout[i]=rad2ref_landsat5(dout[i],doy,sun_elevation,kexo[i]);
- }
- outrast[i][col] = dout[i];
- }
+ dout[i] =
+ dn2rad_landsat5(c_year, c_month, c_day, year, month, day,
+ i + 1, d[i]);
+ if (i == 5) { /*if band 6, process brightness temperature */
+ if (dout[i] <= 0.0) {
+ G_set_d_null_value(&dout[i], 1);
+ }
+ else {
+ dout[i] = tempk_landsat5(dout[i]);
+ }
}
- for(i=0;i<MAXFILES;i++){
- if (G_put_raster_row (outfd[i], outrast[i], out_data_type) < 0)
- G_fatal_error (_("Cannot write to <%s.%i>"),result,i+1);
+ else { /*process reflectance */
+ dout[i] =
+ rad2ref_landsat5(dout[i], doy, sun_elevation,
+ kexo[i]);
}
+ outrast[i][col] = dout[i];
+ }
}
- for (i=0;i<MAXFILES;i++)
- {
- G_free (inrast[i]);
- G_close_cell (infd[i]);
- G_free (outrast[i]);
- G_close_cell (outfd[i]);
+ for (i = 0; i < MAXFILES; i++) {
+ if (G_put_raster_row(outfd[i], outrast[i], out_data_type) < 0)
+ G_fatal_error(_("Cannot write to <%s.%i>"), result, i + 1);
}
+ }
+ for (i = 0; i < MAXFILES; i++) {
+ G_free(inrast[i]);
+ G_close_cell(infd[i]);
+ G_free(outrast[i]);
+ G_close_cell(outfd[i]);
+ }
- return 0;
+ return 0;
}
Modified: grass-addons/gipe/i.dn2full.l5/rad2ref_landsat5.c
===================================================================
--- grass-addons/gipe/i.dn2full.l5/rad2ref_landsat5.c 2008-08-21 23:06:30 UTC (rev 32985)
+++ grass-addons/gipe/i.dn2full.l5/rad2ref_landsat5.c 2008-08-21 23:12:57 UTC (rev 32986)
@@ -1,18 +1,20 @@
-#include<stdio.h>
-#include<math.h>
-#include<stdlib.h>
-
-#define PI 3.1415926
-
-/* Conversion of Radiance to Reflectance for Landsat 5TM */
-
-double rad2ref_landsat5( double radiance, double doy,double sun_elevation, double k_exo )
-{
- double result, ds;
-
- ds = (1+0.01672*sin(2*PI*(doy-93.5)/365));
- result = (radiance/((cos((90-sun_elevation)*PI/180)/(PI*ds*ds))*k_exo));
-
- return result;
-}
-
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+
+#define PI 3.1415926
+
+ /* Conversion of Radiance to Reflectance for Landsat 5TM */
+double rad2ref_landsat5(double radiance, double doy, double sun_elevation,
+ double k_exo)
+{
+ double result, ds;
+
+ ds = (1 + 0.01672 * sin(2 * PI * (doy - 93.5) / 365));
+ result =
+ (radiance /
+ ((cos((90 - sun_elevation) * PI / 180) / (PI * ds * ds)) * k_exo));
+ return result;
+}
+
+
Modified: grass-addons/gipe/i.dn2full.l5/tempk_landsat5.c
===================================================================
--- grass-addons/gipe/i.dn2full.l5/tempk_landsat5.c 2008-08-21 23:06:30 UTC (rev 32985)
+++ grass-addons/gipe/i.dn2full.l5/tempk_landsat5.c 2008-08-21 23:12:57 UTC (rev 32986)
@@ -1,18 +1,17 @@
-#include<stdio.h>
-#include<math.h>
-#include<stdlib.h>
-
-/*
- * Surface temperature for Landsat 5TM
- * Schneider and Mauser, 1996
- */
-
-double tempk_landsat5( double l6 )
-{
- double result;
-
- result = 1260.56 / (log ((607.76 / (l6)) + 1.0)) ;
-
- return result;
-}
-
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+
+ /*
+ * Surface temperature for Landsat 5TM
+ * Schneider and Mauser, 1996
+ */
+double tempk_landsat5(double l6)
+{
+ double result;
+
+ result = 1260.56 / (log((607.76 / (l6)) + 1.0));
+ return result;
+}
+
+
Modified: grass-addons/gipe/i.dn2full.l7/date2doy.c
===================================================================
--- grass-addons/gipe/i.dn2full.l7/date2doy.c 2008-08-21 23:06:30 UTC (rev 32985)
+++ grass-addons/gipe/i.dn2full.l7/date2doy.c 2008-08-21 23:12:57 UTC (rev 32986)
@@ -1,66 +1,70 @@
+
/*********************************************/
-/*This program converts day/month/year to doy*/
+/*This program converts day/month/year to doy */
+
/*********************************************/
+
/*********************************************/
int date2doy(int day, int month, int year)
{
- int leap=0;
- int day_month_tot=0;
- int doy;
+ int leap = 0;
- doy=0;
+ int day_month_tot = 0;
-/*printf("Date is %i/%i/%i\n", day, month, year);*/
+ int doy;
- if (month == 1) {
- day_month_tot = 0;
- }
- else if (month == 2) {
- day_month_tot = 31;
- }
- else if (month == 3) {
- day_month_tot = 59;
- }
- else if (month == 4) {
- day_month_tot = 90;
- }
- else if (month == 5) {
- day_month_tot = 120;
- }
- else if (month == 6) {
- day_month_tot = 151;
- }
- else if (month == 7) {
- day_month_tot = 181;
- }
- else if (month == 8) {
- day_month_tot = 212;
- }
- else if (month == 9) {
- day_month_tot = 243;
- }
- else if (month == 10) {
- day_month_tot = 273;
- }
- else if (month == 11) {
- day_month_tot = 304;
- }
- else if (month == 12) {
- day_month_tot = 334;
- }
-
- /* Leap year if dividing by 4 leads % 0.0*/
-
- if (year/4*4 == year) {
- leap = 1;
- }
-
- doy=day_month_tot+day;
- if(doy>59){
- doy=day_month_tot+day+leap;
- }
+ doy = 0;
- return(doy);
+ /*printf("Date is %i/%i/%i\n", day, month, year); */
+
+ if (month == 1) {
+ day_month_tot = 0;
+ }
+ else if (month == 2) {
+ day_month_tot = 31;
+ }
+ else if (month == 3) {
+ day_month_tot = 59;
+ }
+ else if (month == 4) {
+ day_month_tot = 90;
+ }
+ else if (month == 5) {
+ day_month_tot = 120;
+ }
+ else if (month == 6) {
+ day_month_tot = 151;
+ }
+ else if (month == 7) {
+ day_month_tot = 181;
+ }
+ else if (month == 8) {
+ day_month_tot = 212;
+ }
+ else if (month == 9) {
+ day_month_tot = 243;
+ }
+ else if (month == 10) {
+ day_month_tot = 273;
+ }
+ else if (month == 11) {
+ day_month_tot = 304;
+ }
+ else if (month == 12) {
+ day_month_tot = 334;
+ }
+
+ /* Leap year if dividing by 4 leads % 0.0 */
+
+ if (year / 4 * 4 == year) {
+ leap = 1;
+ }
+
+ doy = day_month_tot + day;
+ if (doy > 59) {
+ doy = day_month_tot + day + leap;
+ }
+
+ return (doy);
}
-
Modified: grass-addons/gipe/i.dn2full.l7/main.c
===================================================================
--- grass-addons/gipe/i.dn2full.l7/main.c 2008-08-21 23:06:30 UTC (rev 32985)
+++ grass-addons/gipe/i.dn2full.l7/main.c 2008-08-21 23:12:57 UTC (rev 32986)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: i.dn2full.l7
@@ -22,391 +23,451 @@
#define MAXFILES 9
-/*sun exo-atmospheric irradiance*/
+/*sun exo-atmospheric irradiance */
#define KEXO1 1969.0
#define KEXO2 1840.0
#define KEXO3 1551.0
#define KEXO4 1044.0
#define KEXO5 225.7
#define KEXO7 82.07
-#define KEXO8 1385.64 /*to find the real value in the internet*/
+#define KEXO8 1385.64 /*to find the real value in the internet */
#define PI 3.1415926
-int l7_in_read(char *metfName, char *b1, char *b2, char *b3, char *b4, char *b5, char *b61, char *b62, char *b7, char *b8, double *lmin,double *lmax,double *qcalmin,double *qcalmax,double *sun_elevation, double *sun_azimuth,int *day, int *month, int *year);
+int l7_in_read(char *metfName, char *b1, char *b2, char *b3, char *b4,
+ char *b5, char *b61, char *b62, char *b7, char *b8,
+ double *lmin, double *lmax, double *qcalmin, double *qcalmax,
+ double *sun_elevation, double *sun_azimuth, int *day,
+ int *month, int *year);
int date2doy(int day, int month, int year);
-double dn2rad_landsat7( double Lmin, double LMax, double QCalMax, double QCalmin, int DN );
-double rad2ref_landsat7( double radiance, double doy,double sun_elevation, double k_exo );
-double tempk_landsat7( double l6 );
-int
-main(int argc, char *argv[])
+double dn2rad_landsat7(double Lmin, double LMax, double QCalMax,
+ double QCalmin, int DN);
+double rad2ref_landsat7(double radiance, double doy, double sun_elevation,
+ double k_exo);
+double tempk_landsat7(double l6);
+
+int main(int argc, char *argv[])
{
- struct Cell_head cellhd;/*region+header info*/
- char *mapset; /*mapset name*/
- int nrows, ncols;
- int row,col;
+ struct Cell_head cellhd; /*region+header info */
- struct GModule *module;
- struct Option *input,*output;
-
- struct History history; /*metadata*/
+ char *mapset; /*mapset name */
+
+ int nrows, ncols;
+
+ int row, col;
+
+ struct GModule *module;
+
+ struct Option *input, *output;
+
+ struct History history; /*metadata */
+
/************************************/
- /* FMEO Declarations*****************/
- char history_buf[200];
- char *name; /*input raster name*/
- char *result; /*output raster name*/
- /*Prepare new names for output files*/
- char result1[80], result2[80], result3[80], result4[80];
- char result5[80], result61[80], result62[80], result7[80],result8[80] ;
-
- /*File Descriptors*/
- int infd[MAXFILES];
- int outfd[MAXFILES];
+ /* FMEO Declarations**************** */
+ char history_buf[200];
- char **names;
- char **ptr;
-
- int i=0,j=0;
-
- void *inrast[MAXFILES];
- DCELL *outrast[MAXFILES];
- RASTER_MAP_TYPE in_data_type[MAXFILES];
- RASTER_MAP_TYPE out_data_type = DCELL_TYPE; /* 0=numbers 1=text */
+ char *name; /*input raster name */
- double kexo[MAXFILES];
- /*Metfile*/
- char *metfName; /*met file, header in text format*/
- char b1[80],b2[80],b3[80];
- char b4[80],b5[80],b61[80];
- char b62[80],b7[80],b8[80];/*Load .tif names*/
- double lmin[MAXFILES];
- double lmax[MAXFILES];
- double qcalmin[MAXFILES];
- double qcalmax[MAXFILES];
- double sun_elevation;
- double sun_azimuth;/*not useful here, only for parser()*/
- int day,month,year;
- /*EndofMetfile*/
- int doy;
+ char *result; /*output raster name */
+
+ /*Prepare new names for output files */
+ char result1[80], result2[80], result3[80], result4[80];
+
+ char result5[80], result61[80], result62[80], result7[80], result8[80];
+
+ /*File Descriptors */
+ int infd[MAXFILES];
+
+ int outfd[MAXFILES];
+
+ char **names;
+
+ char **ptr;
+
+ int i = 0, j = 0;
+
+ void *inrast[MAXFILES];
+
+ DCELL *outrast[MAXFILES];
+
+ RASTER_MAP_TYPE in_data_type[MAXFILES];
+
+ RASTER_MAP_TYPE out_data_type = DCELL_TYPE; /* 0=numbers 1=text */
+
+ double kexo[MAXFILES];
+
+ /*Metfile */
+ char *metfName; /*met file, header in text format */
+
+ char b1[80], b2[80], b3[80];
+
+ char b4[80], b5[80], b61[80];
+
+ char b62[80], b7[80], b8[80]; /*Load .tif names */
+
+ double lmin[MAXFILES];
+
+ double lmax[MAXFILES];
+
+ double qcalmin[MAXFILES];
+
+ double qcalmax[MAXFILES];
+
+ double sun_elevation;
+
+ double sun_azimuth; /*not useful here, only for parser() */
+
+ int day, month, year;
+
+ /*EndofMetfile */
+ int doy;
+
/************************************/
- G_gisinit(argv[0]);
+ G_gisinit(argv[0]);
- module = G_define_module();
- module->keywords = _("DN, reflectance, temperature, import");
- module->description =
- _("Calculates Top of Atmosphere Reflectance/Temperature from Landsat 7 DN.\n");
+ module = G_define_module();
+ module->keywords = _("DN, reflectance, temperature, import");
+ module->description =
+ _("Calculates Top of Atmosphere Reflectance/Temperature from Landsat 7 DN.\n");
- /* Define the different options */
- input = G_define_option() ;
- input->key = _("metfile");
- input->type = TYPE_STRING;
- input->required = YES;
- input->gisprompt = _("old_file,file,file");
- input->description= _("Landsat 7ETM+ Header File (.met)");
+ /* Define the different options */
+ input = G_define_option();
+ input->key = _("metfile");
+ input->type = TYPE_STRING;
+ input->required = YES;
+ input->gisprompt = _("old_file,file,file");
+ input->description = _("Landsat 7ETM+ Header File (.met)");
- output = G_define_standard_option(G_OPT_R_OUTPUT) ;
- output->key = _("output");
- output->description= _("Base name of the output layers (will add .x)");
+ output = G_define_standard_option(G_OPT_R_OUTPUT);
+ output->key = _("output");
+ output->description = _("Base name of the output layers (will add .x)");
/********************/
- if (G_parser(argc, argv))
- exit (-1);
-
- metfName = input->answer;
- result = output->answer;
+ if (G_parser(argc, argv))
+ exit(-1);
+ metfName = input->answer;
+ result = output->answer;
+
/******************************************/
- /*Fetch parameters for DN2Rad2Ref correction*/
- l7_in_read(metfName,b1,b2,b3,b4,b5,b61,b62,b7,b8,lmin,lmax,qcalmin,qcalmax,&sun_elevation,&sun_azimuth,&day,&month,&year);
+ /*Fetch parameters for DN2Rad2Ref correction */
+ l7_in_read(metfName, b1, b2, b3, b4, b5, b61, b62, b7, b8, lmin, lmax,
+ qcalmin, qcalmax, &sun_elevation, &sun_azimuth, &day, &month,
+ &year);
+
/********************/
- /*Prepare the ouput file names */
+ /*Prepare the ouput file names */
+
/********************/
- snprintf(result1, 80, "%s%s",result,".1");
- snprintf(result2, 80, "%s%s",result,".2");
- snprintf(result3, 80, "%s%s",result,".3");
- snprintf(result4, 80, "%s%s",result,".4");
- snprintf(result5, 80, "%s%s",result,".5");
- snprintf(result61, 80, "%s%s",result,".61");
- snprintf(result62, 80, "%s%s",result,".62");
- snprintf(result7, 80, "%s%s",result,".7");
- snprintf(result8, 80, "%s%s",result,".8");
+ snprintf(result1, 80, "%s%s", result, ".1");
+ snprintf(result2, 80, "%s%s", result, ".2");
+ snprintf(result3, 80, "%s%s", result, ".3");
+ snprintf(result4, 80, "%s%s", result, ".4");
+ snprintf(result5, 80, "%s%s", result, ".5");
+ snprintf(result61, 80, "%s%s", result, ".61");
+ snprintf(result62, 80, "%s%s", result, ".62");
+ snprintf(result7, 80, "%s%s", result, ".7");
+ snprintf(result8, 80, "%s%s", result, ".8");
/********************/
- /*Prepare sun exo-atm irradiance*/
+ /*Prepare sun exo-atm irradiance */
+
/********************/
-
- kexo[0]=KEXO1;
- kexo[1]=KEXO2;
- kexo[2]=KEXO3;
- kexo[3]=KEXO4;
- kexo[4]=KEXO5;
- kexo[5]=0.0;/*filling*/
- kexo[6]=0.0;/*filling*/
- kexo[7]=KEXO7;
- kexo[8]=KEXO8;
-
+
+ kexo[0] = KEXO1;
+ kexo[1] = KEXO2;
+ kexo[2] = KEXO3;
+ kexo[3] = KEXO4;
+ kexo[4] = KEXO5;
+ kexo[5] = 0.0; /*filling */
+ kexo[6] = 0.0; /*filling */
+ kexo[7] = KEXO7;
+ kexo[8] = KEXO8;
+
/***************************************************/
- /*Band1*/
- /* find map in mapset */
- mapset = G_find_cell2 (b1, "");
- if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b1);
- }
- if (G_legal_filename (b1) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b1);
- }
- infd[0] = G_open_cell_old (b1, mapset);
- /* Allocate input buffer */
- in_data_type[0] = G_raster_map_type(b1, mapset);
- if( (infd[0] = G_open_cell_old(b1,mapset)) < 0){
- G_fatal_error(_("Cannot open cell file [%s]"), b1);
- }
- if( (G_get_cellhd(b1,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b1);
- }
- inrast[0] = G_allocate_raster_buf(in_data_type[0]);
+ /*Band1 */
+ /* find map in mapset */
+ mapset = G_find_cell2(b1, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), b1);
+ }
+ if (G_legal_filename(b1) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), b1);
+ }
+ infd[0] = G_open_cell_old(b1, mapset);
+ /* Allocate input buffer */
+ in_data_type[0] = G_raster_map_type(b1, mapset);
+ if ((infd[0] = G_open_cell_old(b1, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), b1);
+ }
+ if ((G_get_cellhd(b1, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), b1);
+ }
+ inrast[0] = G_allocate_raster_buf(in_data_type[0]);
+
/***************************************************/
+
/***************************************************/
- /*Band2*/
- /* find map in mapset */
- mapset = G_find_cell2 (b2, "");
- if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b2);
- }
- if (G_legal_filename (b2) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b2);
- }
- infd[1] = G_open_cell_old (b2, mapset);
- /* Allocate input buffer */
- in_data_type[1] = G_raster_map_type(b2, mapset);
- if( (infd[1] = G_open_cell_old(b2,mapset)) < 0){
- G_fatal_error(_("Cannot open cell file [%s]"), b2);
- }
- if( (G_get_cellhd(b2,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b2);
- }
- inrast[1] = G_allocate_raster_buf(in_data_type[1]);
+ /*Band2 */
+ /* find map in mapset */
+ mapset = G_find_cell2(b2, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), b2);
+ }
+ if (G_legal_filename(b2) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), b2);
+ }
+ infd[1] = G_open_cell_old(b2, mapset);
+ /* Allocate input buffer */
+ in_data_type[1] = G_raster_map_type(b2, mapset);
+ if ((infd[1] = G_open_cell_old(b2, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), b2);
+ }
+ if ((G_get_cellhd(b2, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), b2);
+ }
+ inrast[1] = G_allocate_raster_buf(in_data_type[1]);
+
/***************************************************/
+
/***************************************************/
- /*Band3*/
- /* find map in mapset */
- mapset = G_find_cell2 (b3, "");
- if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b3);
- }
- if (G_legal_filename (b3) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b3);
- }
- infd[2] = G_open_cell_old (b3, mapset);
- /* Allocate input buffer */
- in_data_type[2] = G_raster_map_type(b3, mapset);
- if( (infd[2] = G_open_cell_old(b3,mapset)) < 0){
- G_fatal_error(_("Cannot open cell file [%s]"), b3);
- }
- if( (G_get_cellhd(b3,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b3);
- }
- inrast[2] = G_allocate_raster_buf(in_data_type[2]);
+ /*Band3 */
+ /* find map in mapset */
+ mapset = G_find_cell2(b3, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), b3);
+ }
+ if (G_legal_filename(b3) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), b3);
+ }
+ infd[2] = G_open_cell_old(b3, mapset);
+ /* Allocate input buffer */
+ in_data_type[2] = G_raster_map_type(b3, mapset);
+ if ((infd[2] = G_open_cell_old(b3, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), b3);
+ }
+ if ((G_get_cellhd(b3, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), b3);
+ }
+ inrast[2] = G_allocate_raster_buf(in_data_type[2]);
+
/***************************************************/
+
/***************************************************/
- /*Band4*/
- /* find map in mapset */
- mapset = G_find_cell2 (b4, "");
- if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b4);
- }
- if (G_legal_filename (b4) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b4);
- }
- infd[3] = G_open_cell_old (b4, mapset);
- /* Allocate input buffer */
- in_data_type[3] = G_raster_map_type(b4, mapset);
- if( (infd[3] = G_open_cell_old(b4,mapset)) < 0){
- G_fatal_error(_("Cannot open cell file [%s]"), b4);
- }
- if( (G_get_cellhd(b4,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b4);
- }
- inrast[3] = G_allocate_raster_buf(in_data_type[3]);
+ /*Band4 */
+ /* find map in mapset */
+ mapset = G_find_cell2(b4, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), b4);
+ }
+ if (G_legal_filename(b4) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), b4);
+ }
+ infd[3] = G_open_cell_old(b4, mapset);
+ /* Allocate input buffer */
+ in_data_type[3] = G_raster_map_type(b4, mapset);
+ if ((infd[3] = G_open_cell_old(b4, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), b4);
+ }
+ if ((G_get_cellhd(b4, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), b4);
+ }
+ inrast[3] = G_allocate_raster_buf(in_data_type[3]);
+
/***************************************************/
+
/***************************************************/
- /*Band5*/
- /* find map in mapset */
- mapset = G_find_cell2 (b5, "");
- if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b5);
- }
- if (G_legal_filename (b5) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b5);
- }
- infd[4] = G_open_cell_old (b5, mapset);
- /* Allocate input buffer */
- in_data_type[4] = G_raster_map_type(b5, mapset);
- if( (infd[4] = G_open_cell_old(b5,mapset)) < 0){
- G_fatal_error(_("Cannot open cell file [%s]"), b5);
- }
- if( (G_get_cellhd(b5,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b5);
- }
- inrast[4] = G_allocate_raster_buf(in_data_type[4]);
+ /*Band5 */
+ /* find map in mapset */
+ mapset = G_find_cell2(b5, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), b5);
+ }
+ if (G_legal_filename(b5) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), b5);
+ }
+ infd[4] = G_open_cell_old(b5, mapset);
+ /* Allocate input buffer */
+ in_data_type[4] = G_raster_map_type(b5, mapset);
+ if ((infd[4] = G_open_cell_old(b5, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), b5);
+ }
+ if ((G_get_cellhd(b5, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), b5);
+ }
+ inrast[4] = G_allocate_raster_buf(in_data_type[4]);
+
/***************************************************/
+
/***************************************************/
- /*Band61*/
- /* find map in mapset */
- mapset = G_find_cell2 (b61, "");
- if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b61);
- }
- if (G_legal_filename (b61) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b61);
- }
- infd[5] = G_open_cell_old (b61, mapset);
- /* Allocate input buffer */
- in_data_type[5] = G_raster_map_type(b61, mapset);
- if( (infd[5] = G_open_cell_old(b61,mapset)) < 0){
- G_fatal_error(_("Cannot open cell file [%s]"), b61);
- }
- if( (G_get_cellhd(b61,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b61);
- }
- inrast[5] = G_allocate_raster_buf(in_data_type[5]);
+ /*Band61 */
+ /* find map in mapset */
+ mapset = G_find_cell2(b61, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), b61);
+ }
+ if (G_legal_filename(b61) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), b61);
+ }
+ infd[5] = G_open_cell_old(b61, mapset);
+ /* Allocate input buffer */
+ in_data_type[5] = G_raster_map_type(b61, mapset);
+ if ((infd[5] = G_open_cell_old(b61, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), b61);
+ }
+ if ((G_get_cellhd(b61, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), b61);
+ }
+ inrast[5] = G_allocate_raster_buf(in_data_type[5]);
+
/***************************************************/
+
/***************************************************/
- /*Band62*/
- /* find map in mapset */
- mapset = G_find_cell2 (b62, "");
- if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b62);
- }
- if (G_legal_filename (b62) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b62);
- }
- infd[6] = G_open_cell_old (b62, mapset);
- /* Allocate input buffer */
- in_data_type[6] = G_raster_map_type(b62, mapset);
- if( (infd[6] = G_open_cell_old(b62,mapset)) < 0){
- G_fatal_error(_("Cannot open cell file [%s]"), b62);
- }
- if( (G_get_cellhd(b62,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b62);
- }
- inrast[6] = G_allocate_raster_buf(in_data_type[6]);
+ /*Band62 */
+ /* find map in mapset */
+ mapset = G_find_cell2(b62, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), b62);
+ }
+ if (G_legal_filename(b62) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), b62);
+ }
+ infd[6] = G_open_cell_old(b62, mapset);
+ /* Allocate input buffer */
+ in_data_type[6] = G_raster_map_type(b62, mapset);
+ if ((infd[6] = G_open_cell_old(b62, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), b62);
+ }
+ if ((G_get_cellhd(b62, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), b62);
+ }
+ inrast[6] = G_allocate_raster_buf(in_data_type[6]);
+
/***************************************************/
+
/***************************************************/
- /*Band7*/
- /* find map in mapset */
- mapset = G_find_cell2 (b7, "");
- if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b7);
- }
- if (G_legal_filename (b7) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b7);
- }
- infd[7] = G_open_cell_old (b7, mapset);
- /* Allocate input buffer */
- in_data_type[7] = G_raster_map_type(b7, mapset);
- if( (infd[7] = G_open_cell_old(b7,mapset)) < 0){
- G_fatal_error(_("Cannot open cell file [%s]"), b7);
- }
- if( (G_get_cellhd(b7,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b7);
- }
- inrast[7] = G_allocate_raster_buf(in_data_type[7]);
+ /*Band7 */
+ /* find map in mapset */
+ mapset = G_find_cell2(b7, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), b7);
+ }
+ if (G_legal_filename(b7) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), b7);
+ }
+ infd[7] = G_open_cell_old(b7, mapset);
+ /* Allocate input buffer */
+ in_data_type[7] = G_raster_map_type(b7, mapset);
+ if ((infd[7] = G_open_cell_old(b7, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), b7);
+ }
+ if ((G_get_cellhd(b7, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), b7);
+ }
+ inrast[7] = G_allocate_raster_buf(in_data_type[7]);
+
/***************************************************/
+
/***************************************************/
- /*Band8*/
- /* find map in mapset */
- mapset = G_find_cell2 (b8, "");
- if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b8);
- }
- if (G_legal_filename (b8) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b8);
- }
- infd[8] = G_open_cell_old (b8, mapset);
- /* Allocate input buffer */
- in_data_type[8] = G_raster_map_type(b8, mapset);
- if( (infd[8] = G_open_cell_old(b8,mapset)) < 0){
- G_fatal_error(_("Cannot open cell file [%s]"), b8);
- }
- if( (G_get_cellhd(b8,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b8);
- }
- inrast[8] = G_allocate_raster_buf(in_data_type[8]);
+ /*Band8 */
+ /* find map in mapset */
+ mapset = G_find_cell2(b8, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), b8);
+ }
+ if (G_legal_filename(b8) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), b8);
+ }
+ infd[8] = G_open_cell_old(b8, mapset);
+ /* Allocate input buffer */
+ in_data_type[8] = G_raster_map_type(b8, mapset);
+ if ((infd[8] = G_open_cell_old(b8, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), b8);
+ }
+ if ((G_get_cellhd(b8, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), b8);
+ }
+ inrast[8] = G_allocate_raster_buf(in_data_type[8]);
+
/***************************************************/
+
/***************************************************/
- /* Allocate output buffer, use input map data_type */
- nrows = G_window_rows();
- ncols = G_window_cols();
- out_data_type=DCELL_TYPE;
- for(i=0;i<MAXFILES;i++){
- outrast[i] = G_allocate_raster_buf(out_data_type);
+ /* Allocate output buffer, use input map data_type */
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ out_data_type = DCELL_TYPE;
+ for (i = 0; i < MAXFILES; i++) {
+ outrast[i] = G_allocate_raster_buf(out_data_type);
+ }
+ if ((outfd[0] = G_open_raster_new(result1, 1)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result1);
+ if ((outfd[1] = G_open_raster_new(result2, 1)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result2);
+ if ((outfd[2] = G_open_raster_new(result3, 1)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result3);
+ if ((outfd[3] = G_open_raster_new(result4, 1)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result4);
+ if ((outfd[4] = G_open_raster_new(result5, 1)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result5);
+ if ((outfd[5] = G_open_raster_new(result61, 1)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result61);
+ if ((outfd[6] = G_open_raster_new(result62, 1)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result62);
+ if ((outfd[7] = G_open_raster_new(result7, 1)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result7);
+ if ((outfd[8] = G_open_raster_new(result8, 1)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result8);
+ /* Process pixels */
+ DCELL dout[MAXFILES];
+
+ DCELL d[MAXFILES];
+
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 2);
+ /* read input map */
+ for (i = 0; i < MAXFILES; i++) {
+ if ((G_get_raster_row(infd[i], inrast[i], row, in_data_type[i])) <
+ 0) {
+ G_fatal_error(_("Could not read from <%s>"), i);
+ }
}
- if ( (outfd[0] = G_open_raster_new (result1,1)) < 0)
- G_fatal_error (_("Could not open <%s>"),result1);
- if ( (outfd[1] = G_open_raster_new (result2,1)) < 0)
- G_fatal_error (_("Could not open <%s>"),result2);
- if ( (outfd[2] = G_open_raster_new (result3,1)) < 0)
- G_fatal_error (_("Could not open <%s>"),result3);
- if ( (outfd[3] = G_open_raster_new (result4,1)) < 0)
- G_fatal_error (_("Could not open <%s>"),result4);
- if ( (outfd[4] = G_open_raster_new (result5,1)) < 0)
- G_fatal_error (_("Could not open <%s>"),result5);
- if ( (outfd[5] = G_open_raster_new (result61,1)) < 0)
- G_fatal_error (_("Could not open <%s>"),result61);
- if ( (outfd[6] = G_open_raster_new (result62,1)) < 0)
- G_fatal_error (_("Could not open <%s>"),result62);
- if ( (outfd[7] = G_open_raster_new (result7,1)) < 0)
- G_fatal_error (_("Could not open <%s>"),result7);
- if ( (outfd[8] = G_open_raster_new (result8,1)) < 0)
- G_fatal_error (_("Could not open <%s>"),result8);
- /* Process pixels */
- DCELL dout[MAXFILES];
- DCELL d[MAXFILES];
- for (row = 0; row < nrows; row++)
- {
- G_percent (row, nrows, 2);
- /* read input map */
- for (i=0;i<MAXFILES;i++)
- {
- if((G_get_raster_row(infd[i],inrast[i],row,in_data_type[i])) < 0){
- G_fatal_error (_("Could not read from <%s>"),i);
- }
+ /*process the data */
+ for (col = 0; col < ncols; col++) {
+ for (i = 0; i < MAXFILES; i++) {
+ d[i] = (double)((CELL *) inrast[i])[col];
+ dout[i] =
+ dn2rad_landsat7(lmin[i], lmax[i], qcalmax[i], qcalmin[i],
+ d[i]);
+ if (i == 5 || i == 6) { /*if band 61/62, process brightness temperature */
+ if (dout[i] <= 0.0) {
+ G_set_d_null_value(&outrast[i][col], 1);
+ }
+ else {
+ dout[i] = tempk_landsat7(dout[i]);
+ }
}
- /*process the data */
- for (col=0; col < ncols; col++)
- {
- for(i=0;i<MAXFILES;i++)
- {
- d[i] = (double) ((CELL *) inrast[i])[col];
- dout[i]=dn2rad_landsat7(lmin[i],lmax[i],qcalmax[i],qcalmin[i],d[i]);
- if(i==5||i==6){/*if band 61/62, process brightness temperature*/
- if(dout[i]<=0.0){
- G_set_d_null_value(&outrast[i][col],1);
- }else{
- dout[i]=tempk_landsat7(dout[i]);
- }
- }else{/*process reflectance*/
- dout[i]=rad2ref_landsat7(dout[i],doy,sun_elevation,kexo[i]);
- }
- outrast[i][col] = dout[i];
- }
+ else { /*process reflectance */
+ dout[i] =
+ rad2ref_landsat7(dout[i], doy, sun_elevation,
+ kexo[i]);
}
- for(i=0;i<MAXFILES;i++){
- if (G_put_raster_row (outfd[i], outrast[i], out_data_type) < 0)
- G_fatal_error (_("Cannot write to <%s.%i>"),result,i+1);
- }
+ outrast[i][col] = dout[i];
+ }
}
- for (i=0;i<MAXFILES;i++)
- {
- G_free (inrast[i]);
- G_close_cell (infd[i]);
- G_free (outrast[i]);
- G_close_cell (outfd[i]);
+ for (i = 0; i < MAXFILES; i++) {
+ if (G_put_raster_row(outfd[i], outrast[i], out_data_type) < 0)
+ G_fatal_error(_("Cannot write to <%s.%i>"), result, i + 1);
}
- return 0;
+ }
+ for (i = 0; i < MAXFILES; i++) {
+ G_free(inrast[i]);
+ G_close_cell(infd[i]);
+ G_free(outrast[i]);
+ G_close_cell(outfd[i]);
+ }
+ return 0;
}
Modified: grass-addons/gipe/i.dn2full.l7/rad2ref_landsat7.c
===================================================================
--- grass-addons/gipe/i.dn2full.l7/rad2ref_landsat7.c 2008-08-21 23:06:30 UTC (rev 32985)
+++ grass-addons/gipe/i.dn2full.l7/rad2ref_landsat7.c 2008-08-21 23:12:57 UTC (rev 32986)
@@ -1,18 +1,20 @@
-#include<stdio.h>
-#include<math.h>
-#include<stdlib.h>
-
-#define PI 3.1415926
-
-/* Conversion of Radiance to Reflectance for Landsat 7ETM+ */
-
-double rad2ref_landsat7( double radiance, double doy,double sun_elevation, double k_exo )
-{
- double result, ds;
-
- ds = (1+0.01672*sin(2*PI*(doy-93.5)/365));
- result = (radiance/((cos((90-sun_elevation)*PI/180)/(PI*ds*ds))*k_exo));
-
- return result;
-}
-
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+
+#define PI 3.1415926
+
+ /* Conversion of Radiance to Reflectance for Landsat 7ETM+ */
+double rad2ref_landsat7(double radiance, double doy, double sun_elevation,
+ double k_exo)
+{
+ double result, ds;
+
+ ds = (1 + 0.01672 * sin(2 * PI * (doy - 93.5) / 365));
+ result =
+ (radiance /
+ ((cos((90 - sun_elevation) * PI / 180) / (PI * ds * ds)) * k_exo));
+ return result;
+}
+
+
Modified: grass-addons/gipe/i.dn2full.l7/tempk_landsat7.c
===================================================================
--- grass-addons/gipe/i.dn2full.l7/tempk_landsat7.c 2008-08-21 23:06:30 UTC (rev 32985)
+++ grass-addons/gipe/i.dn2full.l7/tempk_landsat7.c 2008-08-21 23:12:57 UTC (rev 32986)
@@ -1,15 +1,14 @@
-#include<stdio.h>
-#include<math.h>
-#include<stdlib.h>
-
-/* Surface temperature for Landsat 7ETM+ */
-
-double tempk_landsat7( double l6 )
-{
- double result;
-
- result = 1282.71 / (log ((666.09 / (l6)) + 1.0)) ;
-
- return result;
-}
-
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+
+ /* Surface temperature for Landsat 7ETM+ */
+double tempk_landsat7(double l6)
+{
+ double result;
+
+ result = 1282.71 / (log((666.09 / (l6)) + 1.0));
+ return result;
+}
+
+
More information about the grass-commit
mailing list