home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / gmt_customio.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-05-25  |  37.2 KB  |  1,166 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)gmt_customio.c    2.9  5/25/95
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /*
  8.  *
  9.  *    G M T _ C U S T O M I O . C   R O U T I N E S
  10.  *
  11.  * Takes care of all custom format gridfile input/output.  The
  12.  * industrious user may supply his/her own code to read specific data
  13.  * formats.  Four functions must be supplied, and they must conform
  14.  * to the GMT standard and return the same arguments as the generic
  15.  * GMT grdio functions.  See gmt_cdf.c for details.
  16.  *
  17.  * To add another data format:
  18.  *
  19.  *    1. Write the four required routines (see below).
  20.  *    2. increment parameter N_GRD_FORMATS in file gmt_grdio.h
  21.  *    3. Append another entry in the gmt_customio.h file.
  22.  *    4. Provide another entry in the lib/gmtformats.d file
  23.  *
  24.  * Author:    Paul Wessel
  25.  * Date:    9-SEP-1992
  26.  * Modified:    27-JUN-1995
  27.  * Version:    3.0
  28.  *
  29.  * Functions include:
  30.  *
  31.  *    *_read_grd_info :    Read header from file
  32.  *    *_read_grd :        Read header and data set from file
  33.  *    *_write_grd_info :    Update header in existing file
  34.  *    *_write_grd :        Write header and data set to new file
  35.  *
  36.  * where * is a prefix specific to a particular data format
  37.  *
  38.  * NOTE:  1. GMT assumes that read_grd_info has been called before calls
  39.  *         to read_grd.
  40.  *      2. Some formats may permit pipes to be used.  In that case GMT
  41.  *         expects the filename to be "=" (the equal sign).  It is the
  42.  *         responsibility of the custom routines to test for "=" and
  43.  *         give error messages if piping is not supported for this format
  44.  *         (e.g., netcdf uses fseek and can therefore not use pipes; other
  45.  *         formats may have similar limitations)
  46.  * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
  47.  
  48. #include "gmt.h"
  49.  
  50. int grdio_init () {
  51.     int id;
  52.  
  53.     /* FORMAT # 0: Default GMT netcdf-based grdio  LEAVE THIS AS IS! */
  54.     
  55.     gmtio_readinfo[0]  = (PFI) cdf_read_grd_info;
  56.     gmtio_writeinfo[0] = (PFI) cdf_write_grd_info;
  57.     gmtio_readgrd[0]   = (PFI) cdf_read_grd;
  58.     gmtio_writegrd[0]  = (PFI) cdf_write_grd;
  59.     
  60.     /*
  61.      * ----------------------------------------------
  62.      * ADD CUSTOM FORMATS BELOW AS THEY ARE NEEDED */
  63.     
  64.     
  65.     /* FORMAT # 1: GMT native binary (float) grdio */
  66.     
  67.     id = 1;
  68.     gmtio_readinfo[id]  = (PFI) bin_read_grd_info;
  69.     gmtio_writeinfo[id] = (PFI) bin_read_grd_info;
  70.     gmtio_readgrd[id]   = (PFI) bin_read_grd;
  71.     gmtio_writegrd[id]  = (PFI) bin_write_grd;
  72.     
  73.     /* FORMAT # 2: GMT native binary (short) grdio */
  74.     
  75.     id = 2;
  76.     gmtio_readinfo[id]  = (PFI) short_read_grd_info;
  77.     gmtio_writeinfo[id] = (PFI) short_read_grd_info;
  78.     gmtio_readgrd[id]   = (PFI) short_read_grd;
  79.     gmtio_writegrd[id]  = (PFI) short_write_grd;
  80.     
  81.     /* FORMAT # 3: SUN 8-bit standard rasterfile grdio */
  82.     
  83.     id = 3;
  84.     gmtio_readinfo[id]  = (PFI) ras_read_grd_info;
  85.     gmtio_writeinfo[id] = (PFI) ras_read_grd_info;
  86.     gmtio_readgrd[id]   = (PFI) ras_read_grd;
  87.     gmtio_writegrd[id]  = (PFI) ras_write_grd;
  88.     
  89.     /* FORMAT # 4: */
  90.     id = 4;
  91.  
  92. }
  93.  
  94. /* CUSTOM I/O FUNCTIONS FOR GRIDDED DATA FILES */
  95.  
  96. /*-----------------------------------------------------------
  97.  * Format # :    1
  98.  * Type :    Native binary (float) C file
  99.  * Prefix :    bin_
  100.  * Author :    Paul Wessel, SOEST
  101.  * Date :    27-JUN-1994
  102.  * Purpose:    The native binary output format is used
  103.  *        primarily for piped i/o between programs
  104.  *        that otherwise would use temporary, large
  105.  *        intermediate grdfiles.  Note that not all
  106.  *        programs can support piping (they may have
  107.  *        to re-read the file or accept more than one
  108.  *        grdfile).
  109.  * Functions :    bin_read_grd_info, bin_write_grd_info,
  110.  *        bin_read_grd, bin_write_grd
  111.  *-----------------------------------------------------------*/
  112.  
  113. int bin_read_grd_info (file, header)
  114. char *file;
  115. struct GRD_HEADER *header; {
  116.     FILE *fp;
  117.     
  118.     if (!strcmp (file, "="))
  119.         fp = stdin;
  120.     else if ((fp = fopen (file, "r")) == NULL) {
  121.         fprintf (stderr, "GMT Fatal Error: Could not open file %s!\n", file);
  122.         exit (-1);
  123.     }
  124.     
  125.     if (fread ((char *)header, sizeof (struct GRD_HEADER), 1, fp) != 1) {
  126.         fprintf (stderr, "GMT Fatal Error: Error reading file %s!\n", file);
  127.         exit (-1);
  128.     }
  129.  
  130.     if (fp != stdin) fclose (fp);
  131.     
  132.     return (FALSE);
  133. }
  134.  
  135. int bin_write_grd_info (file, header)
  136. char *file;
  137. struct GRD_HEADER *header; {
  138.     FILE *fp;
  139.     
  140.     if (!strcmp (file, "="))
  141.         fp = stdout;
  142.     else if ((fp = fopen (file, "w")) == NULL) {
  143.         fprintf (stderr, "GMT Fatal Error: Could not create file %s!\n", file);
  144.         exit (-1);
  145.     }
  146.     
  147.     if (fwrite ((char *)header, sizeof (struct GRD_HEADER), 1, fp) != 1) {
  148.         fprintf (stderr, "GMT Fatal Error: Error reading file %s!\n", file);
  149.         exit (-1);
  150.     }
  151.  
  152.     if (fp != stdout) fclose (fp);
  153.     
  154.     return (FALSE);
  155. }
  156.  
  157. int bin_read_grd (file, header, grid, w, e, s, n, pad, complex)
  158. char *file;
  159. struct GRD_HEADER *header;
  160. float *grid;
  161. double w, e, s, n;    /* Sub-region to extract  [Use entire file if 0,0,0,0] */
  162. int pad[];        /* padding (# of empty rows/columns) on w, e, s, n of grid, respectively */
  163. BOOLEAN complex;    /* TRUE if array is to hold real and imaginary parts */
  164. {
  165.     int first_col, last_col, first_row, last_row, nm, one_or_zero;
  166.     int i, j, j2, ij, width_in, width_out, height_in, i_0_out, inc = 1;
  167.     FILE *fp;
  168.     BOOLEAN piping = FALSE, geo = FALSE;
  169.     float *tmp;
  170.     double off, half_or_zero, x, small;
  171.     
  172.     if (!strcmp (file, "=")) {
  173.         fp = stdin;
  174.         piping = TRUE;
  175.     }
  176.     else if ((fp = fopen (file, "r")) != NULL)    /* Skip header */
  177.         fseek (fp, (long) sizeof (struct GRD_HEADER), 0);
  178.     else {
  179.         fprintf (stderr, "GMT Fatal Error: Could not open file %s!\n", file);
  180.         exit (-1);
  181.     }
  182.     
  183.     if (w == 0.0 && e == 0.0) {    /* Get entire file and return */
  184.         nm = header->nx * header->ny;
  185.         if (fread ((char *)grid, sizeof (float), nm, fp) != nm) {
  186.             fprintf (stderr, "GMT Fatal Error: Error reading file %s!\n", file);
  187.             exit (-1);
  188.         }
  189.         return (FALSE);
  190.     }
  191.  
  192.     /* Must deal with a subregion */
  193.  
  194.     if (w < header->x_min || e > header->x_max) geo = TRUE;    /* Dealing with periodic grid */
  195.     
  196.     one_or_zero = (header->node_offset) ? 0 : 1;
  197.     off = (header->node_offset) ? 0.0 : 0.5;
  198.     
  199.     /* Get dimension of subregion */
  200.     
  201.     width_in  = rint ((e - w) / header->x_inc) + one_or_zero;
  202.     height_in = rint ((n - s) / header->y_inc) + one_or_zero;
  203.     
  204.     /* Get first and last row and column numbers */
  205.     
  206.     small = 0.1 * header->x_inc;
  207.     first_col = floor ((w - header->x_min + small) / header->x_inc);
  208.     last_col  = ceil ((e - header->x_min - small) / header->x_inc) - 1 + one_or_zero;
  209.     small = 0.1 * header->y_inc;
  210.     first_row = floor ((header->y_max - n + small) / header->y_inc);
  211.     last_row  = ceil ((header->y_max - s - small) / header->y_inc) - 1 + one_or_zero;
  212.  
  213.     if ((last_col - first_col + 1) > width_in) last_col--;
  214.     if ((last_row - first_row + 1) > height_in) last_row--;
  215.     if ((last_col - first_col + 1) > width_in) first_col++;
  216.     if ((last_row - first_row + 1) > height_in) first_row++;
  217.     
  218.     width_out = width_in;        /* Width of output array */
  219.     if (pad[0] > 0) width_out += pad[0];
  220.     if (pad[1] > 0) width_out += pad[1];
  221.     
  222.     i_0_out = pad[0];        /* Edge offset in output */
  223.     if (complex) {    /* Need twice as much space and load every 2nd cell */
  224.         width_out *= 2;
  225.         i_0_out *= 2;
  226.         inc = 2;
  227.     }
  228.     
  229.     if (geo) {    /* Must rollover in longitudes */
  230.         int *k;
  231.         tmp = (float *) memory (CNULL, header->nx, sizeof (float), "bin_read_grd");
  232.         k = (int *) memory (CNULL, width_in, sizeof (int), "bin_read_grd");
  233.         
  234.         half_or_zero = (header->node_offset) ? 0.5 : 0.0;
  235.         small = 0.1 * header->x_inc;    /* Anything smaller than 0.5 dx will do */
  236.         for (i = 0; i < width_in; i++) {
  237.             x = w + (i + half_or_zero) * header->x_inc;
  238.             if ((header->x_min - x) > small)
  239.                 x += 360.0;
  240.             else if ((x - header->x_max) > small)
  241.                 x -= 360.0;
  242.             k[i] = (int) floor (((x - header->x_min) / header->x_inc) + off);
  243.         }
  244.         if (piping)    /* Skip data by reading it */
  245.             for (j = 0; j < first_row; j++) fread ((char *) tmp, sizeof (float), header->nx, fp);
  246.         else /* Simply seek by it */
  247.             fseek (fp, (long) (first_row * header->nx * sizeof (float)), 1);
  248.             
  249.         for (j = first_row, j2 = 0; j <= last_row; j++, j2++) {
  250.             fread ((char *) tmp, sizeof (float), header->nx, fp);    /* Get one row */
  251.             ij = (j2 + pad[3]) * width_out + i_0_out;
  252.             for (i = 0; i < width_in; i++) grid[ij+i*inc] = tmp[k[i]];
  253.         }
  254.         if (piping)    /* Skip data by reading it */
  255.             for (j = last_row + 1; j < header->ny; j++) fread ((char *) tmp, sizeof (float), header->nx, fp);
  256.             
  257.         free ((char *)k);
  258.     }
  259.     else {    /* A bit easier here */
  260.         tmp = (float *) memory (CNULL, width_in, sizeof (float), "bin_read_grd");
  261.         if (piping)    /* Skip data by reading it */
  262.             for (j = 0; j < first_row; j++) fread ((char *) tmp, sizeof (float), header->nx, fp);
  263.         else /* Simply seek by it */
  264.             fseek (fp, (long) (first_row * header->nx * sizeof (float)), 1);
  265.         for (j = first_row, j2 = 0; j <= last_row; j++, j2++) {
  266.             ij = (j2 + pad[3]) * width_out + i_0_out;
  267.             fread ((char *) tmp, sizeof (float), header->nx, fp);
  268.             if (complex)
  269.                 for (i = 0; i < width_in; i++) grid[ij+2*i] = tmp[first_col+i];
  270.             else
  271.                 for (i = 0; i < width_in; i++) grid[ij+i] = tmp[first_col+i];
  272.         }
  273.         if (piping)    /* Skip data by reading it */
  274.             for (j = last_row + 1; j < header->ny; j++) fread ((char *) tmp, sizeof (float), header->nx, fp);
  275.     }
  276.     
  277.     header->nx = width_in;
  278.     header->ny = height_in;
  279.     header->x_min = w;
  280.     header->x_max = e;
  281.     header->y_min = s;
  282.     header->y_max = n;
  283.  
  284.     header->z_min = MAXDOUBLE;    header->z_max = -MAXDOUBLE;
  285.     for (j = 0; j < header->ny; j++) {
  286.         for (i = 0; i < header->nx; i++) {
  287.             j2 = (complex) ? 2 * (i + pad[0]) : i + pad[0];
  288.             ij = (j + pad[3]) * width_out + j2;
  289.             if (bad_float ((double)grid[ij])) continue;
  290.             header->z_min = MIN (header->z_min, grid[ij]);
  291.             header->z_max = MAX (header->z_max, grid[ij]);
  292.         }
  293.     }
  294.     if (fp != stdin) fclose (fp);
  295.     
  296.     free ((char *)tmp);
  297.     return (FALSE);
  298. }
  299.  
  300. int bin_write_grd (file, header, grid, w, e, s, n, pad, complex)
  301. char *file;
  302. struct GRD_HEADER *header;
  303. float *grid;
  304. double w, e, s, n;    /* Sub-region to write out */
  305. int pad[];        /* padding on w, e, s, n of grid, respectively */
  306. BOOLEAN complex;    /* TRUE if array is to hold real and imaginary parts */
  307. {
  308.     int i, i2, nm, inc = 1, *k;
  309.     int j, ij, j2, width_in, width_out, height_out, one_or_zero;
  310.     int first_col, last_col, first_row, last_row;
  311.     
  312.     BOOLEAN geo = FALSE;
  313.     double off, half_or_zero, small, x;
  314.     
  315.     float *tmp;
  316.     
  317.     FILE *fp;
  318.     
  319.     if (!strcmp (file, "=")) {
  320.         fp = stdout;
  321.     }
  322.     else if ((fp = fopen (file, "w")) == NULL) {
  323.         fprintf (stderr, "GMT Fatal Error: Could not create file %s!\n", file);
  324.         exit (-1);
  325.     }
  326.     
  327.     if (w == 0.0 && e == 0.0) {    /* Write entire file and return */
  328.         /* Find min/max of data */
  329.         
  330.         nm = header->nx * header->ny;
  331.         header->z_min = MAXDOUBLE;    header->z_max = -MAXDOUBLE;
  332.         
  333.         for (i = 0; i < nm; i++) {
  334.                 ij = (complex) ? 2 * i : i;
  335.                 if (bad_float ((double)grid[ij])) continue;
  336.                 header->z_min = MIN (header->z_min, grid[ij]);
  337.                 header->z_max = MAX (header->z_max, grid[ij]);
  338.         }
  339.         
  340.         if (fwrite ((char *)header, sizeof (struct GRD_HEADER), 1, fp) != 1) {
  341.             fprintf (stderr, "GMT Fatal Error: Error writing file %s!\n", file);
  342.             exit (-1);
  343.         }
  344.  
  345.         nm = header->nx * header->ny;
  346.         if (fwrite ((char *)grid, sizeof (float), nm, fp) != nm) {
  347.             fprintf (stderr, "GMT Fatal Error: Error writing file %s!\n", file);
  348.             exit (-1);
  349.         }
  350.         return (FALSE);
  351.     }
  352.  
  353.     if (w < header->x_min || e > header->x_max) geo = TRUE;    /* Dealing with periodic grid */
  354.     
  355.     one_or_zero = (header->node_offset) ? 0 : 1;
  356.     off = (header->node_offset) ? 0.0 : 0.5;
  357.     
  358.     /* Get dimension of subregion to write */
  359.     
  360.     width_out  = rint ((e - w) / header->x_inc) + one_or_zero;
  361.     height_out = rint ((n - s) / header->y_inc) + one_or_zero;
  362.     
  363.     /* Get first and last row and column numbers */
  364.     
  365.     small = 0.1 * header->x_inc;
  366.     first_col = floor ((w - header->x_min + small) / header->x_inc);
  367.     last_col  = ceil ((e - header->x_min - small) / header->x_inc) -1 + one_or_zero;
  368.     small = 0.1 * header->y_inc;
  369.     first_row = floor ((header->y_max - n + small) / header->y_inc);
  370.     last_row  = ceil ((header->y_max - s - small) / header->y_inc) -1 + one_or_zero;
  371.     
  372.     if ((last_col - first_col + 1) > width_out) last_col--;
  373.     if ((last_row - first_row + 1) > height_out) last_row--;
  374.     if ((last_col - first_col + 1) > width_out) first_col++;
  375.     if ((last_row - first_row + 1) > height_out) first_row++;
  376.  
  377.     width_in = width_out;        /* Physical width of input array */
  378.     if (pad[0] > 0) width_in += pad[0];
  379.     if (pad[1] > 0) width_in += pad[1];
  380.     
  381.     if (complex) inc = 2;
  382.  
  383.     header->x_min = w;
  384.     header->x_max = e;
  385.     header->y_min = s;
  386.     header->y_max = n;
  387.     
  388.     /* Find xmin/zmax */
  389.     
  390.     header->z_min = MAXDOUBLE;    header->z_max = -MAXDOUBLE;
  391.     for (j = first_row, j2 = pad[3]; j <= last_row; j++, j2++) {
  392.         for (i = first_col, i2 = pad[0]; i <= last_col; i++, i2++) {
  393.             ij = (j2 * width_in + i2) * inc;
  394.             if (bad_float ((double)grid[ij])) continue;
  395.             header->z_min = MIN (header->z_min, grid[ij]);
  396.             header->z_max = MAX (header->z_max, grid[ij]);
  397.         }
  398.     }
  399.     
  400.     /* store header information and array */
  401.     
  402.     if (fwrite ((char *)header, sizeof (struct GRD_HEADER), 1, fp) != 1) {
  403.         fprintf (stderr, "GMT Fatal Error: Error writing file %s!\n", file);
  404.         exit (-1);
  405.     }
  406.  
  407.     if (geo) {
  408.         tmp = (float *) memory (CNULL, width_in, sizeof (float), "bin_write_grd");
  409.         k = (int *) memory (CNULL, width_out, sizeof (int), "bin_write_grd");
  410.         
  411.         half_or_zero = (header->node_offset) ? 0.5 : 0.0;
  412.         small = 0.1 * header->x_inc;    /* Anything smaller than 0.5 dx will do */
  413.         for (i = 0; i < width_out; i++) {
  414.             x = w + (i + half_or_zero) * header->x_inc;
  415.             if ((header->x_min - x) > small)
  416.                 x += 360.0;
  417.             else if ((x - header->x_max) > small)
  418.                 x -= 360.0;
  419.             k[i] = (int) floor (((x - header->x_min) / header->x_inc) + off);
  420.         }
  421.         i2 = first_col + pad[0];
  422.         for (j = 0, j2 = first_row + pad[3]; j < height_out; j++, j2++) {
  423.             ij = j2 * width_in + i2;
  424.             for (i = 0; i < width_out; i++) tmp[i] = grid[inc * (ij+k[i])];
  425.             fwrite ((char *)tmp, sizeof (float), width_out, fp);
  426.         }
  427.         free ((char *)k);
  428.     }
  429.     else {
  430.         if (complex) tmp = (float *) memory (CNULL, width_in, sizeof (float), "bin_write_grd");
  431.         i2 = first_col + pad[0];
  432.         for (j = 0, j2 = first_row + pad[3]; j < height_out; j++, j2++) {
  433.             ij = j2 * width_in + i2;
  434.             if (complex) {
  435.                 for (i = 0; i < width_out; i++) tmp[i] = grid[2*(ij+i)];
  436.                 fwrite ((char *)tmp, sizeof (float), width_out, fp);
  437.             }
  438.             else
  439.                 fwrite ((char *)&grid[ij], sizeof (float), width_out, fp);
  440.         }
  441.     }
  442.     if (fp != stdout) fclose (fp);
  443.     
  444.     if (complex || geo) free ((char *)tmp);
  445.     
  446.     return (FALSE);
  447.  
  448. }
  449.  
  450. /*-----------------------------------------------------------
  451.  * Format # :    2
  452.  * Type :    Native binary (short) C file
  453.  * Prefix :    bin_
  454.  * Author :    Paul Wessel, SOEST
  455.  * Date :    27-JUN-1994
  456.  * Purpose:    The native binary output format is used
  457.  *        primarily for piped i/o between programs
  458.  *        that otherwise would use temporary, large
  459.  *        intermediate grdfiles.  Note that not all
  460.  *        programs can support piping (they may have
  461.  *        to re-read the file or accept more than one
  462.  *        grdfile).  Datasets with limited range may
  463.  *        be stored using 2-byte integers which will
  464.  *        reduce storage space and improve throughput.
  465.  * Functions :    short_read_grd_info, short_write_grd_info,
  466.  *        short_read_grd, short_write_grd
  467.  *-----------------------------------------------------------*/
  468.  
  469. int short_read_grd_info (file, header)
  470. char *file;
  471. struct GRD_HEADER *header; {    
  472.     return (bin_read_grd_info (file, header));
  473. }
  474.  
  475. int short_write_grd_info (file, header)
  476. char *file;
  477. struct GRD_HEADER *header; {
  478.     return (bin_write_grd_info (file, header));
  479. }
  480.  
  481. int short_read_grd (file, header, grid, w, e, s, n, pad, complex)
  482. char *file;
  483. struct GRD_HEADER *header;
  484. float *grid;
  485. double w, e, s, n;    /* Sub-region to extract  [Use entire file if 0,0,0,0] */
  486. int pad[];        /* padding (# of empty rows/columns) on w, e, s, n of grid, respectively */
  487. BOOLEAN complex;    /* TRUE if array is to hold real and imaginary parts */
  488. {
  489.     int first_col, last_col, first_row, last_row, nm, kk, one_or_zero;
  490.     int i, j, j2, ij, width_in, width_out, height_in, i_0_out, inc = 1;
  491.     FILE *fp;
  492.     BOOLEAN piping = FALSE, geo = FALSE, check = FALSE;
  493.     short *tmp;
  494.     double off, half_or_zero, x, small;
  495.     
  496.     if (complex) {
  497.         fprintf (stderr, "GMT Fatal Error: short int grdfile %s cannot hold complex data!\n", file);
  498.         exit (-1);
  499.     }
  500.     if (!strcmp (file, "=")) {
  501.         fp = stdin;
  502.         piping = TRUE;
  503.     }
  504.     else if ((fp = fopen (file, "r")) != NULL)    /* Skip header */
  505.         fseek (fp, (long) sizeof (struct GRD_HEADER), 0);
  506.     else {
  507.         fprintf (stderr, "GMT Fatal Error: Could not open file %s!\n", file);
  508.         exit (-1);
  509.     }
  510.     
  511.     check = !bad_float (grd_in_nan_value);
  512.     
  513.     if (w == 0.0 && e == 0.0) {    /* Get entire file and return */
  514.         nm = header->nx * header->ny;
  515.         tmp = (short *) memory (CNULL, nm, sizeof (short), "short_read_grd");
  516.         if (fread ((char *)tmp, sizeof (short), nm, fp) != nm) {
  517.             fprintf (stderr, "GMT Fatal Error: Error reading file %s!\n", file);
  518.             exit (-1);
  519.         }
  520.         for (i = 0; i < nm; i++) {
  521.             grid[i] = tmp[i];
  522.             if (check && grid[i] == grd_in_nan_value) grid[i] = gmt_NaN;
  523.         }
  524.         free ((char *)tmp);
  525.         return (FALSE);
  526.     }
  527.  
  528.     /* Must deal with a subregion */
  529.  
  530.     if (w < header->x_min || e > header->x_max) geo = TRUE;    /* Dealing with periodic grid */
  531.     
  532.     one_or_zero = (header->node_offset) ? 0 : 1;
  533.     off = (header->node_offset) ? 0.0 : 0.5;
  534.     
  535.     /* Get dimension of subregion */
  536.     
  537.     width_in  = rint ((e - w) / header->x_inc) + one_or_zero;
  538.     height_in = rint ((n - s) / header->y_inc) + one_or_zero;
  539.     
  540.     /* Get first and last row and column numbers */
  541.     
  542.     small = 0.1 * header->x_inc;
  543.     first_col = floor ((w - header->x_min + small) / header->x_inc);
  544.     last_col  = ceil ((e - header->x_min - small) / header->x_inc) - 1 + one_or_zero;
  545.     small = 0.1 * header->y_inc;
  546.     first_row = floor ((header->y_max - n + small) / header->y_inc);
  547.     last_row  = ceil ((header->y_max - s - small) / header->y_inc) - 1 + one_or_zero;
  548.  
  549.     if ((last_col - first_col + 1) > width_in) last_col--;
  550.     if ((last_row - first_row + 1) > height_in) last_row--;
  551.     if ((last_col - first_col + 1) > width_in) first_col++;
  552.     if ((last_row - first_row + 1) > height_in) first_row++;
  553.     
  554.     width_out = width_in;        /* Width of output array */
  555.     if (pad[0] > 0) width_out += pad[0];
  556.     if (pad[1] > 0) width_out += pad[1];
  557.     
  558.     i_0_out = pad[0];        /* Edge offset in output */
  559.     
  560.     if (geo) {    /* Must rollover in longitudes */
  561.         int *k;
  562.         tmp = (short *) memory (CNULL, header->nx, sizeof (short), "short_read_grd");
  563.         k = (int *) memory (CNULL, width_in, sizeof (int), "short_read_grd");
  564.         
  565.         half_or_zero = (header->node_offset) ? 0.5 : 0.0;
  566.         small = 0.1 * header->x_inc;    /* Anything smaller than 0.5 dx will do */
  567.         for (i = 0; i < width_in; i++) {
  568.             x = w + (i + half_or_zero) * header->x_inc;
  569.             if ((header->x_min - x) > small)
  570.                 x += 360.0;
  571.             else if ((x - header->x_max) > small)
  572.                 x -= 360.0;
  573.             k[i] = (int) floor (((x - header->x_min) / header->x_inc) + off);
  574.         }
  575.         if (piping)    /* Skip data by reading it */
  576.             for (j = 0; j < first_row; j++) fread ((char *) tmp, sizeof (short), header->nx, fp);
  577.         else /* Simply seek by it */
  578.             fseek (fp, (long) (first_row * header->nx * sizeof (short)), 1);
  579.             
  580.         for (j = first_row, j2 = 0; j <= last_row; j++, j2++) {
  581.             fread ((char *) tmp, sizeof (short), header->nx, fp);    /* Get one row */
  582.             ij = (j2 + pad[3]) * width_out + i_0_out;
  583.             for (i = 0; i < width_in; i++) {
  584.                 kk = ij+i*inc;
  585.                 grid[kk] = tmp[k[i]];
  586.                 if (check && grid[kk] == grd_in_nan_value) grid[kk] = gmt_NaN;
  587.             }
  588.         }
  589.         if (piping)    /* Skip data by reading it */
  590.             for (j = last_row + 1; j < header->ny; j++) fread ((char *) tmp, sizeof (short), header->nx, fp);
  591.             
  592.         free ((char *)k);
  593.     }
  594.     else {    /* A bit easier here */
  595.         tmp = (short *) memory (CNULL, width_in, sizeof (short), "short_read_grd");
  596.         if (piping)    /* Skip data by reading it */
  597.             for (j = 0; j < first_row; j++) fread ((char *) tmp, sizeof (short), header->nx, fp);
  598.         else /* Simply seek by it */
  599.             fseek (fp, (long) (first_row * header->nx * sizeof (short)), 1);
  600.         for (j = first_row, j2 = 0; j <= last_row; j++, j2++) {
  601.             ij = (j2 + pad[3]) * width_out + i_0_out;
  602.             fread ((char *) tmp, sizeof (short), header->nx, fp);
  603.             for (i = 0; i < width_in; i++) {
  604.                 kk = ij+i;
  605.                 grid[kk] = tmp[first_col+i];
  606.                 if (check && grid[kk] == grd_in_nan_value) grid[kk] = gmt_NaN;
  607.             }
  608.         }
  609.         if (piping)    /* Skip data by reading it */
  610.             for (j = last_row + 1; j < header->ny; j++) fread ((char *) tmp, sizeof (short), header->nx, fp);
  611.     }
  612.     
  613.     header->nx = width_in;
  614.     header->ny = height_in;
  615.     header->x_min = w;
  616.     header->x_max = e;
  617.     header->y_min = s;
  618.     header->y_max = n;
  619.  
  620.     header->z_min = MAXDOUBLE;    header->z_max = -MAXDOUBLE;
  621.     for (j = 0; j < header->ny; j++) {
  622.         for (i = 0; i < header->nx; i++) {
  623.             ij = (j + pad[3]) * width_out + i + pad[0];
  624.             if (bad_float ((double) grid[ij])) continue;
  625.             header->z_min = MIN (header->z_min, grid[ij]);
  626.             header->z_max = MAX (header->z_max, grid[ij]);
  627.         }
  628.     }
  629.     if (fp != stdin) fclose (fp);
  630.     
  631.     free ((char *)tmp);
  632.     return (FALSE);
  633. }
  634.  
  635. int short_write_grd (file, header, grid, w, e, s, n, pad, complex)
  636. char *file;
  637. struct GRD_HEADER *header;
  638. float *grid;
  639. double w, e, s, n;    /* Sub-region to write out */
  640. int pad[];        /* padding on w, e, s, n of grid, respectively */
  641. BOOLEAN complex;    /* TRUE if array is to hold real and imaginary parts */
  642. {
  643.     int i, i2, nm, kk, *k;
  644.     int j, ij, j2, width_in, width_out, height_out, one_or_zero;
  645.     int first_col, last_col, first_row, last_row, check = FALSE;
  646.     
  647.     BOOLEAN geo = FALSE;
  648.     double off, half_or_zero, small, x;
  649.     
  650.     short *tmp;
  651.     
  652.     FILE *fp;
  653.     
  654.     if (complex) {
  655.         fprintf (stderr, "GMT Fatal Error: short int grdfile %s cannot hold complex data!\n", file);
  656.         exit (-1);
  657.     }
  658.     if (!strcmp (file, "=")) {
  659.         fp = stdout;
  660.     }
  661.     else if ((fp = fopen (file, "w")) == NULL) {
  662.         fprintf (stderr, "GMT Fatal Error: Could not create file %s!\n", file);
  663.         exit (-1);
  664.     }
  665.     
  666.     check = !bad_float (grd_out_nan_value);
  667.  
  668.     if (w == 0.0 && e == 0.0) {    /* Write entire file and return */
  669.         /* Find min/max of data */
  670.         
  671.         nm = header->nx * header->ny;
  672.         header->z_min = MAXDOUBLE;    header->z_max = -MAXDOUBLE;
  673.         for (i = 0; i < nm; i++) {
  674.                 ij = (complex) ? 2 * i : i;
  675.                 if (bad_float ((double)grid[ij])) continue;
  676.                 header->z_min = MIN (header->z_min, grid[ij]);
  677.                 header->z_max = MAX (header->z_max, grid[ij]);
  678.         }
  679.         if (fwrite ((char *)header, sizeof (struct GRD_HEADER), 1, fp) != 1) {
  680.             fprintf (stderr, "GMT Fatal Error: Error writing file %s!\n", file);
  681.             exit (-1);
  682.         }
  683.  
  684.         tmp = (short *) memory (CNULL, nm, sizeof (short), "short_read_grd");
  685.         for (i = 0; i < nm; i++) {
  686.             if (check && bad_float ((double)grid[i])) grid[i] = grd_out_nan_value;
  687.             tmp[i] = rint ((double)grid[i]);
  688.         }
  689.         if (fwrite ((char *)tmp, sizeof (short), nm, fp) != nm) {
  690.             fprintf (stderr, "GMT Fatal Error: Error writing file %s!\n", file);
  691.             exit (-1);
  692.         }
  693.         free ((char *)tmp);
  694.         return (FALSE);
  695.     }
  696.  
  697.     if (w < header->x_min || e > header->x_max) geo = TRUE;    /* Dealing with periodic grid */
  698.     
  699.     one_or_zero = (header->node_offset) ? 0 : 1;
  700.     off = (header->node_offset) ? 0.0 : 0.5;
  701.     
  702.     /* Get dimension of subregion to write */
  703.     
  704.     width_out  = rint ((e - w) / header->x_inc) + one_or_zero;
  705.     height_out = rint ((n - s) / header->y_inc) + one_or_zero;
  706.     
  707.     /* Get first and last row and column numbers */
  708.     
  709.     small = 0.1 * header->x_inc;
  710.     first_col = floor ((w - header->x_min + small) / header->x_inc);
  711.     last_col  = ceil ((e - header->x_min - small) / header->x_inc) -1 + one_or_zero;
  712.     small = 0.1 * header->y_inc;
  713.     first_row = floor ((header->y_max - n + small) / header->y_inc);
  714.     last_row  = ceil ((header->y_max - s - small) / header->y_inc) -1 + one_or_zero;
  715.     
  716.     if ((last_col - first_col + 1) > width_out) last_col--;
  717.     if ((last_row - first_row + 1) > height_out) last_row--;
  718.     if ((last_col - first_col + 1) > width_out) first_col++;
  719.     if ((last_row - first_row + 1) > height_out) first_row++;
  720.  
  721.     width_in = width_out;        /* Physical width of input array */
  722.     if (pad[0] > 0) width_in += pad[0];
  723.     if (pad[1] > 0) width_in += pad[1];
  724.     
  725.     header->x_min = w;
  726.     header->x_max = e;
  727.     header->y_min = s;
  728.     header->y_max = n;
  729.     
  730.     /* Find xmin/zmax */
  731.     
  732.     header->z_min = MAXDOUBLE;    header->z_max = -MAXDOUBLE;
  733.     for (j = first_row, j2 = pad[3]; j <= last_row; j++, j2++) {
  734.         for (i = first_col, i2 = pad[0]; i <= last_col; i++, i2++) {
  735.             ij = j2 * width_in + i2;
  736.             if (bad_float ((double)grid[ij])) continue;
  737.             header->z_min = MIN (header->z_min, grid[ij]);
  738.             header->z_max = MAX (header->z_max, grid[ij]);
  739.         }
  740.     }
  741.     
  742.     /* store header information and array */
  743.     
  744.     if (fwrite ((char *)header, sizeof (struct GRD_HEADER), 1, fp) != 1) {
  745.         fprintf (stderr, "GMT Fatal Error: Error writing file %s!\n", file);
  746.         exit (-1);
  747.     }
  748.  
  749.     tmp = (short *) memory (CNULL, width_in, sizeof (short), "short_write_grd");
  750.     if (geo) {
  751.         k = (int *) memory (CNULL, width_out, sizeof (int), "short_write_grd");
  752.         
  753.         half_or_zero = (header->node_offset) ? 0.5 : 0.0;
  754.         small = 0.1 * header->x_inc;    /* Anything smaller than 0.5 dx will do */
  755.         for (i = 0; i < width_out; i++) {
  756.             x = w + (i + half_or_zero) * header->x_inc;
  757.             if ((header->x_min - x) > small)
  758.                 x += 360.0;
  759.             else if ((x - header->x_max) > small)
  760.                 x -= 360.0;
  761.             k[i] = (int) floor (((x - header->x_min) / header->x_inc) + off);
  762.         }
  763.         i2 = first_col + pad[0];
  764.         for (j = 0, j2 = first_row + pad[3]; j < height_out; j++, j2++) {
  765.             ij = j2 * width_in + i2;
  766.             for (i = 0; i < width_out; i++) {
  767.                 kk = ij+k[i];
  768.                 if (check && bad_float ((double)grid[kk])) grid[kk] = grd_out_nan_value;
  769.                 tmp[i] = rint ((double)grid[kk]);
  770.             }
  771.             fwrite ((char *)tmp, sizeof (short), width_out, fp);
  772.         }
  773.         free ((char *)k);
  774.     }
  775.     else {
  776.         i2 = first_col + pad[0];
  777.         for (j = 0, j2 = first_row + pad[3]; j < height_out; j++, j2++) {
  778.             ij = j2 * width_in + i2;
  779.             for (i = 0; i < width_out; i++) {
  780.                 kk = ij+i;
  781.                 if (check && bad_float ((double)grid[kk])) grid[kk] = grd_out_nan_value;
  782.                 tmp[i] = rint ((double)grid[kk]);
  783.             }
  784.             fwrite ((char *)tmp, sizeof (short), width_out, fp);
  785.         }
  786.     }
  787.     if (fp != stdout) fclose (fp);
  788.     
  789.     free ((char *)tmp);
  790.     
  791.     return (FALSE);
  792.  
  793. }
  794.  
  795. /*-----------------------------------------------------------
  796.  * Format # :    3
  797.  * Type :    Standard 8-bit Sun rasterfiles
  798.  * Prefix :    ras_
  799.  * Author :    Paul Wessel, SOEST
  800.  * Date :    4-DEC-1994
  801.  * Purpose:    Rasterfiles may often be used to store
  802.  *        datasets of limited dynamic range where
  803.  *        8-bits is all that is needed.  Since the
  804.  *        usual information like w,e,s,n is not part
  805.  *        of a rasterfile's header, we assign default
  806.  *        values as follows:
  807.  *            w = s = 0.
  808.  *            e = ras_width;
  809.  *            n = ras_height;
  810.  *            dx = dy = 1
  811.  *        Such files are always pixel registrered
  812.  *            
  813.  * Functions :    ras_read_grd_info, ras_write_grd_info,
  814.  *        ras_read_grd, ras_write_grd
  815.  *-----------------------------------------------------------*/
  816.  
  817. struct rasterfile {
  818.     int  ras_magic;
  819.     int  ras_width;
  820.     int  ras_height;
  821.     int  ras_depth;
  822.     int  ras_length;
  823.     int  ras_type;
  824.     int  ras_maptype;
  825.     int  ras_maplength;
  826. };
  827.  
  828. int ras_read_grd_info (file, header)
  829. char *file;
  830. struct GRD_HEADER *header; {
  831.     FILE *fp;
  832.     struct rasterfile h;
  833.     unsigned char u;
  834.     int i;
  835.     
  836.     if (!strcmp (file, "="))
  837.         fp = stdin;
  838.     else if ((fp = fopen (file, "r")) == NULL) {
  839.         fprintf (stderr, "GMT Fatal Error: Could not open file %s!\n", file);
  840.         exit (-1);
  841.     }
  842.     
  843.     if (fread ((char *)&h, sizeof (struct rasterfile), 1, fp) != 1) {
  844.         fprintf (stderr, "GMT Fatal Error: Error reading file %s!\n", file);
  845.         exit (-1);
  846.     }
  847.     if (h.ras_type != 1 || h.ras_depth != 8) {
  848.         fprintf (stderr, "GMT Fatal Error: file %s not 8-bit standard Sun rasterfile!\n", file);
  849.         exit (-1);
  850.     }
  851.  
  852.     for (i = 0; i < h.ras_maplength; i++) fread ((char *)&u, sizeof (unsigned char *), 1, fp);    /* Skip colormap */
  853.     
  854.     if (fp != stdin) fclose (fp);
  855.     
  856.     grd_init (header, 0, (char **)NULL, FALSE);
  857.     
  858.     header->x_min = header->y_min = 0.0;
  859.     header->x_max = header->nx = h.ras_width;
  860.     header->y_max = header->ny = h.ras_height;
  861.     header->x_inc = header->y_inc = 1.0;
  862.     header->node_offset = 1;
  863.     
  864.     return (FALSE);
  865. }
  866.  
  867. int ras_write_grd_info (file, header)
  868. char *file;
  869. struct GRD_HEADER *header; {
  870.     FILE *fp;
  871.     struct rasterfile h;
  872.     
  873.     if (!strcmp (file, "="))
  874.         fp = stdout;
  875.     else if ((fp = fopen (file, "w")) == NULL) {
  876.         fprintf (stderr, "GMT Fatal Error: Could not create file %s!\n", file);
  877.         exit (-1);
  878.     }
  879.     
  880.     h.ras_magic = 0x59a66a95;
  881.     h.ras_width = header->nx;
  882.     h.ras_height = header->ny;
  883.     h.ras_depth = 8;
  884.     h.ras_length = header->ny * (int) ceil (header->nx/2.0) * 2;
  885.     h.ras_type = 1;
  886.     h.ras_maptype = 0;
  887.     h.ras_maplength = 0;
  888.     
  889.     if (fwrite ((char *)&h, sizeof (struct rasterfile), 1, fp) != 1) {
  890.         fprintf (stderr, "GMT Fatal Error: Error reading file %s!\n", file);
  891.         exit (-1);
  892.     }
  893.  
  894.     if (fp != stdout) fclose (fp);
  895.     
  896.     return (FALSE);
  897. }
  898.  
  899. int ras_read_grd (file, header, grid, w, e, s, n, pad, complex)
  900. char *file;
  901. struct GRD_HEADER *header;
  902. float *grid;
  903. double w, e, s, n;    /* Sub-region to extract  [Use entire file if 0,0,0,0] */
  904. int pad[];        /* padding (# of empty rows/columns) on w, e, s, n of grid, respectively */
  905. BOOLEAN complex;    /* Must be FALSE for rasterfiles */
  906. {
  907.     int first_col, last_col, first_row, last_row, one_or_zero;
  908.     int i, j, j2, ij, width_in, width_out, height_in, i_0_out, n2;
  909.     FILE *fp;
  910.     BOOLEAN piping = FALSE, check;
  911.     unsigned char *tmp;
  912.     double small;
  913.     struct rasterfile h;
  914.     
  915.     if (complex) {
  916.         fprintf (stderr, "GMT Fatal Error: Rasterfile %s cannot hold complex data!\n", file);
  917.         exit (-1);
  918.     }
  919.     if (!strcmp (file, "=")) {
  920.         fp = stdin;
  921.         piping = TRUE;
  922.     }
  923.     else if ((fp = fopen (file, "r")) != NULL) {    /* Skip header */
  924.         if (fread ((char *)&h, sizeof (struct rasterfile), 1, fp) != 1) {
  925.             fprintf (stderr, "GMT Fatal Error: Error reading file %s!\n", file);
  926.             exit (-1);
  927.         }
  928.         if (h.ras_maplength) fseek (fp, (long) h.ras_maplength, 1);
  929.     }
  930.     else {
  931.         fprintf (stderr, "GMT Fatal Error: Could not open file %s!\n", file);
  932.         exit (-1);
  933.     }
  934.     
  935.     n2 = (int) ceil (header->nx / 2.0) * 2;
  936.     tmp = (unsigned char *) memory (CNULL, n2, sizeof (unsigned char), "ras_read_grd");
  937.  
  938.     header->z_min = MAXDOUBLE;    header->z_max = -MAXDOUBLE;
  939.  
  940.     check = !bad_float (grd_in_nan_value);
  941.  
  942.     if (w == 0.0 && e == 0.0) {    /* Get entire file and return */
  943.         for (j = ij = 0; j < header->ny; j++) {
  944.             if (fread ((char *)tmp, sizeof (unsigned char), n2, fp) != n2) {
  945.                 fprintf (stderr, "GMT Fatal Error: Error reading file %s!\n", file);
  946.                 exit (-1);
  947.             }
  948.             for (i = 0; i < header->nx; i++, ij++) {
  949.                 grid[ij] = tmp[i];
  950.                 if (check && grid[ij] == grd_in_nan_value) grid[ij] = gmt_NaN;
  951.                 if (bad_float ((double)grid[ij])) continue;
  952.                 header->z_min = MIN (header->z_min, grid[ij]);
  953.                 header->z_max = MAX (header->z_max, grid[ij]);
  954.             }
  955.         }
  956.         free ((char *)tmp);
  957.         return (FALSE);
  958.     }
  959.  
  960.     /* Must deal with a subregion */
  961.  
  962.     one_or_zero = (header->node_offset) ? 0 : 1;
  963.     
  964.     /* Get dimension of subregion */
  965.     
  966.     width_in  = rint ((e - w) / header->x_inc) + one_or_zero;
  967.     height_in = rint ((n - s) / header->y_inc) + one_or_zero;
  968.     
  969.     /* Get first and last row and column numbers */
  970.     
  971.     small = 0.1 * header->x_inc;
  972.     first_col = floor ((w - header->x_min + small) / header->x_inc);
  973.     last_col  = ceil ((e - header->x_min - small) / header->x_inc) - 1 + one_or_zero;
  974.     small = 0.1 * header->y_inc;
  975.     first_row = floor ((header->y_max - n + small) / header->y_inc);
  976.     last_row  = ceil ((header->y_max - s - small) / header->y_inc) - 1 + one_or_zero;
  977.  
  978.     if ((last_col - first_col + 1) > width_in) last_col--;
  979.     if ((last_row - first_row + 1) > height_in) last_row--;
  980.     if ((last_col - first_col + 1) > width_in) first_col++;
  981.     if ((last_row - first_row + 1) > height_in) first_row++;
  982.     
  983.     width_out = width_in;        /* Width of output array */
  984.     if (pad[0] > 0) width_out += pad[0];
  985.     if (pad[1] > 0) width_out += pad[1];
  986.     
  987.     i_0_out = pad[0];        /* Edge offset in output */
  988.     
  989.     if (piping)    /* Skip data by reading it */
  990.         for (j = 0; j < first_row; j++) fread ((char *) tmp, sizeof (unsigned char), n2, fp);
  991.     else /* Simply seek by it */
  992.         fseek (fp, (long) (first_row * n2 * sizeof (unsigned char)), 1);
  993.     for (j = first_row, j2 = 0; j <= last_row; j++, j2++) {
  994.         ij = (j2 + pad[3]) * width_out + i_0_out;
  995.         fread ((char *) tmp, sizeof (unsigned char), n2, fp);
  996.         for (i = 0; i < width_in; i++, ij++) {
  997.             grid[ij] = tmp[first_col+i];
  998.             if (check && grid[ij] == grd_in_nan_value) grid[ij] = gmt_NaN;
  999.             if (bad_float ((double)grid[ij])) continue;
  1000.             header->z_min = MIN (header->z_min, grid[ij]);
  1001.             header->z_max = MAX (header->z_max, grid[ij]);
  1002.         }
  1003.     }
  1004.     if (piping)    /* Skip data by reading it */
  1005.         for (j = last_row + 1; j < header->ny; j++) fread ((char *) tmp, sizeof (unsigned char), n2, fp);
  1006.     
  1007.     header->nx = width_in;
  1008.     header->ny = height_in;
  1009.     header->x_min = w;
  1010.     header->x_max = e;
  1011.     header->y_min = s;
  1012.     header->y_max = n;
  1013.  
  1014.     if (fp != stdin) fclose (fp);
  1015.     
  1016.     free ((char *)tmp);
  1017.     return (FALSE);
  1018. }
  1019.  
  1020. int ras_write_grd (file, header, grid, w, e, s, n, pad, complex)
  1021. char *file;
  1022. struct GRD_HEADER *header;
  1023. float *grid;
  1024. double w, e, s, n;    /* Sub-region to write out */
  1025. int pad[];        /* padding on w, e, s, n of grid, respectively */
  1026. BOOLEAN complex;    /* TRUE if array is to hold real and imaginary parts */
  1027. {
  1028.     int i, i2, k, nm;
  1029.     int j, ij, j2, width_in, width_out, height_out, one_or_zero, n2;
  1030.     int first_col, last_col, first_row, last_row;
  1031.     
  1032.     BOOLEAN check;
  1033.     double small;
  1034.     
  1035.     unsigned char *tmp;
  1036.     
  1037.     FILE *fp;
  1038.     
  1039.     struct rasterfile h;
  1040.     
  1041.     if (complex) {
  1042.         fprintf (stderr, "GMT Fatal Error: Rasterfile %s cannot hold complex data!\n", file);
  1043.         exit (-1);
  1044.     }
  1045.     if (!strcmp (file, "=")) {
  1046.         fp = stdout;
  1047.     }
  1048.     else if ((fp = fopen (file, "w")) == NULL) {
  1049.         fprintf (stderr, "GMT Fatal Error: Could not create file %s!\n", file);
  1050.         exit (-1);
  1051.     }
  1052.     
  1053.     h.ras_magic = 0x59a66a95;
  1054.     h.ras_width = header->nx;
  1055.     h.ras_height = header->ny;
  1056.     h.ras_depth = 8;
  1057.     h.ras_length = header->ny * (int) ceil (header->nx/2.0) * 2;
  1058.     h.ras_type = 1;
  1059.     h.ras_maptype = 0;
  1060.     h.ras_maplength = 0;
  1061.     
  1062.     n2 = (int) ceil (header->nx / 2.0) * 2;
  1063.     tmp = (unsigned char *) memory (CNULL, n2, sizeof (unsigned char), "ras_write_grd");
  1064.     
  1065.     check = !bad_float (grd_out_nan_value);
  1066.  
  1067.     if (w == 0.0 && e == 0.0) {    /* Write entire file and return */
  1068.         nm = header->nx * header->ny;
  1069.         header->z_min = MAXDOUBLE;    header->z_max = -MAXDOUBLE;
  1070.         
  1071.         for (i = 0; i < nm; i++) {
  1072.                 ij = (complex) ? 2 * i : i;
  1073.                 if (bad_float ((double)grid[ij])) continue;
  1074.                 header->z_min = MIN (header->z_min, grid[ij]);
  1075.                 header->z_max = MAX (header->z_max, grid[ij]);
  1076.         }
  1077.         
  1078.         if (fwrite ((char *)&h, sizeof (struct rasterfile), 1, fp) != 1) {
  1079.             fprintf (stderr, "GMT Fatal Error: Error reading file %s!\n", file);
  1080.             exit (-1);
  1081.         }
  1082.  
  1083.         for (j = ij = 0; j < header->ny; j++) {
  1084.             for (i = 0; i < header->nx; i++, ij++) {
  1085.                 if (check && bad_float ((double)grid[ij])) grid[ij] = grd_out_nan_value;
  1086.                 tmp[i] = grid[ij];
  1087.             }
  1088.             if (fwrite ((char *)tmp, sizeof (unsigned char), n2, fp) != n2) {
  1089.                 fprintf (stderr, "GMT Fatal Error: Error writing file %s!\n", file);
  1090.                 exit (-1);
  1091.             }
  1092.         }
  1093.         return (FALSE);
  1094.     }
  1095.  
  1096.     one_or_zero = (header->node_offset) ? 0 : 1;
  1097.     
  1098.     /* Get dimension of subregion to write */
  1099.     
  1100.     width_out  = rint ((e - w) / header->x_inc) + one_or_zero;
  1101.     height_out = rint ((n - s) / header->y_inc) + one_or_zero;
  1102.     
  1103.     /* Get first and last row and column numbers */
  1104.     
  1105.     small = 0.1 * header->x_inc;
  1106.     first_col = floor ((w - header->x_min + small) / header->x_inc);
  1107.     last_col  = ceil ((e - header->x_min - small) / header->x_inc) -1 + one_or_zero;
  1108.     small = 0.1 * header->y_inc;
  1109.     first_row = floor ((header->y_max - n + small) / header->y_inc);
  1110.     last_row  = ceil ((header->y_max - s - small) / header->y_inc) -1 + one_or_zero;
  1111.     
  1112.     if ((last_col - first_col + 1) > width_out) last_col--;
  1113.     if ((last_row - first_row + 1) > height_out) last_row--;
  1114.     if ((last_col - first_col + 1) > width_out) first_col++;
  1115.     if ((last_row - first_row + 1) > height_out) first_row++;
  1116.  
  1117.     width_in = width_out;        /* Physical width of input array */
  1118.     if (pad[0] > 0) width_in += pad[0];
  1119.     if (pad[1] > 0) width_in += pad[1];
  1120.     
  1121.     header->x_min = w;
  1122.     header->x_max = e;
  1123.     header->y_min = s;
  1124.     header->y_max = n;
  1125.     
  1126.     h.ras_width = header->nx;
  1127.     h.ras_height = header->ny;
  1128.     h.ras_length = header->ny * (int) ceil (header->nx/2.0) * 2;
  1129.     
  1130.     /* store header information and array */
  1131.     
  1132.     if (fwrite ((char *)&h, sizeof (struct rasterfile), 1, fp) != 1) {
  1133.         fprintf (stderr, "GMT Fatal Error: Error writing file %s!\n", file);
  1134.         exit (-1);
  1135.     }
  1136.  
  1137.     i2 = first_col + pad[0];
  1138.     for (j = 0, j2 = first_row + pad[3]; j < height_out; j++, j2++) {
  1139.         ij = j2 * width_in + i2;
  1140.         for (i = 0; i < width_out; i++) {
  1141.             k = ij + i;
  1142.             if (check && bad_float ((double)grid[k])) grid[k] = grd_out_nan_value;
  1143.             tmp[i] = grid[k];
  1144.         }
  1145.         fwrite ((char *)tmp, sizeof (unsigned char), width_out, fp);
  1146.     }
  1147.     if (fp != stdout) fclose (fp);
  1148.     
  1149.     free ((char *)tmp);
  1150.     
  1151.     return (FALSE);
  1152.  
  1153. }
  1154. /* Add custom code here */
  1155.  
  1156. /*-----------------------------------------------------------
  1157.  * Format # :    4
  1158.  * Type :    
  1159.  * Prefix :    
  1160.  * Author :    
  1161.  * Date :    
  1162.  * Purpose:    
  1163.  * Functions :    
  1164.  *        
  1165.  *-----------------------------------------------------------*/
  1166.