[GRASS-SVN] r73244 - in grass-addons/grass7/imagery: . i.rh
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Sep 4 01:37:38 PDT 2018
Author: ychemin
Date: 2018-09-04 01:37:38 -0700 (Tue, 04 Sep 2018)
New Revision: 73244
Added:
grass-addons/grass7/imagery/i.rh/
grass-addons/grass7/imagery/i.rh/Makefile
grass-addons/grass7/imagery/i.rh/i.rh.html
grass-addons/grass7/imagery/i.rh/main.c
grass-addons/grass7/imagery/i.rh/rh.c
Modified:
grass-addons/grass7/imagery/Makefile
Log:
added i.rh
Modified: grass-addons/grass7/imagery/Makefile
===================================================================
--- grass-addons/grass7/imagery/Makefile 2018-09-04 05:57:32 UTC (rev 73243)
+++ grass-addons/grass7/imagery/Makefile 2018-09-04 08:37:38 UTC (rev 73244)
@@ -21,6 +21,7 @@
i.modis \
i.points.auto \
i.pysptools.unmix \
+ i.rh \
i.rotate \
i.segment.gsoc \
i.segment.hierarchical \
Added: grass-addons/grass7/imagery/i.rh/Makefile
===================================================================
--- grass-addons/grass7/imagery/i.rh/Makefile (rev 0)
+++ grass-addons/grass7/imagery/i.rh/Makefile 2018-09-04 08:37:38 UTC (rev 73244)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.rh
+
+LIBES = $(GISLIB) $(RASTERLIB)
+DEPENDENCIES = $(GISDEP) $(RASTERDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/grass7/imagery/i.rh/i.rh.html
===================================================================
--- grass-addons/grass7/imagery/i.rh/i.rh.html (rev 0)
+++ grass-addons/grass7/imagery/i.rh/i.rh.html 2018-09-04 08:37:38 UTC (rev 73244)
@@ -0,0 +1,30 @@
+<h2>DESCRIPTION</h2>
+
+<em>i.rh</em> calculates the relative humidity
+
+-s will return saturated water vapour instead.
+
+-a will return actual water vapour instead.
+
+<h2>NOTES</h2>
+
+
+<h2>REFERENCES</h2>
+https://www.researchgate.net/publication/227247013_High-resolution_Surface_Relative_Humidity_Computation_Using_MODIS_Image_in_Peninsular_Malaysia/
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="i.vi.html">i.vi</a><br>
+<a href="i.albedo.html">i.albedo</a><br>
+<a href="i.wi.html">i.vi</a><br>
+</em>
+
+
+<h2>AUTHOR</h2>
+
+Yann Chemin, Plumergat, France<br>
+
+
+<p>
+<i>Last changed: $Date: 2017-01-19 15:39:35 +0000 (Thu, 19 Jan 2017) $</i>
Added: grass-addons/grass7/imagery/i.rh/main.c
===================================================================
--- grass-addons/grass7/imagery/i.rh/main.c (rev 0)
+++ grass-addons/grass7/imagery/i.rh/main.c 2018-09-04 08:37:38 UTC (rev 73244)
@@ -0,0 +1,223 @@
+
+/****************************************************************************
+ *
+ * MODULE: i.rh
+ * AUTHOR(S): Yann Chemin - yann.chemin at gmail.com
+ * PURPOSE: Calculates relative humidity
+ *
+ * COPYRIGHT: (C) 2017-2019 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+double rh(double PW,double Pa,double Ta,double dem);
+double esat(double tamean);
+double eact(double esat,double rh);
+
+int main(int argc, char *argv[])
+{
+ struct Cell_head cellhd; /*region+header info */
+ char *mapset=""; /*mapset name */
+ int nrows, ncols;
+ int row, col;
+ struct GModule *module;
+ struct Option *input1, *input2, *input3, *input4, *output1;
+ struct Flag *flag1, *flag2;
+ struct History history; /*metadata */
+ char *name; /*input raster name */
+ char *result1; /*output raster name */
+ int infd_pw, infd_pa, infd_ta, infd_dem;
+ int outfd1;
+ char *pw, *pa, *ta, *dem;
+ int i = 0, j = 0;
+ void *inrast_pw, *inrast_pa, *inrast_ta, *inrast_dem;
+ CELL * outrast1;
+ RASTER_MAP_TYPE data_type_output = CELL_TYPE;
+ RASTER_MAP_TYPE data_type_pw;
+ RASTER_MAP_TYPE data_type_pa;
+ RASTER_MAP_TYPE data_type_ta;
+ RASTER_MAP_TYPE data_type_dem;
+ /************************************/
+ G_gisinit(argv[0]);
+ module = G_define_module();
+ G_add_keyword(_("imagery"));
+ G_add_keyword(_("atmosphere"));
+ G_add_keyword(_("humidity"));
+ G_add_keyword(_("water vapour"));
+ G_add_keyword(_("precipitable water"));
+ module->description = _("Water in atmosphere: relative humidity, water vapour (saturated, actual)");
+
+ /* Define the different options */
+ input1 = G_define_standard_option(G_OPT_R_INPUT);
+ input1->key = "pw";
+ input1->description = _("Name of the Precipitable Water layer");
+ input2 = G_define_standard_option(G_OPT_R_INPUT);
+ input2->key = "pa";
+ input2->description = _("Name of the Atmospheric Pressure layer");
+ input3 = G_define_standard_option(G_OPT_R_INPUT);
+ input3->key = "ta";
+ input3->description = _("Name of the mean daily air temperature layer");
+ input4 = G_define_standard_option(G_OPT_R_INPUT);
+ input4->key = "dem";
+ input4->description = _("Name of the elevation layer");
+
+ output1 = G_define_standard_option(G_OPT_R_OUTPUT);
+ output1->description = _("Name of the output layer");
+
+ flag1 = G_define_flag();
+ flag1->key = 's';
+ flag1->description = _("Output the saturated water vapour (esat)");
+
+ flag2 = G_define_flag();
+ flag2->key = 'a';
+ flag2->description = _("Output the actual water vapour (eact)");
+
+ /********************/
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+ pw = input1->answer;
+ pa = input2->answer;
+ ta = input3->answer;
+ dem = input4->answer;
+ result1 = output1->answer;
+ if(flag1->answer && flag2->answer)
+ G_fatal_error(_("Cannot output both esat and eact, please choose"));
+ /***************************************************/
+ data_type_pw = Rast_map_type(pw, mapset);
+ if ((infd_pw = Rast_open_old(pw, mapset)) < 0)
+ G_fatal_error(_("Cannot open cell file [%s]"), pw);
+ Rast_get_cellhd(pw, mapset, &cellhd);
+ inrast_pw = Rast_allocate_buf(data_type_pw);
+ /***************************************************/
+ data_type_pa = Rast_map_type(pa, mapset);
+ if ((infd_pa = Rast_open_old(pa, mapset)) < 0)
+ G_fatal_error(_("Cannot open cell file [%s]"), pa);
+ Rast_get_cellhd(pa, mapset, &cellhd);
+ inrast_pa = Rast_allocate_buf(data_type_pa);
+ /***************************************************/
+ data_type_ta = Rast_map_type(ta, mapset);
+ if ((infd_ta = Rast_open_old(ta, mapset)) < 0)
+ G_fatal_error(_("Cannot open cell file [%s]"), ta);
+ Rast_get_cellhd(ta, mapset, &cellhd);
+ inrast_ta = Rast_allocate_buf(data_type_ta);
+ /***************************************************/
+ data_type_dem = Rast_map_type(dem, mapset);
+ if ((infd_dem = Rast_open_old(dem, mapset)) < 0)
+ G_fatal_error(_("Cannot open cell file [%s]"), dem);
+ Rast_get_cellhd(dem, mapset, &cellhd);
+ inrast_dem = Rast_allocate_buf(data_type_dem);
+ /***************************************************/
+ G_debug(3, "number of rows %d", cellhd.rows);
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+ outrast1 = Rast_allocate_buf(data_type_output);
+ /* Create New raster files */
+ if ((outfd1 = Rast_open_new(result1, data_type_output)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result1);
+ /* Process pixels */
+ for (row = 0; row < nrows; row++)
+ {
+ CELL d;
+ DCELL d_pw;
+ DCELL d_pa;
+ DCELL d_ta;
+ DCELL d_dem;
+ G_percent(row, nrows, 2);
+ /* read input maps */
+ Rast_get_row(infd_pw, inrast_pw, row, data_type_pw);
+ Rast_get_row(infd_pa, inrast_pa, row, data_type_pa);
+ Rast_get_row(infd_ta, inrast_ta, row, data_type_ta);
+ Rast_get_row(infd_dem, inrast_dem, row, data_type_dem);
+ /*process the data */
+ for (col = 0; col < ncols; col++)
+ {
+ switch (data_type_pw) {
+ case CELL_TYPE:
+ d_pw = (double)((CELL *) inrast_pw)[col];
+ break;
+ case FCELL_TYPE:
+ d_pw = (double)((FCELL *) inrast_pw)[col];
+ break;
+ case DCELL_TYPE:
+ d_pw = ((DCELL *) inrast_pw)[col];
+ break;
+ }
+ switch (data_type_pa) {
+ case CELL_TYPE:
+ d_pa = (double)((CELL *) inrast_pa)[col];
+ break;
+ case FCELL_TYPE:
+ d_pa = (double)((FCELL *) inrast_pa)[col];
+ break;
+ case DCELL_TYPE:
+ d_pa = ((DCELL *) inrast_pa)[col];
+ break;
+ }
+ switch (data_type_ta) {
+ case CELL_TYPE:
+ d_ta = (double)((CELL *) inrast_ta)[col];
+ break;
+ case FCELL_TYPE:
+ d_ta = (double)((FCELL *) inrast_ta)[col];
+ break;
+ case DCELL_TYPE:
+ d_ta = ((DCELL *) inrast_ta)[col];
+ break;
+ }
+ switch (data_type_dem) {
+ case CELL_TYPE:
+ d_dem = (double)((CELL *) inrast_dem)[col];
+ break;
+ case FCELL_TYPE:
+ d_dem = (double)((FCELL *) inrast_dem)[col];
+ break;
+ case DCELL_TYPE:
+ d_dem = ((DCELL *) inrast_dem)[col];
+ break;
+ }
+ if (Rast_is_d_null_value(&d_pw) ||
+ Rast_is_d_null_value(&d_pa) ||
+ Rast_is_d_null_value(&d_ta) ||
+ Rast_is_d_null_value(&d_dem)) {
+ Rast_set_c_null_value(&outrast1[col], 1);
+ }
+ else {
+ d_rh = rh(pw,pa,ta,dem);
+ d_esat = esat(ta);
+ d_eact = eact(d_esat,d_rh);
+ if(flag1->answer) d=d_esat;
+ else if(flag2->answer) d=d_eact;
+ else d=d_rh;
+ outrast1[col] = d;
+ }
+ }
+ Rast_put_row(outfd1, outrast1, data_type_output);
+ }
+ free(inrast_pw);
+ free(inrast_pa);
+ free(inrast_ta);
+ free(inrast_dem);
+ free(outrast1);
+ Rast_close(infd_pw);
+ Rast_close(infd_pa);
+ Rast_close(infd_ta);
+ Rast_close(infd_dem);
+ Rast_close(outfd1);
+ Rast_short_history(result1, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(result1, &history);
+ exit(EXIT_SUCCESS);
+}
+
+
Added: grass-addons/grass7/imagery/i.rh/rh.c
===================================================================
--- grass-addons/grass7/imagery/i.rh/rh.c (rev 0)
+++ grass-addons/grass7/imagery/i.rh/rh.c 2018-09-04 08:37:38 UTC (rev 73244)
@@ -0,0 +1,37 @@
+#include<stdio.h>
+#include<math.h>
+
+double rh(double PW,double Pa,double Ta,double dem){
+ /*
+ https{//www.researchgate.net/publication/227247013_High-resolution_Surface_Relative_Humidity_Computation_Using_MODIS_Image_in_Peninsular_Malaysia/
+ PW <- MOD05_L2 product
+ Pa <- MOD07 product
+ Pa <- 1013.3-0.1038*dem
+ Ta <- MOD07 product
+ Ta <- -0.0065*dem+TaMOD07 (if dem>400m)
+
+ https{//ladsweb.modaps.eosdis.nasa.gov/archive/allData/61/MOD05_L2/
+ https{//ladsweb.modaps.eosdis.nasa.gov/archive/allData/61/MOD07_L2/
+ */
+ /*Specific Humidity*/
+ q=0.001*(-0.0762*PW*PW+1.753*PW+12.405);
+ ta=-0.0065*dem+Ta;
+ a=17.2694*ta/(ta+238.3);
+ return(q*Pa/(380*exp(a)));
+}
+
+double esat(double tamean){
+ /*
+ esat{ saturated vapour pressure
+ tamean{ air temperature daily mean
+ */
+ return(610.78*exp(17.2694*tamean/(tamean+238.3));
+}
+double eact(double esat,double rh){
+ /*
+ eact{ actual vapour pressure
+ esat{ saturated vapour pressure
+ rh{ relative humidity
+ */
+ return(0.01*esat*rh);
+}
More information about the grass-commit
mailing list