[GRASS-SVN] r30683 - grass-addons/gipe/i.dn2full.l7
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Mar 21 23:23:01 EDT 2008
Author: ychemin
Date: 2008-03-21 23:23:01 -0400 (Fri, 21 Mar 2008)
New Revision: 30683
Added:
grass-addons/gipe/i.dn2full.l7/dn2rad_landsat5.c
grass-addons/gipe/i.dn2full.l7/l5inread.c
grass-addons/gipe/i.dn2full.l7/rad2ref_landsat5.c
grass-addons/gipe/i.dn2full.l7/tempk_landsat5.c
Removed:
grass-addons/gipe/i.dn2full.l7/dn2rad_landsat7.c
grass-addons/gipe/i.dn2full.l7/l7inread.c
Modified:
grass-addons/gipe/i.dn2full.l7/Makefile
grass-addons/gipe/i.dn2full.l7/main.c
Log:
Updated DN2rad of Landsat 5TM using NLAPS report file
Modified: grass-addons/gipe/i.dn2full.l7/Makefile
===================================================================
--- grass-addons/gipe/i.dn2full.l7/Makefile 2008-03-21 14:38:12 UTC (rev 30682)
+++ grass-addons/gipe/i.dn2full.l7/Makefile 2008-03-22 03:23:01 UTC (rev 30683)
@@ -1,6 +1,6 @@
MODULE_TOPDIR = ../..
-PGM = i.dn2full.l7
+PGM = i.dn2full.l5
LIBES = $(GISLIB)
DEPENDENCIES = $(GISDEP)
Added: grass-addons/gipe/i.dn2full.l7/dn2rad_landsat5.c
===================================================================
--- grass-addons/gipe/i.dn2full.l7/dn2rad_landsat5.c (rev 0)
+++ grass-addons/gipe/i.dn2full.l7/dn2rad_landsat5.c 2008-03-22 03:23:01 UTC (rev 30683)
@@ -0,0 +1,163 @@
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+
+// Conversion of DN to Radiance for Landsat 5TM
+// yann.chemin at gmail.com - Yann Chemin - GPL>=2, 2008.
+
+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;
+}
+
Deleted: grass-addons/gipe/i.dn2full.l7/dn2rad_landsat7.c
===================================================================
--- grass-addons/gipe/i.dn2full.l7/dn2rad_landsat7.c 2008-03-21 14:38:12 UTC (rev 30682)
+++ grass-addons/gipe/i.dn2full.l7/dn2rad_landsat7.c 2008-03-22 03:23:01 UTC (rev 30683)
@@ -1,19 +0,0 @@
-#include<stdio.h>
-#include<math.h>
-#include<stdlib.h>
-
-// Conversion of DN to Radiance for Landsat 7ETM+
-// http://ltpwww.gsfc.nasa.gov/IAS/handbook/handbook_htmls/chapter11/chapter11.html#section11.3
-// ychemin at gmail.com - Yann Chemin - LGPL, Copylefted, 2006.
-
-double dn2rad_landsat7( double Lmin, double LMax, double QCalMax, double QCalmin, int DN )
-{
- double result, gain, offset;
-
- gain = (LMax-Lmin)/(QCalMax-QCalmin);
- offset = Lmin;
- result = gain * ((double) DN - QCalmin) + offset ;
-
- return result;
-}
-
Added: grass-addons/gipe/i.dn2full.l7/l5inread.c
===================================================================
--- grass-addons/gipe/i.dn2full.l7/l5inread.c (rev 0)
+++ grass-addons/gipe/i.dn2full.l7/l5inread.c 2008-03-22 03:23:01 UTC (rev 30683)
@@ -0,0 +1,94 @@
+#include<stdio.h>
+#include<stdlib.h>
+#include<string.h>
+#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)
+{
+ FILE *f;
+ char s[1000], *ptr;
+ 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;
+}
Deleted: grass-addons/gipe/i.dn2full.l7/l7inread.c
===================================================================
--- grass-addons/gipe/i.dn2full.l7/l7inread.c 2008-03-21 14:38:12 UTC (rev 30682)
+++ grass-addons/gipe/i.dn2full.l7/l7inread.c 2008-03-22 03:23:01 UTC (rev 30683)
@@ -1,399 +0,0 @@
-#include<stdio.h>
-#include<stdlib.h>
-#include<string.h>
-#include <grass/gis.h>
-#include <grass/glocale.h>
-
-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)
-{
- FILE *f;
- char s[1000], *ptr;
- int i=0;
- char *p;
-
- char band1[80],band2[80],band3[80];
- char band4[80],band5[80],band61[80];
- char band62[80],band7[80],band8[80];//Load .tif names
-
- if((f=fopen(metfName,"r"))==NULL){
- G_fatal_error (_(".met file [%s] not found, are you in the proper directory?"), metfName);
- }
-
- while (fgets(s,1000,f)!=NULL)
- {
-// printf("2%s\n",s);
- ptr = strstr(s, "ACQUISITION_DATE");
- if (ptr != NULL)
- {
-// printf("2\t");
- p = strtok(ptr, " =");
- p = strtok(NULL, " =-");
- *year = atoi(p);
-// printf("3\t");
- p = strtok(NULL, "-");
- *month = atoi(p);
-// printf("4\t");
- p = strtok(NULL, "-");
- *day = atoi(p);
-// printf("5\n");
- }
- ptr = strstr(s, "BAND1_FILE_NAME");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =\"");
- // printf("l7inread.c:p=%s\n",p);
- snprintf(band1, 80, "%s", p);
- // printf("l7inread.c:band1=%s\n",band1);
- }
- ptr = strstr(s, "BAND2_FILE_NAME");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =\"");
- snprintf(band2, 80, "%s", p);
- }
- ptr = strstr(s, "BAND3_FILE_NAME");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =\"");
- snprintf(band3, 80, "%s", p);
- }
- ptr = strstr(s, "BAND4_FILE_NAME");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =\"");
- snprintf(band4, 80, "%s", p);
- }
- ptr = strstr(s, "BAND5_FILE_NAME");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =\"");
- snprintf(band5, 80, "%s", p);
- }
- ptr = strstr(s, "BAND61_FILE_NAME");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =\"");
- snprintf(band61, 80, "%s", p);
- }
- ptr = strstr(s, "BAND62_FILE_NAME");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =\"");
- snprintf(band62, 80, "%s", p);
- }
- ptr = strstr(s, "BAND7_FILE_NAME");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =\"");
- snprintf(band7, 80, "%s", p);
- }
- ptr = strstr(s, "BAND8_FILE_NAME");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =\"");
- snprintf(band8, 80, "%s", p);
- }
- ptr = strstr(s, "SUN_AZIMUTH");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- *sun_azimuth = atof(p);
-
- }
- ptr = strstr(s, "SUN_ELEVATION");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- *sun_elevation = atof(p);
-
- }
- ptr = strstr(s, "\tLMAX_BAND1");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmax[0] = atof(p);
- }
- ptr = strstr(s, "\tLMIN_BAND1");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmin[0] = atof(p);
- }
- ptr = strstr(s, "\tLMAX_BAND2");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmax[1] = atof(p);
- }
- ptr = strstr(s, "\tLMIN_BAND2");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmin[1] = atof(p);
- }
- ptr = strstr(s, "\tLMAX_BAND3");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmax[2] = atof(p);
- }
- ptr = strstr(s, "\tLMIN_BAND3");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmin[2] = atof(p);
- }
- ptr = strstr(s, "\tLMAX_BAND4");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmax[3] = atof(p);
- }
- ptr = strstr(s, "\tLMIN_BAND4");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmin[3] = atof(p);
- }
- ptr = strstr(s, "\tLMAX_BAND5");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmax[4] = atof(p);
- }
- ptr = strstr(s, "\tLMIN_BAND5");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmin[4] = atof(p);
- }
- ptr = strstr(s, "\tLMAX_BAND61");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmax[5] = atof(p);
- }
- ptr = strstr(s, "\tLMIN_BAND61");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmin[5] = atof(p);
- }
- ptr = strstr(s, "\tLMAX_BAND62");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmax[6] = atof(p);
- }
- ptr = strstr(s, "\tLMIN_BAND62");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmin[6] = atof(p);
- }
- ptr = strstr(s, "\tLMAX_BAND7");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmax[7] = atof(p);
- }
- ptr = strstr(s, "\tLMIN_BAND7");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmin[7] = atof(p);
- }
- ptr = strstr(s, "\tLMAX_BAND8");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmax[8] = atof(p);
- }
- ptr = strstr(s, "\tLMIN_BAND8");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- lmin[8] = atof(p);
- }
- ptr = strstr(s, "QCALMAX_BAND1");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmax[0] = atof(p);
- }
- ptr = strstr(s, "QCALMIN_BAND1");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmin[0] = atof(p);
- }
- ptr = strstr(s, "QCALMAX_BAND2");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmax[1] = atof(p);
- }
- ptr = strstr(s, "QCALMIN_BAND2");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmin[1] = atof(p);
- }
- ptr = strstr(s, "QCALMAX_BAND3");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmax[2] = atof(p);
- }
- ptr = strstr(s, "QCALMIN_BAND3");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmin[2] = atof(p);
- }
- ptr = strstr(s, "QCALMAX_BAND4");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmax[3] = atof(p);
- }
- ptr = strstr(s, "QCALMIN_BAND4");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmin[3] = atof(p);
- }
- ptr = strstr(s, "QCALMAX_BAND5");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmax[4] = atof(p);
- }
- ptr = strstr(s, "QCALMIN_BAND5");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmin[4] = atof(p);
- }
- ptr = strstr(s, "QCALMAX_BAND61");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmax[5] = atof(p);
- }
- ptr = strstr(s, "QCALMIN_BAND61");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmin[5] = atof(p);
- }
- ptr = strstr(s, "QCALMAX_BAND62");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmax[6] = atof(p);
- }
- ptr = strstr(s, "QCALMIN_BAND62");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmin[6] = atof(p);
- }
- ptr = strstr(s, "QCALMAX_BAND7");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmax[7] = atof(p);
- }
- ptr = strstr(s, "QCALMIN_BAND7");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmin[7] = atof(p);
- }
- ptr = strstr(s, "QCALMAX_BAND8");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmax[8] = atof(p);
- }
- ptr = strstr(s, "QCALMIN_BAND8");
- if (ptr != NULL)
- {
- p = strtok(ptr, " =");
- p = strtok(NULL, " =");
- qcalmin[8] = atof(p);
- }
- }
-
- snprintf(b1, 80, "%s", band1);
- snprintf(b2, 80, "%s", band2);
- snprintf(b3, 80, "%s", band3);
- snprintf(b4, 80, "%s", band4);
- snprintf(b5, 80, "%s", band5);
- snprintf(b61, 80, "%s", band61);
- snprintf(b62, 80, "%s", band62);
- snprintf(b7, 80, "%s", band7);
- snprintf(b8, 80, "%s", band8);
-// 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);
-
-// for (i=0;i<MAXFILES;i++){
-// printf("lmin[%i]=%f\n",i,lmin[i]);
-// printf("lmax[%i]=%f\n",i,lmax[i]);
-// printf("qcalmin[%i]=%f\n",i,qcalmin[i]);
-// printf("qcalmax[%i]=%f\n",i,qcalmax[i]);
-// }
- (void)fclose(f);
- return;
-}
Modified: grass-addons/gipe/i.dn2full.l7/main.c
===================================================================
--- grass-addons/gipe/i.dn2full.l7/main.c 2008-03-21 14:38:12 UTC (rev 30682)
+++ grass-addons/gipe/i.dn2full.l7/main.c 2008-03-22 03:23:01 UTC (rev 30683)
@@ -1,8 +1,8 @@
/****************************************************************************
*
- * MODULE: i.dn2full.l7
+ * MODULE: i.dn2full.l5
* AUTHOR(S): Yann Chemin - ychemin at gmail.com
- * PURPOSE: Calculate TOA Reflectance&Temperature for Landsat7 from DN.
+ * PURPOSE: Calculate TOA Reflectance&Temperature for Landsat5 from DN.
*
* COPYRIGHT: (C) 2002 by the GRASS Development Team
*
@@ -20,7 +20,7 @@
#include <grass/gis.h>
#include <grass/glocale.h>
-#define MAXFILES 9
+#define MAXFILES 7
//sun exo-atmospheric irradiance
#define KEXO1 1969.0
@@ -28,17 +28,16 @@
#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 KEXO6 82.07
#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 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_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 );
+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[])
@@ -48,7 +47,6 @@
int nrows, ncols;
int row,col;
- int verbose=1;
struct GModule *module;
struct Option *input,*output;
@@ -62,7 +60,7 @@
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] ;
+ char result5[80], result6[80], result7[80];
//File Descriptors
int infd[MAXFILES];
@@ -80,19 +78,21 @@
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];
+ 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]);
@@ -100,56 +100,59 @@
module = G_define_module();
module->keywords = _("DN, reflectance, temperature, import");
module->description =
- _("Calculates Top of Atmosphere Reflectance/Temperature from Landsat 7 DN.\n");
+ _("Calculates Top of Atmosphere Reflectance/Temperature from Landsat 5 DN.\n");
/* Define the different options */
input = G_define_option() ;
- input->key = _("metfile");
+ input->key = _("file");
input->type = TYPE_STRING;
input->required = YES;
input->gisprompt = _("old_file,file,file");
- input->description= _("Landsat 7ETM+ Header File (.met)");
+ input->description= _("Landsat 5TM NLAPS processing report File (.txt)");
- output = G_define_option() ;
+ output = G_define_standard_option(G_OPT_R_OUTPUT) ;
output->key = _("output");
- output->type = TYPE_STRING;
- output->required = YES;
- output->gisprompt = _("new,cell,raster");
output->description= _("Base name of the output layers (will add .x)");
-
- /* Define the different flags */
-
- flag1 = G_define_flag() ;
- flag1->key =_('q');
- flag1->description =_("Quiet");
-
/********************/
if (G_parser(argc, argv))
exit (-1);
metfName = input->answer;
result = output->answer;
-
- verbose = (!flag1->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);
+ l5_in_read(metfName,&l5path,&l5row,&latitude,&longitude,&sun_elevation,&sun_azimuth,&c_year,&c_month,&c_day,&day,&month,&year,&decimal_hour);
+ /********************/
+ //Prepare the input file names
+ /********************/
+ doy = date2doy(day,month,year);
+ printf("doy=%i\n",doy);
+ if(year<2000){
+ temp = year - 1900;
+ } else {
+ temp = year - 2000;
+ }
+ if (temp >=10){
+ snprintf(b1, 80, "%s%s%s%s%s%s%s%s","LT5",l5path,"0",l5row,"00",temp,doy,"50_B1.tif");
+ snprintf(b2, 80, "%s%s%s%s%s%s%s%s","LT5",l5path,"0",l5row,"00",temp,doy,"50_B2.tif");
+ snprintf(b3, 80, "%s%s%s%s%s%s%s%s","LT5",l5path,"0",l5row,"00",temp,doy,"50_B3.tif");
+ snprintf(b4, 80, "%s%s%s%s%s%s%s%s","LT5",l5path,"0",l5row,"00",temp,doy,"50_B4.tif");
+ snprintf(b5, 80, "%s%s%s%s%s%s%s%s","LT5",l5path,"0",l5row,"00",temp,doy,"50_B5.tif");
+ snprintf(b6, 80, "%s%s%s%s%s%s%s%s","LT5",l5path,"0",l5row,"00",temp,doy,"50_B6.tif");
+ snprintf(b7, 80, "%s%s%s%s%s%s%s%s","LT5",l5path,"0",l5row,"00",temp,doy,"50_B7.tif");
+ } else {
+ snprintf(b1, 80, "%s%s%s%s%s%s%s%s","LT5",l5path,"0",l5row,"000",temp,doy,"50_B1.tif");
+ snprintf(b2, 80, "%s%s%s%s%s%s%s%s","LT5",l5path,"0",l5row,"000",temp,doy,"50_B2.tif");
+ snprintf(b3, 80, "%s%s%s%s%s%s%s%s","LT5",l5path,"0",l5row,"000",temp,doy,"50_B3.tif");
+ snprintf(b4, 80, "%s%s%s%s%s%s%s%s","LT5",l5path,"0",l5row,"000",temp,doy,"50_B4.tif");
+ snprintf(b5, 80, "%s%s%s%s%s%s%s%s","LT5",l5path,"0",l5row,"000",temp,doy,"50_B5.tif");
+ snprintf(b6, 80, "%s%s%s%s%s%s%s%s","LT5",l5path,"0",l5row,"000",temp,doy,"50_B6.tif");
+ snprintf(b7, 80, "%s%s%s%s%s%s%s%s","LT5",l5path,"0",l5row,"000",temp,doy,"50_B7.tif");
+ }
// printf("%f/%f/%i-%i-%i\n",sun_elevation,sun_azimuth,day,month,year);
// for(i=0;i<MAXFILES;i++){
// printf("%i=>%f, %f, %f, %f\n",i,lmin[i],lmax[i],qcalmin[i],qcalmax[i]);
// }
-// doy = date2doy(day,month,year);
-// printf("doy=%i\n",doy);
- /********************/
-// printf("b1=%s\n",b1);
-// printf("b2=%s\n",b2);
-// printf("b3=%s\n",b3);
-// printf("b4=%s\n",b4);
-// printf("b5=%s\n",b5);
-// printf("b61=%s\n",b61);
-// printf("b62=%s\n",b62);
-// printf("b7=%s\n",b7);
-// printf("b8=%s\n",b8);
// exit(1);
/********************/
//Prepare the ouput file names
@@ -166,14 +169,10 @@
// printf("%s\n",result4);
snprintf(result5, 80, "%s%s",result,".5");
// printf("%s\n",result5);
- snprintf(result61, 80, "%s%s",result,".61");
-// printf("%s\n",result61);
- snprintf(result62, 80, "%s%s",result,".62");
-// printf("%s\n",result62);
+ snprintf(result6, 80, "%s%s",result,".6");
+// printf("%s\n",result6);
snprintf(result7, 80, "%s%s",result,".7");
// printf("%s\n",result7);
- snprintf(result8, 80, "%s%s",result,".8");
-// printf("%s\n",result8);
/********************/
//Prepare sun exo-atm irradiance
@@ -184,12 +183,8 @@
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[5]=KEXO6;
- //******************************************
/***************************************************/
//Band1
/* find map in mapset */
@@ -296,48 +291,27 @@
inrast[4] = G_allocate_raster_buf(in_data_type[4]);
/***************************************************/
/***************************************************/
- //Band61
+ //Band6
/* find map in mapset */
- mapset = G_find_cell2 (b61, "");
+ mapset = G_find_cell2 (b6, "");
if (mapset == NULL){
- G_fatal_error (_("cell file [%s] not found"), b61);
+ G_fatal_error (_("cell file [%s] not found"), b6);
}
- if (G_legal_filename (b61) < 0){
- G_fatal_error (_("[%s] is an illegal name"), b61);
+ if (G_legal_filename (b6) < 0){
+ G_fatal_error (_("[%s] is an illegal name"), b6);
}
- infd[5] = G_open_cell_old (b61, mapset);
+ infd[5] = G_open_cell_old (b6, 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);
+ 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(b61,mapset,&cellhd)) < 0){
- G_fatal_error(_("Cannot read file header of [%s]"), b61);
+ 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]);
/***************************************************/
/***************************************************/
- //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, "");
@@ -347,39 +321,18 @@
if (G_legal_filename (b7) < 0){
G_fatal_error (_("[%s] is an illegal name"), b7);
}
- infd[7] = G_open_cell_old (b7, mapset);
+ infd[6] = 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){
+ 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[7] = G_allocate_raster_buf(in_data_type[7]);
+ inrast[6] = G_allocate_raster_buf(in_data_type[6]);
/***************************************************/
/***************************************************/
- //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();
@@ -397,22 +350,16 @@
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)
+ 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);
- 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++)
{
- if (verbose){
- G_percent (row, nrows, 2);
- }
+ G_percent (row, nrows, 2);
/* read input map */
for (i=0;i<MAXFILES;i++)
{
@@ -426,15 +373,15 @@
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
+ 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){
dout[i]=-999.990;
}else{
- dout[i]=tempk_landsat7(dout[i]);
+ dout[i]=tempk_landsat5(dout[i]);
}
}else{//process reflectance
- dout[i]=rad2ref_landsat7(dout[i],doy,sun_elevation,kexo[i]);
+ dout[i]=rad2ref_landsat5(dout[i],doy,sun_elevation,kexo[i]);
}
((DCELL *) outrast[i])[col] = dout[i];
}
Added: grass-addons/gipe/i.dn2full.l7/rad2ref_landsat5.c
===================================================================
--- grass-addons/gipe/i.dn2full.l7/rad2ref_landsat5.c (rev 0)
+++ grass-addons/gipe/i.dn2full.l7/rad2ref_landsat5.c 2008-03-22 03:23:01 UTC (rev 30683)
@@ -0,0 +1,19 @@
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+
+#define PI 3.1415926
+
+// Conversion of Radiance to Reflectance for Landsat 5TM
+// ychemin at gmail.com - Yann Chemin - GPL>=2, 2008.
+
+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;
+}
+
Added: grass-addons/gipe/i.dn2full.l7/tempk_landsat5.c
===================================================================
--- grass-addons/gipe/i.dn2full.l7/tempk_landsat5.c (rev 0)
+++ grass-addons/gipe/i.dn2full.l7/tempk_landsat5.c 2008-03-22 03:23:01 UTC (rev 30683)
@@ -0,0 +1,16 @@
+#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;
+}
+
More information about the grass-commit
mailing list