/*-------------------------------------------------------------------- * The MB-system: mbareaclean.c 2/27/2003 * $Id: mbareaclean.c 1829 2010-02-05 02:53:39Z caress $ * * Copyright (c) 2003-2009 by * David W. Caress (caress@mbari.org) * Monterey Bay Aquarium Research Institute * Moss Landing, CA 95039 * and Dale N. Chayes (dale@ldeo.columbia.edu) * Lamont-Doherty Earth Observatory * Palisades, NY 10964 * * See README file for copying and redistribution conditions. *--------------------------------------------------------------------*/ /* * mbareaclean identifies and flags artifacts in swath sonar bathymetry data. * The edit events are output to edit save files which can be applied * to the data by the program mbprocess. These are the same edit save * files created and/or modified by mbclean and mbedit. * The input data are one swath file or a datalist referencing multiple * swath files. An area is specified in longitude and latitude bounds, * along with a bin size in meters. The area is divided into a grid with * square cells of the specified bin size. As the data are read, each of * the soundings that fall within one of the bins is stored. Once all of * data are read, one or more statistical tests are performed on the soundings * within each bin, providing there are a sufficient number of soundings. * The user may specify one or both of the following actions: * 1) Previously unflagged soundings that fail a test are flagged as bad. * 2) Previously flagged soundings that pass all tests are unflagged. * If a sounding's flag status is changed, that flagging action is output * to the edit save file of the swath file containing that sounding. This * program will create edit save files if necessary, or append to those that * already exist. * * Author: D. W. Caress * Date: February 27, 2003 * Amsterdam Airport * * $Log: mbareaclean.c,v $ * Revision 5.13 2008/12/05 17:32:52 caress * Check-in mods 5 December 2008 including contributions from Gordon Keith. * * Revision 5.12 2008/07/10 18:16:33 caress * Proceeding towards 5.1.1beta20. * * Revision 5.10 2007/10/08 16:48:07 caress * State of the code on 8 October 2007. * * Revision 5.9 2006/08/09 22:41:27 caress * Fixed programs that read or write grids so that they do not use the GMT_begin() function; these programs will now work when GMT is built in the default fashion, when GMT is built in the default fashion, with "advisory file locking" enabled. * * Revision 5.8 2006/02/01 07:31:06 caress * Modifications suggested by Gordon Keith * * Revision 5.7 2006/01/27 19:13:04 caress * Version 5.0.8beta2 * * Revision 5.6 2006/01/18 15:17:00 caress * Added stdlib.h include. * * Revision 5.5 2005/11/05 01:07:54 caress * Programs changed to register arrays through mb_register_array() rather than allocating the memory directly with mb_realloc() or mb_malloc(). * * Revision 5.4 2004/12/02 06:39:28 caress * Fixes while supporting Reson 7k data. * * Revision 5.3 2004/09/16 00:57:46 caress * Fixed parsing of bounds argument. * * Revision 5.2 2003/07/26 18:01:22 caress * Changed beamflag handling code. * * Revision 5.1 2003/04/17 21:17:10 caress * Release 5.0.beta30 * * Revision 5.0 2003/03/10 20:47:08 caress * Initial version. * * */ /* standard include files */ #include #include #include #include #include #include #include /* mbio include files */ #include "../../include/mb_status.h" #include "../../include/mb_format.h" #include "../../include/mb_define.h" #include "../../include/mb_io.h" #include "../../include/mb_swap.h" #include "../../include/mb_process.h" #include "../../include/mb_info.h" /* allocation */ #define FILEALLOCNUM 16 #define PINGALLOCNUM 128 #define SNDGALLOCNUM 128 struct mbareaclean_file_struct { char filelist[MB_PATH_MAXLINE]; int file_format; int nping; int nping_alloc; int nnull; int nflag; int ngood; int nunflagged; int nflagged; double *ping_time_d; int *pingmultiplicity; double *ping_altitude; int nsndg; int nsndg_alloc; int sndg_countstart; int beams_bath; struct mbareaclean_sndg_struct *sndg; }; struct mbareaclean_sndg_struct { int sndg_file; int sndg_ping; int sndg_beam; double sndg_depth; double sndg_x; double sndg_y; char sndg_beamflag_org; char sndg_beamflag_esf; char sndg_beamflag; char sndg_edit; }; /* sounding atorage values and arrays */ int nfile = 0; int nfile_alloc = 0; struct mbareaclean_file_struct *files = NULL; int nsndg = 0; int nsndg_alloc = 0; int sndg_countstart = 0; int **gsndg = NULL; int *gsndgnum = NULL; int *gsndgnum_alloc = NULL; struct mbareaclean_sndg_struct *sndg = NULL; /* sounding pointer resolving function */ int getsoundingptr(int verbose, int soundingid, struct mbareaclean_sndg_struct **sndgptr, int *error); int flag_sounding(int verbose, int flag, int output_bad, int output_good, struct mbareaclean_sndg_struct *sndg, int *error); static char rcs_id[] = "$Id: mbareaclean.c 1829 2010-02-05 02:53:39Z caress $"; /*--------------------------------------------------------------------*/ int main (int argc, char **argv) { char program_name[] = "MBAREACLEAN"; char help_message[] = "MBAREACLEAN identifies and flags artifacts in swath bathymetry data"; char usage_message[] = "mbareaclean [-Fformat -Iinfile -Rwest/east/south/north -B -G -Sbinsize \n\t -Mthreshold/nmin -Dthreshold[/nmin[/nmax]] -Ttype -N[-]minbeam/maxbeam]"; extern char *optarg; int errflg = 0; int c; int help = 0; int flag = 0; /* MBIO status variables */ int status; int verbose = 0; int error = MB_ERROR_NO_ERROR; char *message = NULL; /* MBIO read control parameters */ void *mbio_ptr = NULL; void *store_ptr = NULL; int kind; int read_datalist = MB_NO; char read_file[MB_PATH_MAXLINE]; char swathfile[MB_PATH_MAXLINE]; char swathfileread[MB_PATH_MAXLINE]; void *datalist; int look_processed = MB_DATALIST_LOOK_UNSET; int read_data; double file_weight; int format; int formatread; int variable_beams; int traveltime; int beam_flagging; int pings; int lonflip; double bounds[4]; int btime_i[7]; int etime_i[7]; double btime_d; double etime_d; double speedmin; double timegap; struct mb_info_struct mb_info; int time_i[7]; double time_d; int pingsread; double navlon; double navlat; double speed; double heading; double distance; double altitude; double sonardepth; int beams_bath; int beams_amp; int pixels_ss; char *beamflag; char *beamflagorg; int *detect; double *bath; double *amp; double *bathlon; double *bathlat; double *ss; double *sslon; double *sslat; char comment[MB_COMMENT_MAXLINE]; /* mbareaclean control parameters */ int median_filter = MB_NO; double median_filter_threshold = 0.25; int median_filter_nmin = 10; int density_filter = MB_NO; int density_filter_nmax = 0; int plane_fit = MB_NO; double plane_fit_threshold = 0.05; int plane_fit_nmin = 10; int std_dev_filter = MB_NO; double std_dev_threshold = 2.0; int std_dev_nmin = 10; int output_good = MB_NO; int output_bad = MB_NO; int flag_detect = MB_DETECT_AMPLITUDE; int use_detect = MB_NO; int limit_beams = MB_NO; int beam_in = MB_YES; int min_beam = 0; int max_beam = 0; int max_beam_no = 0; double areabounds[4]; int areaboundsset = MB_NO; double binsize = 0.0; int binsizeset = MB_NO; double dx, dy; int nx, ny; double mtodeglon; double mtodeglat; double mean; double std_dev; int detect_status; int detect_error; /* median filter parameters */ int binnum; int binnummax; double *bindepths; double threshold; double median_depth; /* counting parameters */ int files_tot = 0; int pings_tot = 0; int beams_tot = 0; int beams_good_org_tot = 0; int beams_flag_org_tot = 0; int beams_null_org_tot = 0; int pings_file = 0; int beams_file = 0; int beams_good_org_file = 0; int beams_flag_org_file = 0; int beams_null_org_file = 0; /* save file control variables */ int esffile_open = MB_NO; char esffile[MB_PATH_MAXLINE]; struct mb_esf_struct esf; int action; double xx, yy; int flagsounding; double median_depth_low; double median_depth_high; int done; int ix, iy, ib, kgrid; double d1, d2; int i1, i2, n; int i, j; /* get current default values */ status = mb_defaults(verbose,&format,&pings,&lonflip,bounds, btime_i,etime_i,&speedmin,&timegap); /* reset all defaults but the format and lonflip */ strcpy(read_file,"datalist.mb-1"); format = 0; pings = 1; bounds[0] = -360.; bounds[1] = 360.; bounds[2] = -90.; bounds[3] = 90.; btime_i[0] = 1962; btime_i[1] = 2; btime_i[2] = 21; btime_i[3] = 10; btime_i[4] = 30; btime_i[5] = 0; btime_i[6] = 0; etime_i[0] = 2062; etime_i[1] = 2; etime_i[2] = 21; etime_i[3] = 10; etime_i[4] = 30; etime_i[5] = 0; etime_i[6] = 0; speedmin = 0.0; timegap = 1000000000.0; /* process argument list */ while ((c = getopt(argc, argv, "VvHhBbGgD:d:F:f:I:i:M:m:N:n:P:p:S:sT:t::R:r:")) != -1) switch (c) { case 'H': case 'h': help++; break; case 'V': case 'v': verbose++; break; case 'B': case 'b': output_bad = MB_YES; flag++; break; case 'D': case 'd': std_dev_filter = MB_YES; sscanf (optarg,"%lf/%d", &std_dev_threshold,&std_dev_nmin); flag++; break; case 'F': case 'f': sscanf (optarg,"%d", &format); flag++; break; case 'G': case 'g': output_good = MB_YES; flag++; break; case 'I': case 'i': sscanf (optarg,"%s", read_file); flag++; break; case 'M': case 'm': median_filter = MB_YES; n = sscanf (optarg,"%lf/%d/%d", &d1,&i1,&i2); if (n > 0) median_filter_threshold = d1; if (n > 1) median_filter_nmin = i1; if (n > 2) { density_filter = MB_YES; density_filter_nmax = i2; } flag++; break; case 'N': case 'n': limit_beams = MB_YES; sscanf (optarg,"%d/%d", &min_beam,&max_beam_no); if (optarg[0] == '-') { min_beam = -min_beam; beam_in = MB_NO; } if (max_beam_no < 0) max_beam_no = -max_beam_no; max_beam = max_beam_no; if (max_beam < min_beam) max_beam = min_beam; flag++; break; case 'P': case 'p': plane_fit = MB_YES; sscanf (optarg,"%lf", &plane_fit_threshold); n = sscanf (optarg,"%lf/%d/%lf", &d1,&i1,&d2); if (n > 0) plane_fit_threshold = d1; if (n > 1) plane_fit_nmin = i1; flag++; break; case 'R': case 'r': mb_get_bounds(optarg, areabounds); areaboundsset = MB_YES; flag++; break; case 'S': case 's': sscanf (optarg,"%lf", &binsize); binsizeset = MB_YES; flag++; break; case 'T': case 't': use_detect = MB_YES; sscanf (optarg,"%d", &flag_detect); flag++; break; case '?': errflg++; } /* if error flagged then print it and exit */ if (errflg) { fprintf(stderr,"usage: %s\n", usage_message); fprintf(stderr,"\nProgram <%s> Terminated\n", program_name); error = MB_ERROR_BAD_USAGE; exit(error); } /* turn on median filter if nothing specified */ if (median_filter == MB_NO && plane_fit == MB_NO && std_dev_filter == MB_NO) median_filter = MB_YES; /* turn on output bad if nothing specified */ if (output_bad == MB_NO && output_good == MB_NO) output_bad = MB_YES; /* print starting message */ if (verbose == 1 || help) { fprintf(stderr,"\nProgram %s\n",program_name); fprintf(stderr,"Version %s\n",rcs_id); fprintf(stderr,"MB-system Version %s\n",MB_VERSION); } /* print starting debug statements */ if (verbose >= 2) { fprintf(stderr,"\ndbg2 Program <%s>\n",program_name); fprintf(stderr,"dbg2 Version %s\n",rcs_id); fprintf(stderr,"dbg2 MB-system Version %s\n",MB_VERSION); fprintf(stderr,"dbg2 Control Parameters:\n"); fprintf(stderr,"dbg2 verbose: %d\n",verbose); fprintf(stderr,"dbg2 help: %d\n",help); fprintf(stderr,"dbg2 pings: %d\n",pings); fprintf(stderr,"dbg2 lonflip: %d\n",lonflip); fprintf(stderr,"dbg2 bounds[0]: %f\n",bounds[0]); fprintf(stderr,"dbg2 bounds[1]: %f\n",bounds[1]); fprintf(stderr,"dbg2 bounds[2]: %f\n",bounds[2]); fprintf(stderr,"dbg2 bounds[3]: %f\n",bounds[3]); fprintf(stderr,"dbg2 btime_i[0]: %d\n",btime_i[0]); fprintf(stderr,"dbg2 btime_i[1]: %d\n",btime_i[1]); fprintf(stderr,"dbg2 btime_i[2]: %d\n",btime_i[2]); fprintf(stderr,"dbg2 btime_i[3]: %d\n",btime_i[3]); fprintf(stderr,"dbg2 btime_i[4]: %d\n",btime_i[4]); fprintf(stderr,"dbg2 btime_i[5]: %d\n",btime_i[5]); fprintf(stderr,"dbg2 btime_i[6]: %d\n",btime_i[6]); fprintf(stderr,"dbg2 etime_i[0]: %d\n",etime_i[0]); fprintf(stderr,"dbg2 etime_i[1]: %d\n",etime_i[1]); fprintf(stderr,"dbg2 etime_i[2]: %d\n",etime_i[2]); fprintf(stderr,"dbg2 etime_i[3]: %d\n",etime_i[3]); fprintf(stderr,"dbg2 etime_i[4]: %d\n",etime_i[4]); fprintf(stderr,"dbg2 etime_i[5]: %d\n",etime_i[5]); fprintf(stderr,"dbg2 etime_i[6]: %d\n",etime_i[6]); fprintf(stderr,"dbg2 speedmin: %f\n",speedmin); fprintf(stderr,"dbg2 timegap: %f\n",timegap); fprintf(stderr,"dbg2 data format: %d\n",format); fprintf(stderr,"dbg2 input file: %s\n",read_file); fprintf(stderr,"dbg2 median_filter: %d\n",median_filter); fprintf(stderr,"dbg2 median_filter_threshold: %f\n",median_filter_threshold); fprintf(stderr,"dbg2 median_filter_nmin: %d\n",median_filter_nmin); fprintf(stderr,"dbg2 density_filter: %d\n",density_filter); fprintf(stderr,"dbg2 density_filter_nmax: %d\n",density_filter_nmax); fprintf(stderr,"dbg2 plane_fit: %d\n",plane_fit); fprintf(stderr,"dbg2 plane_fit_threshold: %f\n",plane_fit_threshold); fprintf(stderr,"dbg2 plane_fit_nmin: %d\n",plane_fit_nmin); fprintf(stderr,"dbg2 std_dev_filter: %d\n",std_dev_filter); fprintf(stderr,"dbg2 std_dev_threshold: %f\n",std_dev_threshold); fprintf(stderr,"dbg2 std_dev_nmin: %d\n",std_dev_nmin); fprintf(stderr,"dbg2 use_detect: %d\n",use_detect); fprintf(stderr,"dbg2 flag_detect: %d\n",flag_detect); fprintf(stderr,"dbg2 limit_beams: %d\n",limit_beams); fprintf(stderr,"dbg2 beam_in: %d\n",beam_in); fprintf(stderr,"dbg2 min_beam: %d\n",min_beam); fprintf(stderr,"dbg2 max_beam_no %d\n",max_beam_no); fprintf(stderr,"dbg2 output_good: %d\n",output_good); fprintf(stderr,"dbg2 output_bad: %d\n",output_bad); fprintf(stderr,"dbg2 areaboundsset: %d\n",areaboundsset); fprintf(stderr,"dbg2 areabounds[0]: %f\n",areabounds[0]); fprintf(stderr,"dbg2 areabounds[1]: %f\n",areabounds[1]); fprintf(stderr,"dbg2 areabounds[2]: %f\n",areabounds[2]); fprintf(stderr,"dbg2 areabounds[3]: %f\n",areabounds[3]); fprintf(stderr,"dbg2 binsizeset: %d\n",binsizeset); fprintf(stderr,"dbg2 binsize: %f\n",binsize); } /* if help desired then print it and exit */ if (help) { fprintf(stderr,"\n%s\n",help_message); fprintf(stderr,"\nusage: %s\n", usage_message); exit(error); } /* if bounds not set get bounds of input data */ if (areaboundsset == MB_NO) { formatread = format; status = mb_get_info_datalist(verbose, read_file, &formatread, &mb_info, lonflip, &error); areabounds[0] = mb_info.lon_min; areabounds[1] = mb_info.lon_max; areabounds[2] = mb_info.lat_min; areabounds[3] = mb_info.lat_max; if (binsizeset == MB_NO) binsize = 0.2 * mb_info.altitude_max; } /* calculate grid properties */ mb_coor_scale(verbose,0.5*(areabounds[2]+areabounds[3]),&mtodeglon,&mtodeglat); if (binsize <= 0.0) binsize = (areabounds[1] - areabounds[0]) / 101 / mtodeglon; dx = binsize * mtodeglon; dy = binsize * mtodeglat; nx = 1 + (int)((areabounds[1] - areabounds[0]) / dx); ny = 1 + (int)((areabounds[3] - areabounds[2]) / dy); if (nx > 1 && ny > 1) { dx = (areabounds[1] - areabounds[0]) / (nx - 1); dy = (areabounds[3] - areabounds[2]) / (ny - 1); } /* allocate grid arrays */ nsndg = 0; nsndg_alloc = 0; status = mb_mallocd(verbose,__FILE__,__LINE__, nx * ny * sizeof(int *), (void **)&gsndg, &error); if (status == MB_SUCCESS) status = mb_mallocd(verbose,__FILE__,__LINE__, nx * ny * sizeof(int), (void **)&gsndgnum, &error); if (status == MB_SUCCESS) status = mb_mallocd(verbose,__FILE__,__LINE__, nx * ny * sizeof(int), (void **)&gsndgnum_alloc, &error); /* if error initializing memory then quit */ if (error != MB_ERROR_NO_ERROR) { mb_error(verbose,error,&message); fprintf(stderr,"\nMBIO Error allocating data arrays:\n%s\n",message); fprintf(stderr,"\nProgram <%s> Terminated\n", program_name); exit(error); } /* if error initializing memory then quit */ for (i=0;i= 0) { fprintf(stderr,"Area of interest:\n"); fprintf(stderr," Minimum Longitude: %.6f Maximum Longitude: %.6f\n", areabounds[0],areabounds[1]); fprintf(stderr," Minimum Latitude: %.6f Maximum Latitude: %.6f\n", areabounds[2],areabounds[3]); fprintf(stderr," Bin Size: %f\n", binsize); fprintf(stderr," Dimensions: %d %d\n", nx, ny); fprintf(stderr,"Cleaning algorithms:\n"); if (median_filter == MB_YES) { fprintf(stderr," Median filter: ON\n"); fprintf(stderr," Median filter threshold: %f\n", median_filter_threshold); fprintf(stderr," Median filter minimum N: %d\n", median_filter_nmin); } else fprintf(stderr," Median filter: OFF\n"); if (density_filter == MB_YES) { fprintf(stderr," Density filter: ON\n"); fprintf(stderr," Density filter maximum N: %d\n", density_filter_nmax); } else fprintf(stderr," Density filter: OFF\n"); if (plane_fit == MB_YES) { fprintf(stderr," Plane fit: ON\n"); fprintf(stderr," Plane fit threshold: %f\n", median_filter_threshold); fprintf(stderr," Plane fit minimum N: %d\n", median_filter_nmin); } else fprintf(stderr," Plane fit: OFF\n"); if (std_dev_filter == MB_YES) { fprintf(stderr," Standard deviation filter: ON\n"); fprintf(stderr," Standard deviation filter threshold: %f\n", std_dev_threshold); fprintf(stderr," Standard deviation filter minimum N: %d\n", std_dev_nmin); } else fprintf(stderr," Standard deviation filter: OFF\n"); fprintf(stderr,"Restrictions:\n"); if (use_detect) { fprintf(stderr," Only flag if bottom detection algorithn is: "); if (flag_detect == MB_DETECT_UNKNOWN) fprintf(stderr,"UNKNOWN\n"); else if (flag_detect == MB_DETECT_AMPLITUDE) fprintf(stderr,"AMPLITUDE\n"); else if (flag_detect == MB_DETECT_PHASE) fprintf(stderr,"PHASE\n"); else fprintf(stderr,"%d\n", flag_detect); } if (limit_beams == MB_YES) { fprintf(stderr," Only flag if beams "); if (beam_in == MB_YES) fprintf(stderr,"between"); else fprintf(stderr,"outside"); fprintf(stderr," beams %d - %d\n", min_beam, max_beam_no); } else fprintf(stderr," Flag all beams\n"); fprintf(stderr,"Output:\n"); if (output_bad == MB_YES) fprintf(stderr," Flag unflagged soundings identified as bad: ON\n"); else fprintf(stderr," Flag unflagged soundings identified as bad: OFF\n"); if (output_good == MB_YES) fprintf(stderr," Unflag flagged soundings identified as good: ON\n"); else fprintf(stderr," Unflag flagged soundings identified as good: OFF\n"); } /* get format if required */ if (format == 0) mb_get_format(verbose,read_file,NULL,&format,&error); /* determine whether to read one file or a list of files */ if (format < 0) read_datalist = MB_YES; /* open file list */ if (read_datalist == MB_YES) { if ((status = mb_datalist_open(verbose,&datalist, read_file,look_processed,&error)) != MB_SUCCESS) { error = MB_ERROR_OPEN_FAIL; fprintf(stderr,"\nUnable to open data list file: %s\n", read_file); fprintf(stderr,"\nProgram <%s> Terminated\n", program_name); exit(error); } if ((status = mb_datalist_read(verbose,datalist, swathfile,&format,&file_weight,&error)) == MB_SUCCESS) read_data = MB_YES; else read_data = MB_NO; } /* else copy single filename to be read */ else { strcpy(swathfile, read_file); read_data = MB_YES; } /* loop over all files to be read */ while (read_data == MB_YES) { /* check format and get format flags */ if ((status = mb_format_flags(verbose,&format, &variable_beams, &traveltime, &beam_flagging, &error)) != MB_SUCCESS) { mb_error(verbose,error,&message); fprintf(stderr,"\nMBIO Error returned from function regarding input format %d:\n%s\n",format,message); fprintf(stderr,"\nProgram <%s> Terminated\n", program_name); exit(error); } /* check for "fast bathymetry" or "fbt" file */ strcpy(swathfileread, swathfile); formatread = format; if (use_detect == MB_NO) mb_get_fbt(verbose, swathfileread, &formatread, &error); /* initialize reading the input swath sonar file */ if ((status = mb_read_init( verbose,swathfileread,formatread,pings,lonflip,bounds, btime_i,etime_i,speedmin,timegap, &mbio_ptr,&btime_d,&etime_d, &beams_bath,&beams_amp,&pixels_ss,&error)) != MB_SUCCESS) { mb_error(verbose,error,&message); fprintf(stderr,"\nMBIO Error returned from function :\n%s\n",message); fprintf(stderr,"\nMultibeam File <%s> not initialized for reading\n",swathfileread); fprintf(stderr,"\nProgram <%s> Terminated\n", program_name); exit(error); } /* initialize and increment counting variables */ pings_file = 0; beams_file = 0; /* give the statistics */ if (verbose >= 0) { fprintf(stderr,"\nProcessing %s\n",swathfileread); } /* allocate memory for data arrays */ if (error == MB_ERROR_NO_ERROR) status = mb_register_array(verbose, mbio_ptr, MB_MEM_TYPE_BATHYMETRY, sizeof(char), (void **)&beamflag, &error); if (error == MB_ERROR_NO_ERROR) status = mb_register_array(verbose, mbio_ptr, MB_MEM_TYPE_BATHYMETRY, sizeof(char), (void **)&detect, &error); if (error == MB_ERROR_NO_ERROR) status = mb_register_array(verbose, mbio_ptr, MB_MEM_TYPE_BATHYMETRY, sizeof(double), (void **)&bath, &error); if (error == MB_ERROR_NO_ERROR) status = mb_register_array(verbose, mbio_ptr, MB_MEM_TYPE_AMPLITUDE, sizeof(double), (void **)&, &error); if (error == MB_ERROR_NO_ERROR) status = mb_register_array(verbose, mbio_ptr, MB_MEM_TYPE_BATHYMETRY, sizeof(double), (void **)&bathlon, &error); if (error == MB_ERROR_NO_ERROR) status = mb_register_array(verbose, mbio_ptr, MB_MEM_TYPE_BATHYMETRY, sizeof(double), (void **)&bathlat, &error); if (error == MB_ERROR_NO_ERROR) status = mb_register_array(verbose, mbio_ptr, MB_MEM_TYPE_SIDESCAN, sizeof(double), (void **)&ss, &error); if (error == MB_ERROR_NO_ERROR) status = mb_register_array(verbose, mbio_ptr, MB_MEM_TYPE_SIDESCAN, sizeof(double), (void **)&sslon, &error); if (error == MB_ERROR_NO_ERROR) status = mb_register_array(verbose, mbio_ptr, MB_MEM_TYPE_SIDESCAN, sizeof(double), (void **)&sslat, &error); if (error == MB_ERROR_NO_ERROR) status = mb_register_array(verbose, mbio_ptr, MB_MEM_TYPE_BATHYMETRY, sizeof(char), (void **)&beamflagorg, &error); /* if error initializing memory then quit */ if (error != MB_ERROR_NO_ERROR) { mb_error(verbose,error,&message); fprintf(stderr,"\nMBIO Error allocating data arrays:\n%s\n",message); fprintf(stderr,"\nProgram <%s> Terminated\n", program_name); exit(error); } /* update memory for files */ if (nfile >= nfile_alloc) { nfile_alloc += FILEALLOCNUM; status = mb_reallocd(verbose, __FILE__, __LINE__, nfile_alloc * sizeof(struct mbareaclean_file_struct), (void **)&files, &error); /* if error initializing memory then quit */ if (error != MB_ERROR_NO_ERROR) { mb_error(verbose,error,&message); fprintf(stderr,"\nMBIO Error allocating data arrays:\n%s\n",message); fprintf(stderr,"\nProgram <%s> Terminated\n", program_name); exit(error); } } /* initialize current file */ strcpy(files[nfile].filelist, swathfile); files[nfile].file_format = format; files[nfile].nping = 0; files[nfile].nping_alloc = PINGALLOCNUM; files[nfile].nnull = 0; files[nfile].nflag = 0; files[nfile].ngood = 0; files[nfile].nflagged = 0; files[nfile].nunflagged = 0; files[nfile].ping_time_d = NULL; files[nfile].pingmultiplicity = NULL; files[nfile].ping_altitude = NULL; files[nfile].nsndg = 0; files[nfile].nsndg_alloc = SNDGALLOCNUM; files[nfile].sndg_countstart = nsndg; files[nfile].beams_bath = beams_bath; files[nfile].sndg = NULL; status = mb_mallocd(verbose,__FILE__,__LINE__, files[nfile].nping_alloc * sizeof(double), (void **)&(files[nfile].ping_time_d), &error); if (status == MB_SUCCESS) status = mb_mallocd(verbose,__FILE__,__LINE__, files[nfile].nping_alloc * sizeof(int), (void **)&(files[nfile].pingmultiplicity), &error); if (status == MB_SUCCESS) status = mb_mallocd(verbose,__FILE__,__LINE__, files[nfile].nping_alloc * sizeof(double), (void **)&(files[nfile].ping_altitude), &error); if (status == MB_SUCCESS) status = mb_mallocd(verbose,__FILE__,__LINE__, files[nfile].nsndg_alloc * sizeof(struct mbareaclean_sndg_struct), (void **)&(files[nfile].sndg), &error); if (error != MB_ERROR_NO_ERROR) { mb_error(verbose,error,&message); fprintf(stderr,"\nMBIO Error allocating data arrays:\n%s\n",message); fprintf(stderr,"\nProgram <%s> Terminated\n", program_name); exit(error); } nfile++; /* now deal with old edit save file */ if (status == MB_SUCCESS) { /* handle esf edits */ status = mb_esf_load(verbose, swathfile, MB_YES, MB_NO, esffile, &esf, &error); } /* read */ done = MB_NO; files_tot++; pings_file = 0; beams_file = 0; beams_good_org_file = 0; beams_flag_org_file = 0; beams_null_org_file = 0; while (done == MB_NO) { if (verbose > 1) fprintf(stderr,"\n"); /* read next record */ error = MB_ERROR_NO_ERROR; status = mb_read(verbose,mbio_ptr,&kind, &pingsread,time_i,&time_d, &navlon,&navlat, &speed,&heading, &distance,&altitude,&sonardepth, &beams_bath,&beams_amp,&pixels_ss, beamflag,bath,amp,bathlon,bathlat, ss,sslon,sslat, comment,&error); if (verbose >= 2) { fprintf(stderr,"\ndbg2 current data status:\n"); fprintf(stderr,"dbg2 kind: %d\n",kind); fprintf(stderr,"dbg2 status: %d\n",status); } if (status == MB_SUCCESS && kind == MB_DATA_DATA) { for (i=0;i= files[nfile-1].nping_alloc) { files[nfile-1].nping_alloc += PINGALLOCNUM; status = mb_reallocd(verbose, __FILE__, __LINE__, files[nfile-1].nping_alloc * sizeof(double), (void **)&(files[nfile-1].ping_time_d), &error); if (status == MB_SUCCESS) status = mb_reallocd(verbose, __FILE__, __LINE__, files[nfile-1].nping_alloc * sizeof(int), (void **)&(files[nfile-1].pingmultiplicity), &error); if (status == MB_SUCCESS) status = mb_reallocd(verbose, __FILE__, __LINE__, files[nfile-1].nping_alloc * sizeof(double), (void **)&(files[nfile-1].ping_altitude), &error); if (error != MB_ERROR_NO_ERROR) { mb_error(verbose,error,&message); fprintf(stderr,"\nMBIO Error allocating data arrays:\n%s\n",message); fprintf(stderr,"\nProgram <%s> Terminated\n", program_name); exit(error); } } /* store the ping data */ files[nfile-1].ping_time_d[files[nfile-1].nping] = time_d; if (files[nfile-1].nping > 0 && files[nfile-1].ping_time_d[files[nfile-1].nping] == files[nfile-1].ping_time_d[files[nfile-1].nping - 1]) { files[nfile-1].pingmultiplicity[files[nfile-1].nping] = files[nfile-1].pingmultiplicity[files[nfile-1].nping - 1] + 1; } else { files[nfile-1].pingmultiplicity[files[nfile-1].nping] = 0; } files[nfile-1].ping_altitude[files[nfile-1].nping] = altitude; files[nfile-1].nping++; /* check beam range */ if (limit_beams == MB_YES && max_beam_no == 0) max_beam = beams_bath - min_beam; /* now loop over the beams and store the soundings in the grid bins */ for (ib=0;ib= 0 && ix < nx && iy >= 0 && iy < ny) { if (files[nfile-1].nsndg >= files[nfile-1].nsndg_alloc) { files[nfile-1].nsndg_alloc += SNDGALLOCNUM; status = mb_reallocd(verbose, __FILE__, __LINE__, files[nfile-1].nsndg_alloc * sizeof(struct mbareaclean_sndg_struct), (void **)&files[nfile-1].sndg, &error); if (error != MB_ERROR_NO_ERROR) { mb_error(verbose,error,&message); fprintf(stderr,"\nMBIO Error allocating sounding arrays:\n%s\n",message); fprintf(stderr,"\nProgram <%s> Terminated\n", program_name); exit(error); } } /* allocate space for sounding if needed */ if (gsndgnum[kgrid] >= gsndgnum_alloc[kgrid]) { gsndgnum_alloc[kgrid] += SNDGALLOCNUM; status = mb_reallocd(verbose, __FILE__, __LINE__, gsndgnum_alloc[kgrid] * sizeof(int), (void **)&gsndg[kgrid], &error); if (error != MB_ERROR_NO_ERROR) { mb_error(verbose,error,&message); fprintf(stderr,"\nMBIO Error allocating sounding arrays:\n%s\n",message); fprintf(stderr,"\nProgram <%s> Terminated\n", program_name); exit(error); } } /* store sounding data */ sndg = &(files[nfile-1].sndg[files[nfile-1].nsndg]); sndg->sndg_file = nfile - 1; sndg->sndg_ping = files[nfile - 1].nping - 1; sndg->sndg_beam = ib; sndg->sndg_depth = bath[ib]; sndg->sndg_x = bathlon[ib]; sndg->sndg_y = bathlat[ib]; sndg->sndg_beamflag_org = beamflag[ib]; sndg->sndg_beamflag_esf = beamflagorg[ib]; sndg->sndg_beamflag = beamflagorg[ib]; sndg->sndg_edit = MB_YES; if (use_detect && detect[ib] != flag_detect) sndg->sndg_edit = MB_NO; if (limit_beams == MB_YES) { if (min_beam <= ib && ib <= max_beam) { if (beam_in == MB_NO) sndg->sndg_edit = MB_NO; } else { if (beam_in == MB_YES) sndg->sndg_edit = MB_NO; } } /* fprintf(stderr,"beam: %d edit: %d\n", ib, sndg->sndg_edit);*/ files[nfile-1].nsndg++; nsndg++; gsndg[kgrid][gsndgnum[kgrid]] = files[nfile-1].sndg_countstart + files[nfile-1].nsndg - 1; gsndgnum[kgrid]++; /*fprintf(stderr,"NEW sounding:%d:%d file:%d time_d:%f depth:%f\n", nsndg-1,gsndg[kgrid][gsndgnum[kgrid]-1],sndg->sndg_file, files[sndg->sndg_file].ping_time_d[sndg->sndg_ping], sndg->sndg_depth);*/ } } } } else if (error > MB_ERROR_NO_ERROR) { done = MB_YES; } /* process a record */ /* reset counters and data */ } /* close the files */ status = mb_close(verbose,&mbio_ptr,&error); mb_esf_close(verbose, &esf, &error); /* check memory */ if (verbose >= 4) status = mb_memory_list(verbose,&error); /* give the statistics */ if (verbose >= 0) { fprintf(stderr,"pings:%4d beams: %7d good %7d flagged %7d null \n", pings_file,beams_good_org_file ,beams_flag_org_file,beams_null_org_file); } /* figure out whether and what to read next */ if (read_datalist == MB_YES) { if ((status = mb_datalist_read(verbose,datalist, swathfile,&format,&file_weight,&error)) == MB_SUCCESS) read_data = MB_YES; else read_data = MB_NO; } else { read_data = MB_NO; } /* end loop over files in list */ } if (read_datalist == MB_YES) mb_datalist_close(verbose,&datalist,&error); /* loop over grid cells to find maximum number of soundings */ binnummax = 0; for (ix=0;ix Terminated\n", program_name); exit(error); } /* deal with median filter */ if (median_filter == MB_YES) { /* loop over grid cells applying median filter test */ for (ix=0;ixsndg_beamflag)) { bindepths[binnum] = sndg->sndg_depth; binnum++; } /*fprintf(stderr,"ix:%d iy:%d kgrid:%d soundingid:%d beamflag:%d binnum:%d\n", ix,iy,kgrid,gsndg[kgrid][i],sndg->sndg_beamflag,binnum);*/ } /* apply median filter only if there are enough soundings */ if (binnum >= median_filter_nmin) { /* run qsort */ qsort((char *)bindepths,binnum,sizeof(double), (void *)mb_double_compare); median_depth = bindepths[binnum / 2]; if (density_filter == MB_YES && binnum / 2 - density_filter_nmax / 2 >= 0) median_depth_low = bindepths[binnum / 2 + density_filter_nmax / 2]; else median_depth_low = bindepths[0]; if (density_filter == MB_YES && binnum / 2 + density_filter_nmax / 2 < binnum) median_depth_high = bindepths[binnum / 2 + density_filter_nmax / 2]; else median_depth_high = bindepths[binnum-1]; /* if (binnum>0) fprintf(stderr,"bin: %d %d %d pos: %f %f nsoundings:%d median:%f\n", ix,iy,kgrid,xx,yy,binnum,median_depth);*/ /* process the soundings */ for (i=0;isndg_file].ping_altitude[sndg->sndg_ping]); flagsounding = MB_NO; if (fabs(sndg->sndg_depth - median_depth) > threshold) flagsounding = MB_YES; if (density_filter == MB_YES && (sndg->sndg_depth > median_depth_high || sndg->sndg_depth < median_depth_low)) flagsounding = MB_YES; /*fprintf(stderr,"sounding:%d file:%d time_d:%f depth:%f median:%f altitude:%f threshold:%f", gsndg[kgrid][i],sndg->sndg_file, files[sndg->sndg_file].ping_time_d[sndg->sndg_ping], sndg->sndg_depth,median_depth, files[sndg->sndg_file].ping_altitude[sndg->sndg_ping], threshold);*/ flag_sounding(verbose, flagsounding, output_bad, output_good, sndg, &error); } } } } /* deal with standard deviation filter */ if (std_dev_filter == MB_YES) { /* loop over grid cells applying std dev filter test */ for (ix=0;ixsndg_beamflag)) { mean += sndg->sndg_depth; binnum++; } } mean /= binnum; /* get standard deviation */ std_dev = 0.0; for (i=0;isndg_beamflag)) std_dev += (sndg->sndg_depth - mean) * (sndg->sndg_depth - mean); } std_dev = sqrt(std_dev / binnum); threshold = std_dev * std_dev_threshold; if (binnum>0) fprintf(stderr,"bin: %d %d %d pos: %f %f nsoundings:%d / %d mean:%f std_dev:%f\n", ix,iy,kgrid,xx,yy,binnum,gsndgnum[kgrid],mean,std_dev); /* apply standard deviation threshold only if there are enough soundings */ if (binnum >= std_dev_nmin) { /* process the soundings */ for (i=0;isndg_depth - mean) > threshold, output_bad, output_good, sndg, &error); } } } } /* loop over files checking for changed soundings */ for (i=0; i < nfile; i++) { /* open esf file */ status = mb_esf_load(verbose, files[i].filelist, MB_NO, MB_YES, esffile, &esf, &error); if (status == MB_SUCCESS && esf.esffp != NULL) esffile_open = MB_YES; if (status == MB_FAILURE && error == MB_ERROR_OPEN_FAIL) { esffile_open = MB_NO; fprintf(stderr, "\nUnable to open new edit save file %s\n", esf.esffile); } /* loop over all of the soundings */ for (j=0;jsndg_beamflag != sndg->sndg_beamflag_org) { if (mb_beam_ok(sndg->sndg_beamflag)) { action = MBP_EDIT_UNFLAG; } else if (mb_beam_check_flag_manual(sndg->sndg_beamflag)) { action = MBP_EDIT_FLAG; } else if (mb_beam_check_flag_filter(sndg->sndg_beamflag)) { action = MBP_EDIT_FILTER; } mb_esf_save(verbose, &esf, files[i].ping_time_d[sndg->sndg_ping], sndg->sndg_beam + files[i].pingmultiplicity[sndg->sndg_ping] * MB_ESF_MULTIPLICITY_FACTOR, action, &error); } } /* close esf file */ mb_esf_close(verbose, &esf, &error); /* update mbprocess parameter file */ if (esffile_open == MB_YES) { /* update mbprocess parameter file */ status = mb_pr_update_format(verbose, files[i].filelist, MB_YES, files[i].file_format, &error); status = mb_pr_update_edit(verbose, files[i].filelist, MBP_EDIT_ON, esffile, &error); } } /* give the total statistics */ if (verbose >= 0) { fprintf(stderr,"\nMBareaclean Processing Totals:\n"); fprintf(stderr,"-------------------------\n"); fprintf(stderr,"%d total swath data files processed\n",files_tot); fprintf(stderr,"%d total pings processed\n",pings_tot); fprintf(stderr,"%d total soundings processed\n",beams_tot); fprintf(stderr,"-------------------------\n"); for (i=0;i= 4) status = mb_memory_list(verbose,&error); /* print output debug statements */ if (verbose >= 2) { fprintf(stderr,"\ndbg2 Program <%s> completed\n", program_name); fprintf(stderr,"dbg2 Ending status:\n"); fprintf(stderr,"dbg2 status: %d\n",status); } /* end it all */ exit(error); } /*--------------------------------------------------------------------*/ int getsoundingptr(int verbose, int soundingid, struct mbareaclean_sndg_struct **sndgptr, int *error) { /* local variables */ char *function_name = "getsoundingptr"; int status = MB_SUCCESS; int i, j; /* print input debug statements */ if (verbose >= 2) { fprintf(stderr,"\ndbg2 MBIO function <%s> called\n", function_name); fprintf(stderr,"dbg2 Input arguments:\n"); fprintf(stderr,"dbg2 verbose: %d\n",verbose); fprintf(stderr,"dbg2 soundingid: %d\n",soundingid); fprintf(stderr,"dbg2 sndgptr: %lu\n",(size_t)sndgptr); } /* loop over the files until the sounding is found */ *sndgptr = NULL; for (i=0; i < nfile && *sndgptr == NULL; i++) { if (soundingid >= files[i].sndg_countstart && soundingid < files[i].sndg_countstart + files[i].nsndg) { j = soundingid - files[i].sndg_countstart; *sndgptr = &(files[i].sndg[j]); } } /* print output debug statements */ if (verbose >= 2) { fprintf(stderr,"\ndbg2 MBIO function <%s> completed\n", function_name); fprintf(stderr,"dbg2 Return values:\n"); fprintf(stderr,"dbg2 *sndgptr: %lu\n",(size_t)sndgptr); fprintf(stderr,"dbg2 error: %d\n",*error); fprintf(stderr,"dbg2 Return status:\n"); fprintf(stderr,"dbg2 status: %d\n",status); } /* return */ return(status); } /*--------------------------------------------------------------------*/ int flag_sounding(int verbose, int flag, int output_bad, int output_good, struct mbareaclean_sndg_struct *sndg, int *error) { /* local variables */ char *function_name = "flag_sounding"; int status = MB_SUCCESS; /* print input debug statements */ if (verbose >= 2) { fprintf(stderr,"\ndbg2 MBIO function <%s> called\n", function_name); fprintf(stderr,"dbg2 Input arguments:\n"); fprintf(stderr,"dbg2 verbose: %d\n",verbose); fprintf(stderr,"dbg2 flag: %d\n",flag); fprintf(stderr,"dbg2 output_bad: %d\n",output_bad); fprintf(stderr,"dbg2 output_good: %d\n",output_good); fprintf(stderr,"dbg2 sndg->sndg_edit: %d\n",sndg->sndg_edit); fprintf(stderr,"dbg2 sndg->sndg_beam: %d\n",sndg->sndg_beam); fprintf(stderr,"dbg2 sndg->sndg_beamflag: %d\n",sndg->sndg_beamflag); } if (sndg->sndg_edit == MB_YES) { if (output_bad == MB_YES && mb_beam_ok(sndg->sndg_beamflag) && flag ) { sndg->sndg_beamflag = MB_FLAG_FLAG + MB_FLAG_FILTER; files[sndg->sndg_file].nflagged++; } else if (output_good == MB_YES && !mb_beam_ok(sndg->sndg_beamflag) && sndg->sndg_beamflag != MB_FLAG_NULL && !flag) { sndg->sndg_beamflag = MB_FLAG_NONE; files[sndg->sndg_file].nunflagged++; } else if (output_good == MB_YES && !mb_beam_ok(sndg->sndg_beamflag) && sndg->sndg_beamflag != MB_FLAG_NULL && flag) { sndg->sndg_edit = MB_NO; } } /* print output debug statements */ if (verbose >= 2) { fprintf(stderr,"\ndbg2 MBIO function <%s> completed\n", function_name); fprintf(stderr,"dbg2 Return values:\n"); fprintf(stderr,"dbg2 sndg->sndg_edit: %d\n",sndg->sndg_edit); fprintf(stderr,"dbg2 sndg->sndg_beamflag: %d\n",sndg->sndg_beamflag); fprintf(stderr,"dbg2 error: %d\n",*error); fprintf(stderr,"dbg2 Return status:\n"); fprintf(stderr,"dbg2 status: %d\n",status); } /* return */ return(status); } /*--------------------------------------------------------------------*/