[GRASS-SVN] r30855 - in grass-addons/gipe: r.out.vic_met
r.out.vic_soil
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Apr 3 23:14:38 EDT 2008
Author: ychemin
Date: 2008-04-03 23:14:38 -0400 (Thu, 03 Apr 2008)
New Revision: 30855
Modified:
grass-addons/gipe/r.out.vic_met/main.c
grass-addons/gipe/r.out.vic_soil/main.c
Log:
Update
Modified: grass-addons/gipe/r.out.vic_met/main.c
===================================================================
--- grass-addons/gipe/r.out.vic_met/main.c 2008-04-03 15:15:59 UTC (rev 30854)
+++ grass-addons/gipe/r.out.vic_met/main.c 2008-04-04 03:14:38 UTC (rev 30855)
@@ -35,7 +35,8 @@
int verbose=1;
int not_ll=0;//if proj is not lat/long, it will be 1.
struct GModule *module;
- struct Option *input1, *input2, *input3, *output1;
+ struct Option *input1, *input2, *input3, *input4;
+ struct Option *output1, *output2;
struct Flag *flag1;
struct History history; //metadata
@@ -48,9 +49,11 @@
char *tmax_name;// input second time_series raster files names
char *tmin_name;// input third time_series raster files names
char *result1; //output file base name
+ char *result2; //output Soil file name
char result_lat_long[50]; //output file name
//File Descriptors
int infd_prcp[MAXFILES], infd_tmax[MAXFILES], infd_tmin[MAXFILES];
+ int infd;
int i=0,j=0;
double xp, yp;
@@ -62,9 +65,11 @@
void *inrast_prcp[MAXFILES];
void *inrast_tmax[MAXFILES];
void *inrast_tmin[MAXFILES];
+ void *inrast_elevation;
RASTER_MAP_TYPE data_type_inrast_prcp[MAXFILES];
RASTER_MAP_TYPE data_type_inrast_tmax[MAXFILES];
RASTER_MAP_TYPE data_type_inrast_tmin[MAXFILES];
+ RASTER_MAP_TYPE data_type_inrast_elevation;
FILE *f; // output ascii file
int grid_count; // grid cell count
@@ -77,13 +82,21 @@
int nfiles1, nfiles2, nfiles3, nfiles_shortest;// count no. input files
char **test1, **test2, **test3; // test number of input files
char **ptr1, **ptr2, **ptr3; // test number of input files
+ /****************************************/
+ /* Elevation map */
+ char *in;
+ FILE *g; // output soil ascii file
+ int process; // process grid cell switch
+ char *dummy_data1; // dummy data part 1
+ char *dummy_data2; // dummy data part 2
+ double elevation; // average evevation of grid cell
/************************************/
+
G_gisinit(argv[0]);
module = G_define_module();
- module->keywords = _("VIC, hydrology, precipitation, Tmax, Tmin");
- module->description = _("creates a meteorological ascii file \
- from 3 time series maps: precipitation, Tmax and Tmin.");
+ module->keywords = _("VIC, hydrology, precipitation, Tmax, Tmin, soil");
+ module->description = _("creates a meteorological ascii file for each non-null cell from 3 time series maps: precipitation, Tmax and Tmin. Also creates a dummy soil file from elevation.");
/* Define the different options */
input1 = G_define_standard_option(G_OPT_R_INPUT) ;
@@ -101,20 +114,30 @@
input3->multiple = YES;
input3->description=_("Names of the Tmin input maps");
- output1 = G_define_option() ;
+ input4 = G_define_standard_option(G_OPT_R_INPUT) ;
+ input4->key = _("dem");
+ input4->description=_("Name of the elevation input map");
+
+ output1 = G_define_option();
output1->key =_("output");
output1->description =_("Base Name of the output vic meteorological ascii files");
output1->answer =_("data_");
+ output2 = G_define_option() ;
+ output2->key =_("output_soil");
+ output2->description=_("Name of the output vic soil ascii file");
+ output2->answer =_("vic_soil.asc");
+
flag1 = G_define_flag() ;
flag1->key = 'a';
- flag1->description =_("append data if file already exists \
- (useful if adding additional year of data)");
+ flag1->description =_("append meteorological data if file already exists (useful if adding additional year of data)");
/********************/
if (G_parser(argc, argv))
exit (EXIT_FAILURE);
- result1 = output1->answer;
+ result1 = output1->answer;
+ in = input4->answer;
+ result2 = output2->answer;
/************************************************/
/* LOADING TMAX TIME SERIES MAPS */
test1 = input1->answers;
@@ -195,6 +218,18 @@
} else {
nfiles_shortest = nfiles3;
}
+ /************************************************/
+ /* Load Elevation file */
+ mapset = G_find_cell2(in, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), in);
+ }
+ data_type_inrast_elevation = G_raster_map_type(in,mapset);
+ if ( (infd = G_open_cell_old (in,mapset)) < 0)
+ G_fatal_error (_("Cannot open cell file [%s]"), in);
+ if (G_get_cellhd (in, mapset, &cellhd) < 0)
+ G_fatal_error (_("Cannot read file header of [%s])"), in);
+ inrast_elevation = G_allocate_raster_buf(data_type_inrast_elevation);
/***************************************************/
G_debug(3, "number of rows %d",cellhd.rows);
@@ -229,14 +264,27 @@
if ((oproj.pj = pj_latlong_from_proj(iproj.pj)) == NULL)
G_fatal_error(_("Unable to set up lat/long projection parameters"));
}//End of stolen from r.sun :P
-
+
+ /*Initialize grid cell process switch*/
+ g=fopen(result2,"w");
+ /*Initialize grid cell process switch*/
+ /*If =1 then process, if =0 then skip*/
+ process = 1;
+
/*Initialize grid cell count*/
grid_count = 1;
-
+ /*Initialize output file */
+ fprintf(g,"#RUN\tGRID\tLAT\tLNG\tINFILT\tDs\tDs_MAX\tWs\tC\tEXPT_1\tEXPT_2\tEXPT_3\tKsat_1\tKsat_2\tKsat_3\tPHI_1\tPHI_2\tPHI_3\tMOIST_1\tMOIST_2\tMOIST_3\tELEV\tDEPTH_1\tDEPTH_2\tDEPTH_3\tAVG_T\tDP\tBUBLE1\tBUBLE2\tBUBLE3\tQUARZ1\tQUARZ2\tQUARZ3\tBULKDN1\tBULKDN2\tBULKDN3\tPARTDN1\tPARTDN2\tPARTDN3\tOFF_GMT\tWcrFT1\tWcrFT2\tWcrFT3\tWpFT1\tWpFT2\tWpFT3\tZ0_SOIL\tZ0_SNOW\tPRCP\tRESM1\tRESM2\tRESM3\tFS_ACTV\tJULY_TAVG\n");
+
+ /*Initialize dummy data*/
+ dummy_data1 = "0.010\t1.e-4\t3.05\t0.93\t2\t4.0\t4.0\t4.0\t250.0\t250.0\t250.0\t-999\t-999\t-999\t0.4\t0.4\t0.4";
+ dummy_data2 = "0.1\t6.90\t2.000\t14.0\t4.0\t75.0\t75.0\t75.0\t0.24\t0.24\t0.24\t1306\t1367\t1367\t2650\t2650\t2650\t-6\t0.330\t0.330\t0.330\t0.133\t0.133\t0.133\t0.001\t0.010\t500\t0.02\t0.02\t0.02\t1\t18.665\n";
+
for (row = 0; row < nrows; row++){
DCELL d_prcp[MAXFILES];
DCELL d_tmax[MAXFILES];
DCELL d_tmin[MAXFILES];
+ DCELL d_elevation;
G_percent(row,nrows,2);
for(i=0;i<nfiles_shortest;i++){
if(G_get_raster_row(infd_prcp[i],inrast_prcp[i],row,data_type_inrast_prcp[i])<0)
@@ -246,6 +294,9 @@
if(G_get_raster_row(infd_tmin[i],inrast_tmin[i],row,data_type_inrast_tmin[i])<0)
G_fatal_error(_("Could not read from tmin<%d>"),i+1);
}
+ if(G_get_raster_row(infd,inrast_elevation,row,data_type_inrast_elevation)<0)
+ G_fatal_error(_("Could not read from <%s>"),in);
+
for (col=0; col < ncols; col++){
for(i=0;i<nfiles_shortest;i++){
/*Extract prcp time series data*/
@@ -285,7 +336,18 @@
break;
}
}
-
+ /*Extract elevation data*/
+ switch(data_type_inrast_elevation){
+ case CELL_TYPE:
+ d_elevation = (double) ((CELL *) inrast_elevation)[col];
+ break;
+ case FCELL_TYPE:
+ d_elevation = (double) ((FCELL *) inrast_elevation)[col];
+ break;
+ case DCELL_TYPE:
+ d_elevation = ((DCELL *) inrast_elevation)[col];
+ break;
+ }
/*Extract lat/long data*/
latitude = ymax - ( row * stepy );
longitude = xmin + ( col * stepx );
@@ -293,14 +355,17 @@
if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0) {
G_fatal_error(_("Error in pj_do_proj"));
}
+ }else{
+ //Do nothing
}
if(G_is_d_null_value(&prcp[0])||
G_is_d_null_value(&d_tmax[0])||
- G_is_d_null_value(&d_tmin[0])){
+ G_is_d_null_value(&d_tmin[0])||
+ G_is_d_null_value(&d_elevation)){
/* Do nothing */
} else {
/* Make the output .dat file name */
- sprintf(result_lat_long,"%s%.4f%s%.4f%s",result1,latitude,"_",longitude,".dat");
+ sprintf(result_lat_long,"%s%.4f%s%.4f",result1,latitude,"_",longitude);
/*Open new ascii file*/
if (flag1->answer){
/*Initialize grid cell in append mode*/
@@ -314,11 +379,14 @@
fprintf(f,"%.2f %.2f %.2f\n", d_prcp[i], d_tmax[i], d_tmin[i]);
}
fclose(f);
+ /*Print to soil ascii file*/
+ fprintf(g,"%d\t%d\t%6.3f\t%7.3f\t%s\t%7.2f\t%s\n", process, grid_count, latitude, longitude, dummy_data1, d_elevation, dummy_data2);
grid_count=grid_count+1;
}
}
}
G_message(_("Created %d VIC meteorological files"),grid_count);
+ G_message(_("Created %d VIC grid cells soil definitions"),grid_count);
for(i=0;i<nfiles1;i++){
G_free (inrast_prcp[i]);
G_close_cell (infd_prcp[i]);
@@ -331,6 +399,9 @@
G_free (inrast_tmin[i]);
G_close_cell (infd_tmin[i]);
}
+ G_free (inrast_elevation);
+ G_close_cell (infd);
+ fclose(g);
exit(EXIT_SUCCESS);
}
Modified: grass-addons/gipe/r.out.vic_soil/main.c
===================================================================
--- grass-addons/gipe/r.out.vic_soil/main.c 2008-04-03 15:15:59 UTC (rev 30854)
+++ grass-addons/gipe/r.out.vic_soil/main.c 2008-04-04 03:14:38 UTC (rev 30855)
@@ -135,6 +135,7 @@
/* Set up parameters for projection to lat/long if necessary */
if ((G_projection() != PROJECTION_LL)) {
not_ll=1;
+ G_message("projection is not Lat/long, converting...");
struct Key_Value *in_proj_info, *in_unit_info;
if ((in_proj_info = G_get_projinfo()) == NULL)
G_fatal_error(_("Can't get projection info of current location"));
@@ -162,7 +163,7 @@
{
/*Extract lat/long data*/
latitude = ymax - ( row * stepy );
- longitude = xmin + ( col * stepx );
+ longitude = xmin + ( col * stepx ) ;
if(not_ll){
if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0) {
G_fatal_error(_("Error in pj_do_proj"));
@@ -179,14 +180,19 @@
d_elevation = (double) ((FCELL *) inrast_elevation)[col];
break;
case DCELL_TYPE:
- d_elevation = ((DCELL *) inrast_elevation)[col];
+ d_elevation = (double) ((DCELL *) inrast_elevation)[col];
break;
}
- /*Print to ascii file*/
- fprintf(f,"%d\t%d\t%6.3f\t%7.3f\t%s\t%7.2f\t%s\n", process, grid_count, latitude, longitude, dummy_data1, d_elevation, dummy_data2);
- grid_count=grid_count+1;
+ if(G_is_d_null_value(&d_elevation)){
+ /* Do nothing */
+ } else {
+ /*Print to ascii file*/
+ fprintf(f,"%d\t%d\t%6.3f\t%7.3f\t%s\t%7.2f\t%s\n", process, grid_count, latitude, longitude, dummy_data1, d_elevation, dummy_data2);
+ grid_count=grid_count+1;
+ }
}
}
+ G_message("Generated soil info for %d grid cells",grid_count);
G_free (inrast_elevation);
G_close_cell (infd);
fclose(f);
More information about the grass-commit
mailing list