[GRASS-SVN] r54966 - in grass-addons/grass7/imagery: . i.gravity i.rotate
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Feb 7 02:41:38 PST 2013
Author: ychemin
Date: 2013-02-07 02:41:38 -0800 (Thu, 07 Feb 2013)
New Revision: 54966
Added:
grass-addons/grass7/imagery/i.gravity/
grass-addons/grass7/imagery/i.gravity/Makefile
grass-addons/grass7/imagery/i.gravity/gravity.c
grass-addons/grass7/imagery/i.gravity/i.gravity.html
grass-addons/grass7/imagery/i.gravity/main.c
Modified:
grass-addons/grass7/imagery/i.rotate/main.c
Log:
Add i.gravity (Bouguer anomaly)
Added: grass-addons/grass7/imagery/i.gravity/Makefile
===================================================================
--- grass-addons/grass7/imagery/i.gravity/Makefile (rev 0)
+++ grass-addons/grass7/imagery/i.gravity/Makefile 2013-02-07 10:41:38 UTC (rev 54966)
@@ -0,0 +1,11 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.gravity
+
+LIBES = $(RASTERLIB) $(GISLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+
+default: cmd
Added: grass-addons/grass7/imagery/i.gravity/gravity.c
===================================================================
--- grass-addons/grass7/imagery/i.gravity/gravity.c (rev 0)
+++ grass-addons/grass7/imagery/i.gravity/gravity.c 2013-02-07 10:41:38 UTC (rev 54966)
@@ -0,0 +1,48 @@
+#include<stdio.h>
+#include<math.h>
+
+/*# International Gravity Formula (For Latitude correction)*/
+double g_lambda(double lambda){
+ /*International Gravity Formula (mGal)
+ lambda=latitude (dd)*/
+ return (978031.8*(1+0.0053024*sin(lambda)*sin(lambda)-0.0000059*sin(2*lambda)*sin(2*lambda)));
+}
+
+/*# Eotvos Correction*/
+double delta_g_eotvos(double alpha,double lambda,double v){
+ /* Eotvos correction (mGal)
+ v=velocity (kph)
+ lambda=latitude (dd)
+ alpha=direction of travel measured clockwise from North*/
+ return (4.040*v*sin(alpha)*cos(lambda)+0.001211*v*v);
+}
+
+/*# Free air Correction*/
+double free_air(h){
+ /* Free air Correction (mGal)
+ h=height (m)*/
+ return (0.3086*h);
+}
+
+/*# Bouguer Correction*/
+double delta_g_bouguer(double rho,double h){
+ /* Bouguer Correction (mGal)
+ G=gravity constant 6.67398 × 10^-11 [m3 kg-1 s-2]
+ rho=density of slab (Mg/m3)
+ h=height (m)*/
+ return (2*M_PI*6.67398*pow(10,-11)*rho*h);
+}
+
+double bouguer_anomaly(double g_obs,double freeair_corr,double bouguer_corr,double terrain_corr,double latitude_corr,double eotvos_corr){
+ /*
+ bouguer anomaly
+ g_obs=observed gravity
+ freeair_corr=free air correction
+ bouguer_corr=bouguer correction
+ terrain_corr=terrain correction
+ latitude_corr=latitude correction
+ eotvos_corr=eotvos correction
+ */
+ return (g_obs+freeair_corr-bouguer_corr+terrain_corr-latitude_corr+eotvos_corr);
+}
+
Added: grass-addons/grass7/imagery/i.gravity/i.gravity.html
===================================================================
--- grass-addons/grass7/imagery/i.gravity/i.gravity.html (rev 0)
+++ grass-addons/grass7/imagery/i.gravity/i.gravity.html 2013-02-07 10:41:38 UTC (rev 54966)
@@ -0,0 +1,41 @@
+<h2>DESCRIPTION</h2>
+
+<em>i.gravity</em> calculates the Bouguer gravity anomaly (after Mussett and Khan).
+<br>
+ Bouguer anomaly computation require:
+<ul>
+ <li> g_obs=observed gravity
+ <li> freeair_corr=free air correction
+ <li> bouguer_corr=bouguer correction
+ <li> terrain_corr=terrain correction
+ <li> latitude_corr=latitude correction
+ <li> eotvos_corr=eotvos correction
+</ul>
+
+<h2>NOTES</h2>
+
+<ul>
+ <li> International Gravity Formula (mGal), lambda=latitude (dd)
+ <li> Eotvos correction (mGal), v=velocity (kph), lambda=latitude (dd), alpha=direction of travel measured clockwise from North
+ <li> Free air Correction (mGal), h=height (m)
+ <li> Bouguer Correction (mGal), rho=density of slab (Mg/m3), h=height (m)
+ <li>
+</ul>
+
+<p>For more details on the algorithms see [1].
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="i.latlon.html">i.latlon</a><br>
+</em>
+
+<h2>REFERENCES</h2>
+
+ <p>[1] Mussett, A.E. and Khan, M.A, M.A. Looking into the Earth: An Introduction to Geological Geophysics.
+
+<h2>AUTHORS</h2>
+
+Yann Chemin, University of London at Birkbeck, United Kingdom.
+
+<p><i>Last changed: $Date: 2013-02-07 22:24:20 +0100 (Tue, 08 Nov 2011) $</i>
Added: grass-addons/grass7/imagery/i.gravity/main.c
===================================================================
--- grass-addons/grass7/imagery/i.gravity/main.c (rev 0)
+++ grass-addons/grass7/imagery/i.gravity/main.c 2013-02-07 10:41:38 UTC (rev 54966)
@@ -0,0 +1,222 @@
+/****************************************************************************
+ *
+ * MODULE: i.gravity
+ * AUTHOR(S): Yann Chemin - yann.chemin at gmail.com
+ * PURPOSE: Calculates the Bouguer anomaly of a gravity dataset
+ *
+ * COPYRIGHT: (C) 2002-2013 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>
+
+/* International Gravity Formula (For Latitude correction)*/
+double g_lambda(double lambda);
+/* Eotvos Correction*/
+double delta_g_eotvos(double alpha,double lambda,double v);
+/* Free air Correction*/
+double free_air(h);
+/* Bouguer Correction*/
+double delta_g_bouguer(double rho,double h);
+/* bouguer anomaly */
+double bouguer_anomaly(double g_obs,double freeair_corr,double bouguer_corr,double terrain_corr,double latitude_corr,double eotvos_corr);
+
+
+int main(int argc, char *argv[])
+{
+ int nrows, ncols;
+ int row, col;
+ struct GModule *module;
+ struct Option *in0, *in1, *in2, *in3, *in4, *in5, *in6, *output1;
+ struct History history; /*metadata */
+ char *result1; /*output raster name */
+ int infd_g_obs, infd_elevation, infd_latitude;
+ int infd_flight_azimuth, infd_flight_velocity;
+ int infd_rho, infd_terrain_corr;
+ int outfd1;
+ char *g_obs, *elevation, *latitude;
+ char *flight_azimuth, *flight_velocity, *rho, *terrain_corr;
+ void *inrast_g_obs, *inrast_elevation, *inrast_latitude;
+ void *inrast_flight_azimuth, *inrast_flight_velocity;
+ void *inrast_rho, *inrast_terrain_corr;
+
+ DCELL * outrast1;
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ G_add_keyword(_("imagery"));
+ G_add_keyword(_("gravity"));
+ G_add_keyword(_("Bouguer"));
+ G_add_keyword(_("Eotvos"));
+ G_add_keyword(_("Free Air"));
+ module->description =
+ _("Bouguer anomaly computation (full slab).");
+
+ /* Define the different options */
+ in0 = G_define_standard_option(G_OPT_R_INPUT);
+ in0->key = "gravity_obs";
+ in0->description = _("Name of the observed gravity map [mGal]");
+ in0->answer = "gravity_obs";
+
+ in1 = G_define_standard_option(G_OPT_R_INPUT);
+ in1->key = "elevation";
+ in1->description = _("Name of the elevation map [m]");
+ in1->answer = "dem";
+
+ in2 = G_define_standard_option(G_OPT_R_INPUT);
+ in2->key = "latitude";
+ in2->description = _("Name of the latitude map [dd.ddd]");
+ in2->answer = "lat";
+
+ in3 = G_define_standard_option(G_OPT_R_INPUT);
+ in3->key = "flight_azimuth";
+ in3->description = _("Name of the flight azimuth map (clockwise from North) [dd.ddd]");
+ in3->answer = "flight_azimuth";
+
+ in4 = G_define_standard_option(G_OPT_R_INPUT);
+ in4->key = "flight_velocity";
+ in4->description = _("Name of the flight velocity [kph]");
+ in4->answer = "flight_velocity";
+
+ in5 = G_define_standard_option(G_OPT_R_INPUT);
+ in5->key = "slab_density";
+ in5->description = _("Name of the slab density [kg/m3]");
+ in5->answer = "slab_density";
+
+ in6 = G_define_standard_option(G_OPT_R_INPUT);
+ in6->key = "terrain_corr";
+ in6->description = _("Name of the terrain correction map []");
+ in6->answer = "terrain_corr";
+
+ output1 = G_define_standard_option(G_OPT_R_OUTPUT);
+ output1->description =
+ _("Name of the Bouguer gravity anomaly [mGal]");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ g_obs = in0->answer;
+ elevation = in1->answer;
+ latitude = in2->answer;
+ flight_azimuth = in3->answer;
+ flight_velocity = in4->answer;
+ rho = in5->answer;
+ terrain_corr = in6->answer;
+ result1 = output1->answer;
+
+ infd_g_obs = Rast_open_old(g_obs, "");
+ inrast_g_obs = Rast_allocate_d_buf();
+
+ infd_elevation = Rast_open_old(elevation, "");
+ inrast_elevation = Rast_allocate_d_buf();
+
+ infd_latitude = Rast_open_old(latitude, "");
+ inrast_latitude = Rast_allocate_d_buf();
+
+ infd_flight_azimuth = Rast_open_old(flight_azimuth, "");
+ inrast_flight_azimuth = Rast_allocate_d_buf();
+
+ infd_flight_velocity = Rast_open_old(flight_velocity, "");
+ inrast_flight_velocity = Rast_allocate_d_buf();
+
+ infd_rho = Rast_open_old(rho, "");
+ inrast_rho = Rast_allocate_d_buf();
+
+ infd_terrain_corr = Rast_open_old(terrain_corr, "");
+ inrast_terrain_corr = Rast_allocate_d_buf();
+
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+ outrast1 = Rast_allocate_d_buf();
+
+ outfd1 = Rast_open_new(result1, DCELL_TYPE);
+
+ /* Process pixels */
+ for (row = 0; row < nrows; row++)
+ {
+ DCELL d;
+ DCELL d_g_obs;
+ DCELL d_elevation;
+ DCELL d_latitude;
+ DCELL d_flight_azimuth;
+ DCELL d_flight_velocity;
+ DCELL d_rho;
+ DCELL d_terrain_corr;
+ DCELL latitude_corr;
+ DCELL eotvos_corr;
+ DCELL freeair_corr;
+ DCELL bouguer_corr;
+ G_percent(row, nrows, 2);
+
+ /* read in maps */
+ Rast_get_d_row(infd_g_obs,inrast_g_obs,row);
+ Rast_get_d_row(infd_elevation,inrast_elevation,row);
+ Rast_get_d_row(infd_latitude,inrast_latitude,row);
+ Rast_get_d_row(infd_flight_azimuth,inrast_flight_azimuth,row);
+ Rast_get_d_row(infd_flight_velocity,inrast_flight_velocity,row);
+ Rast_get_d_row(infd_rho,inrast_rho,row);
+ Rast_get_d_row(infd_terrain_corr,inrast_terrain_corr,row);
+
+ /*process the data */
+ for (col = 0; col < ncols; col++)
+ {
+ d_g_obs = ((DCELL *) inrast_g_obs)[col];
+ d_elevation = ((DCELL *) inrast_elevation)[col];
+ d_latitude = ((DCELL *) inrast_latitude)[col];
+ d_flight_azimuth = ((DCELL *) inrast_flight_azimuth)[col];
+ d_flight_velocity = ((DCELL *) inrast_flight_velocity)[col];
+ d_rho = ((DCELL *) inrast_rho)[col];
+ d_terrain_corr = ((DCELL *) inrast_terrain_corr)[col];
+ if (Rast_is_d_null_value(&d_g_obs) ||
+ Rast_is_d_null_value(&d_elevation) ||
+ Rast_is_d_null_value(&d_latitude) ||
+ Rast_is_d_null_value(&d_flight_azimuth)||
+ Rast_is_d_null_value(&d_flight_velocity)||
+ Rast_is_d_null_value(&d_rho)||
+ Rast_is_d_null_value(&d_terrain_corr))
+ Rast_set_d_null_value(&outrast1[col], 1);
+ else {
+ /* International Gravity Formula (For Latitude correction)*/
+ latitude_corr=g_lambda(d_latitude);
+ /* Eotvos Correction*/
+ eotvos_corr=delta_g_eotvos(d_flight_azimuth,d_latitude,d_flight_velocity);
+ /* Free air Correction*/
+ freeair_corr=free_air(d_elevation);
+ /* Bouguer Correction*/
+ bouguer_corr=delta_g_bouguer(d_rho,d_elevation);
+ /* bouguer anomaly */
+ d=bouguer_anomaly(d_g_obs,freeair_corr,bouguer_corr,d_terrain_corr,latitude_corr,eotvos_corr);
+ outrast1[col] = d;
+ }
+ }
+ Rast_put_d_row(outfd1,outrast1);
+ }
+ G_free(inrast_elevation);
+ G_free(inrast_latitude);
+ G_free(inrast_flight_azimuth);
+ G_free(inrast_flight_velocity);
+ G_free(inrast_rho);
+ G_free(inrast_terrain_corr);
+ Rast_close(infd_elevation);
+ Rast_close(infd_latitude);
+ Rast_close(infd_flight_azimuth);
+ Rast_close(infd_flight_velocity);
+ Rast_close(infd_rho);
+ Rast_close(infd_terrain_corr);
+ G_free(outrast1);
+ Rast_close(outfd1);
+ Rast_short_history(result1, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(result1, &history);
+ exit(EXIT_SUCCESS);
+}
+
Modified: grass-addons/grass7/imagery/i.rotate/main.c
===================================================================
--- grass-addons/grass7/imagery/i.rotate/main.c 2013-02-07 09:56:17 UTC (rev 54965)
+++ grass-addons/grass7/imagery/i.rotate/main.c 2013-02-07 10:41:38 UTC (rev 54966)
@@ -20,6 +20,7 @@
#include <string.h>
#include <grass/gis.h>
#include <grass/raster.h>
+#include <grass/gprojects.h>
#include <grass/glocale.h>
#include <math.h>
@@ -89,6 +90,17 @@
}
}
DCELL d;
+
+ /* Get projection info for mapset */
+ if ((proj_info = G_get_projinfo()) == NULL)
+ G_fatal_error(_("Unable to get projection info of raster map"));
+
+ if ((unit_info = G_get_projunits()) == NULL)
+ G_fatal_error(_("Unable to get projection units of raster map"));
+
+ if (pj_get_kv(&proj, proj_info, unit_info) < 0)
+ G_fatal_error(_("Unable to get projection key values of raster map"));
+
/* Load input matrix with row&col shift to keep center of image*/
for (row = 0; row < nrows; row++)
{
More information about the grass-commit
mailing list