Clarify what happens in gravfft with NaNs

doc/rst/source/explain_fft.rst_

@@ -29,9 +29,9 @@
     Control extension and tapering of data:
     Use modifiers to control how the extension and tapering are to be performed:
 
-        **+e** extends the grid by imposing edge-point symmetry [Default],
+        **+e** extends the grid by imposing edge-point symmetry [Default].
 
-        **+m** extends the grid by imposing edge mirror symmetry
+        **+m** extends the grid by imposing edge mirror symmetry.
 
         **+n** turns off data extension.
 

doc/rst/source/supplements/potential/gravfft.rst

@@ -17,7 +17,7 @@ Synopsis
 [ |-C|\ *n/wavelength/mean\_depth*/**t**\|\ **b**\|\ **w** ]
 [ |-D|\ *density*\|\ *rhogrid* ]
 [ |-E|\ *n_terms* ]
-[ |-F|\ [**f**\ [**+s**]\|\ **b**\|\ **g**\|\ **v**\|\ **n**\|\ **e**] ]
+[ |-F|\ [**f**\ [**+s**\|\ **z**]\|\ **b**\|\ **g**\|\ **v**\|\ **n**\|\ **e**] ]
 [ |-I|\ **w**\|\ **b**\|\ **c**\|\ **t**\|\ **k** ]
 [ |-N|\ *params* ]
 [ |-Q| ]
@@ -99,16 +99,18 @@ Optional Arguments
 
 **-E**\ *n_terms*
     Number of terms used in Parker expansion (limit is 10, otherwise
-    terms depending on n will blow out the program) [Default = 3]
+    terms depending on n will blow out the program) [Default = 3].
 
 .. _-F:
 
-**-F**\ [**f**\ [**+s**]\|\ **b**\|\ **g**\|\ **v**\|\ **n**\|\ **e**]
+**-F**\ [**f**\ [**+s**\|\ **z**]\|\ **b**\|\ **g**\|\ **v**\|\ **n**\|\ **e**]
     Specify desired geopotential field: compute geoid rather than gravity
 
        **f** = Free-air anomalies (mGal) [Default].  Append **+s** to add
        in the slab implied when removing the mean value from the topography.
-       This requires zero topography to mean no mass anomaly.
+       This requires zero topography to mean no mass anomaly. Alternatively,
+       to force the far-field to be exactly zero (i.e., the corner nodes of
+       the grid), select **+z** instead.
 
        **b** = Bouguer gravity anomalies (mGal).
 
@@ -146,7 +148,7 @@ Optional Arguments
     Writes out a grid with the flexural topography (with z positive up)
     whose average depth was set by **-Z**\ *zm* and model parameters by |-T|
     (and output by |-G|). That is the "gravimetric Moho". |-Q|
-    implicitly sets **-N+h**
+    implicitly sets **-N+h**.
 
 .. _-S:
 
@@ -167,7 +169,7 @@ Optional Arguments
     is > 1e10 it will be interpreted as the flexural rigidity (by default it is
     computed from *te* and Young modulus). Optionally, append *+m* to write a grid
     with the Moho's geopotential effect (see |-F|) from model selected by |-T|.
-    If *te* = 0 then the Airy response is returned. **-T+m** implicitly sets **-N+h**
+    If *te* = 0 then the Airy response is returned. **-T+m** implicitly sets **-N+h**.
 
 .. _-W:
 
@@ -207,6 +209,13 @@ meters, select |SYN_OPT-f|. If the data are close to either pole, you should
 consider projecting the grid file onto a rectangular coordinate system
 using :doc:`grdproject </grdproject>`.
 
+Handling of Grids with NaNs
+---------------------------
+
+Since we cannot take FFTs of 2-D grids that contain NaNs, we perform simple substitutions.
+If any of the input grids contain NaNs they will be replaced with zeros. In contrast, if **-D**
+passes a grid with density contrasts then we replace any NaNs with the minimum density in the grid.
+
 Data Detrending
 ---------------
 

src/potential/gravfft.c

@@ -75,10 +75,11 @@ struct GRAVFFT_CTRL {
 		bool active;
 		unsigned int n_terms;
 	} E;
-	struct GRAVFFT_F {	/* -F[f[+s]|b|g|e|n|v] */
+	struct GRAVFFT_F {	/* -F[f[+s|z]|b|g|e|n|v] */
 		bool active;
 		bool slab;
 		bool bouguer;
+		bool zero;
 		unsigned int mode;
 	} F;
 	struct GRAVFFT_G {	/* -G<outfile> */
@@ -273,6 +274,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVFFT_CTRL *Ctrl, struct GMT_OP
 					case 'f': default:
 						Ctrl->F.mode = GRAVFFT_FAA; 	   /* FAA */
 					 	if (opt->arg[1] == '+' && (opt->arg[2] == 's' || opt->arg[2] == '\0')) Ctrl->F.slab = true;
+					 	else if (opt->arg[1] == '+' && (opt->arg[2] == 'z' || opt->arg[2] == '\0')) Ctrl->F.zero = true;
 						break;
 				}
 				break;
@@ -420,7 +422,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
 	if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
 	GMT_Usage (API, 0, "usage: %s %s [<ingrid2>] -G%s [-C<n/wavelength/mean_depth/tbw>] "
-		"[-D<density>] [-E<n_terms>] [-F[f[+s]|b|g|v|n|e]] [-I<cbktw>] [-N%s] [-Q] [-S] "
+		"[-D<density>] [-E<n_terms>] [-F[f[+s|z]|b|g|v|n|e]] [-I<cbktw>] [-N%s] [-Q] [-S] "
 		"[-T<te/rl/rm/rw>[/<ri>][+m]] [%s] [-W<wd>[k]] [-Z<zm>[/<zl>]] [-fg] [%s]\n",
 		name, GMT_INGRID, GMT_OUTGRID, GMT_FFT_OPT, GMT_V_OPT, GMT_PAR_OPT);
 
@@ -444,11 +446,12 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 		"Append a co-registered density grid for a variable density contrast or a constant <density>.");
 	GMT_Usage (API, 1, "\n-E<n_terms>");
 	GMT_Usage (API, -2, "Number of terms used in Parker's expansion [Default = 3].");
-	GMT_Usage (API, 1, "\n-F[f[+s]|b|g|v|n|e]");
+	GMT_Usage (API, 1, "\n-F[f[+s|z]|b|g|v|n|e]");
 	GMT_Usage (API, -2, "Specify desired geopotential field:");
 	GMT_Usage (API, 3, "b: Bouguer anomalies (mGal).");
-	GMT_Usage (API, 3, "f: Free-air anomalies (mGal) [Default].");
+	GMT_Usage (API, 3, "f: Free-air anomalies (mGal) [Default]. Optionally choose one modifier:");
 	GMT_Usage (API, 4, "+s Adjust for implied slab correction [none].");
+	GMT_Usage (API, 4, "+z Adjust the field so the far-field is exactly zero [no shift].");
 	GMT_Usage (API, 3, "g: Geoid anomalies (m).");
 	GMT_Usage (API, 3, "v: Vertical Gravity Gradient (VGG; 1 Eovtos = 0.1 mGal/km).");
 	GMT_Usage (API, 3, "e: East deflections of the vertical (micro-radian).");
@@ -599,7 +602,12 @@ EXTERN_MSC int GMT_gravfft (void *V_API, int mode, void *args) {
 		if (!gmt_grd_domains_match (GMT, Orig[0], Rho, "surface and density")) {
 			Return (GMT_RUNTIME_ERROR);
 		}
-		for (m = 0; m < Rho->header->size; m++) if (gmt_M_is_fnan (Rho->data[m])) Rho->data[m] = Rho->header->z_min;	/* Replace any NaNs with the minimum density */
+		HH = gmt_get_H_hidden (Rho->header);
+		if (HH->has_NaNs == GMT_GRID_HAS_NANS) {
+			uint64_t n_nan = 0;
+			for (m = 0; m < Rho->header->size; m++) if (gmt_M_is_fnan (Rho->data[m])) Rho->data[m] = Rho->header->z_min, n_nan++;	/* Replace any NaNs with the minimum density */
+			GMT_Report (API, GMT_MSG_INFORMATION, "Density grid %s had %" PRIu64 " NaNs; these have been replaced by the minimum density %lg\n", Ctrl->In.file[k], n_nan, Rho->header->z_min);
+		}
 	}
 
 	/* Grids are compatible. Initialize FFT structs, grid headers, read data, and check for NaNs */
@@ -610,12 +618,13 @@ EXTERN_MSC int GMT_gravfft (void *V_API, int mode, void *args) {
                                       GMT_GRID_IS_COMPLEX_REAL, NULL, Ctrl->In.file[k], Orig[k])) == NULL)	/* Get data only */
 			Return (API->error);
 		HH = gmt_get_H_hidden (Orig[k]->header);
-		if (HH->has_NaNs == GMT_GRID_HAS_NANS) {
-			GMT_Report (API, GMT_MSG_ERROR, "Input grid %s has NaNs, cannot be used with FFTs\n", Ctrl->In.file[k]);
-			Return (GMT_RUNTIME_ERROR);
-		}
 		/* Note: If input grid(s) are read-only then we must duplicate them; otherwise Grid[k] points to Orig[k] */
 		(void) gmt_set_outgrid (GMT, Ctrl->In.file[k], false, 0, Orig[k], &Grid[k]);
+		if (HH->has_NaNs == GMT_GRID_HAS_NANS) {
+			uint64_t n_nan = 0;
+			for (m = 0; m < Grid[k]->header->size; m++) if (gmt_M_is_fnan (Grid[k]->data[m])) Grid[k]->data[m] = 0.0, n_nan++;	/* Replace any NaNs with zero */
+			GMT_Report (API, GMT_MSG_INFORMATION, "Input grid %s had %" PRIu64 " NaNs; these have been replaced with 0\n", Ctrl->In.file[k], n_nan);
+		}
 	}
 
 	/* From here we address the first grid via Grid[0] and the 2nd grid (if given) as Grid[1];
@@ -775,6 +784,12 @@ EXTERN_MSC int GMT_gravfft (void *V_API, int mode, void *args) {
 				else
 					for (m = 0; m < Grid[0]->header->size; m++) Grid[0]->data[m] += slab_gravity;
 			}
+			else if (Ctrl->F.zero) {	/* Shift so average FAA at the four corners is zero */
+				double far_field = 0.25 * (Grid[0]->data[gmt_M_ijp(Grid[0]->header, 0, 0)] + Grid[0]->data[gmt_M_ijp(Grid[0]->header, Grid[0]->header->n_rows-1, 0)]
+					+ Grid[0]->data[gmt_M_ijp(Grid[0]->header, 0, Grid[0]->header->n_columns-1)] + Grid[0]->data[gmt_M_ijp(Grid[0]->header, Grid[0]->header->n_rows-1, Grid[0]->header->n_columns-1)]);
+				GMT_Report (API, GMT_MSG_INFORMATION, "Subtract %g mGal from predicted FAA grid to force far-field to be zero\n", far_field);
+					for (m = 0; m < Grid[0]->header->size; m++) Grid[0]->data[m] -= far_field;
+			}
 			break;
 		case GRAVFFT_GEOID:
 			strcpy (Grid[0]->header->title, "Geoid anomalies");