/*-------------------------------------------------------------------- * The MB-system: mbgrdtiff.c 5/30/93 * $Id: mbgrdtiff.c 2230 2015-02-19 00:39:00Z caress $ * * Copyright (c) 1999-2015 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. *--------------------------------------------------------------------*/ /* * mbgrdtiff generates a TIFF image from a GMT grid. The * image generation is similar to that of the GMT program * grdimage. In particular, the color map is applied from * a GMT CPT file, and shading overlay grids may be applied. * The output TIFF file contains information allowing * the ArcView and ArcInfo GIS packages to import the image * as a geographically located coverage. The image is 8 bits * per pixel if the color map is a grayscale, and 24 bits * per pixel otherwise. */ /*--------------------------------------------------------------------*/ /* * The geotiff file variant produced by mbgrdtiff * has the structure given below. The image width in * pixels is denoted as nx and the image height as ny. * The images may be grayscale (1 byte per pixel) or * color (3 rgb bytes per pixel). * * Byte Size Value Meaning * ----------------------------------------------------------------- * * ------- Header -------- * * 0 2 'MM' Big-endian byte order used in this file * 2 2 42 Tiff identifier * 4 4 8 Byte offset of image file directory (IFD) * * ------- IFD ----------- * * 8 2 18 Number of entries in IFD * * 10 2 254 Tag: NewSubfileType * 12 2 4 Type: long (4 byte unsigned int) * 14 4 1 Count: one value * 18 4 0 Value: basic image, nothing fancy * * 22 2 256 Tag: ImageWidth * 24 2 4 Type: long (4 byte unsigned int) * 26 4 1 Count: one value * 30 4 nx Value: image width in pixels * * 34 2 257 Tag: ImageLength * 36 2 4 Type: long (4 byte unsigned int) * 38 4 1 Count: one value * 42 4 ny Value: image length (height) in pixels * * 46 2 258 Tag: BitsPerSample * 48 2 3 Type: short (2 byte unsigned int) * Color case: * 50 4 3 Count: one value for each of r,g,b (8,8,8) * 54 4 256 Value: offset to three BitsPerSample values * Grayscale case: * 50 4 1 Count: one value for grayscale * 54 4 8 Value: one byte per sample * * 58 2 259 Tag: Compression * 60 2 3 Type: short (2 byte unsigned int) * 62 4 1 Count: one value * 66 4 1 Value: no compression * * 70 2 262 Tag: PhotometricInterpretation * 72 2 3 Type: short (2 byte unsigned int) * 74 4 1 Count: one value * 78 4 2 Value: RGB image * * 82 2 273 Tag: StripOffsets * 84 2 4 Type: long (4 byte unsigned int) * 86 4 1 Count: entire image is one strip * 90 4 512 Value: offset to image data * * 94 2 277 Tag: SamplesPerPixel * 96 2 3 Type: short (2 byte unsigned int) * 98 4 1 Count: one value * 102 4 3 Value: 3 samples per pixel for RGB * * 106 2 278 Tag: RowsPerStrip * 108 2 4 Type: long (4 byte unsigned int) * 110 4 1 Count: one value * 114 4 ny Value: ImageLength value (all rows in image) * * 118 2 279 Tag: StripByteCounts * 120 2 4 Type: long (4 byte unsigned int) * 122 4 1 Count: one value * Color case: * 126 4 3*nx*ny Value: 3 * ImageWidth * ImageLength * Grayscale case: * 126 4 nx*ny Value: ImageWidth * ImageLength * * 130 2 282 Tag: XResolution * 132 2 5 Type: rational (2 ints: numerator, denominator) * 134 4 1 Count: one value * 138 4 264 Value: offset to fraction representing dpi * * 142 2 283 Tag: YResolution * 144 2 4 Type: rational (2 ints: numerator, denominator) * 146 4 1 Count: one value * 150 4 272 Value: offset to fraction representing dpi * * 154 2 296 Tag: ResolutionUnit * 156 2 3 Type: short (2 byte unsigned int) * 158 4 1 Count: one value * 162 4 2 Value: Inches * * 166 2 33550 Tag: ModelPixelScaleTag * 168 2 12 Type: double (IEEE double precision) * 170 4 3 Count: 3 for scalex,scaley,scalez where scalez=0 * 174 4 280 Value: offset to values * * 178 2 33922 Tag: ModelTiepointTag * 180 2 12 Type: double (IEEE double precision) * 182 4 6 Count: 6 for i,j,k,x,y,z where k=z=0 * 186 4 304 Value: offset to values * * 190 2 34735 Tag: GeoKeyDirectoryTag * 192 2 3 Type: short (2 byte unsigned int) * 194 4 20 Count: 20 * 198 4 352 Value: offset to values * * 190 2 34736 Tag: GeoDoubleParamsTag * 192 2 3 Type: double (IEEE double precision) * 194 4 ndouble Count: ndouble * 198 4 400 Value: offset to values * * 190 2 34737 Tag: GeoAsciiParamsTag * 192 2 2 Type: ASCII * 194 4 nascii Count: nascii * 198 4 448 Value: offset to values * * ------- Values -------- * * 256 2 8 BitsPerSample R (Color case only) * 258 2 8 BitsPerSample G (Color case only) * 260 2 8 BitsPerSample B (Color case only) * * 264 4 nx XResolution numerator (ImageWidth) * 268 4 4 XResolution denominator (4 inches) * * 272 4 nx YResolution numerator (ImageWidth) * 276 4 4 YResolution denominator (4 inches) * * 280 8 scalex ModelPixelScaleTag scalex * 288 8 scaley ModelPixelScaleTag scaley * 296 8 scalez ModelPixelScaleTag scalez * * ------- If GTModelTypeGeoKey = ModelTypeGeographic then: * 304 8 0 ModelTiePointTag i * 312 8 0 ModelTiePointTag j * 320 8 0 ModelTiePointTag k * 328 8 minlon ModelTiePointTag minimum longitude * 336 8 maxlat ModelTiePointTag minimum latitude * 344 8 0 ModelTiePointTag minimum z * * ------- Else if GTModelTypeGeoKey = ModelTypeProjected then: * 304 8 0 ModelTiePointTag i * 312 8 0 ModelTiePointTag j * 320 8 0 ModelTiePointTag k * 328 8 minX ModelTiePointTag minimum easting * 336 8 maxY ModelTiePointTag minimum northing * 344 8 0 ModelTiePointTag minimum z * * 352 2 1 GeoKeyDirectoryTag KeyDirectoryVersion * 354 2 0 GeoKeyDirectoryTag KeyRevision * 356 2 2 GeoKeyDirectoryTag MinorRevision * 358 2 4 GeoKeyDirectoryTag NumberOfKeys * * 360 2 1024 GeoKeyDirectoryTag KeyId: GTModelTypeGeoKey * 362 2 0 GeoKeyDirectoryTag TiffTagLocation * 364 2 1 GeoKeyDirectoryTag Count * 366 2 1 GeoKeyDirectoryTag ModelType Value (ModelTypeProjected=1, ModelTypeGeographic=2) * * 368 2 1025 GeoKeyDirectoryTag KeyId: GTRasterTypeGeoKey * 370 2 0 GeoKeyDirectoryTag TiffTagLocation * 372 2 1 GeoKeyDirectoryTag Count * 374 2 2 GeoKeyDirectoryTag RasterType Value (RasterPixelIsArea=1, RasterPixelIsPoint=2) * * 368 2 1026 GeoKeyDirectoryTag KeyId: GTCitationGeoKey * 370 2 34737 GeoKeyDirectoryTag TiffTagLocation * 372 2 nascii GeoKeyDirectoryTag Count * 374 2 0 GeoKeyDirectoryTag Value_Offset * * ------- If GTModelTypeGeoKey = ModelTypeGeographic then: * 376 2 2048 GeoKeyDirectoryTag KeyId: GeographicTypeGeoKey * 378 2 0 GeoKeyDirectoryTag TiffTagLocation * 380 2 1 GeoKeyDirectoryTag Count * 382 2 4030 GeoKeyDirectoryTag GeographicType Value (GCSE_WGS84=4030) * * ------- Else if GTModelTypeGeoKey = ModelTypeProjected then: * 376 2 3072 GeoKeyDirectoryTag KeyId: ProjectedCSTypeGeoKey * 378 2 0 GeoKeyDirectoryTag TiffTagLocation * 380 2 1 GeoKeyDirectoryTag Count * 382 2 32660 GeoKeyDirectoryTag ProjectedCSType Value (PCS_WGS84_UTM_zone_60N=32660) * * 400 8 ndouble GeoDoubleParamsTag * * 448 8 nascii GeoAsciiParamsTag * * ------- Image --------- * * 512 3*nx*ny Image in RGB bytes * * */ /*--------------------------------------------------------------------*/ /* MBIO include files */ #include "mb_status.h" #include "mb_format.h" #include "mb_define.h" /* TIFF 6.0 and geoTIFF tag array */ #define TIFF_HEADER_SIZE 1024 #define IMAGE_OFFSET TIFF_HEADER_SIZE #define TIFF_COMMENT_MAXLINE 64 #define NUMBER_TAGS 18 #define NewSubfileType 254 #define ImageWidth 256 #define ImageLength 257 #define BitsPerSample 258 #define Compression 259 #define PhotometricInterpretation 262 #define StripOffsets 273 #define SamplesPerPixel 277 #define RowsPerStrip 278 #define StripByteCounts 279 #define XResolution 282 #define YResolution 283 #define ResolutionUnit 296 #define ModelPixelScaleTag 33550 #define ModelTiepointTag 33922 #define GeoKeyDirectoryTag 34735 #define GeoDoubleParamsTag 34736 #define GeoAsciiParamsTag 34737 #define GTModelTypeGeoKey 1024 #define GTRasterTypeGeoKey 1025 #define GTCitationGeoKey 1026 #define GeographicTypeGeoKey 2048 #define ProjectedCSTypeGeoKey 3072 #define RasterPixelIsArea 1 #define RasterPixelIsPoint 2 #define ModelTypeProjected 1 #define ModelTypeGeographic 2 #define GCS_WGS_84 4326 unsigned short tiff_tag[] = { NewSubfileType, ImageWidth, ImageLength, BitsPerSample, Compression, PhotometricInterpretation, StripOffsets, SamplesPerPixel, RowsPerStrip, StripByteCounts, XResolution, YResolution, ResolutionUnit, ModelPixelScaleTag, ModelTiepointTag, GeoKeyDirectoryTag, GeoDoubleParamsTag, GeoAsciiParamsTag }; unsigned short tiff_type[] = { 4, /* NewSubfileType */ 4, /* ImageWidth */ 4, /* ImageLength */ 3, /* BitsPerSample */ 3, /* Compression */ 3, /* PhotometricInterpretation */ 4, /* StripOffsets */ 3, /* SamplesPerPixel */ 4, /* RowsPerStrip */ 4, /* StripByteCounts */ 5, /* XResolution */ 5, /* YResolution */ 3, /* ResolutionUnit */ 12, /* ModelPixelScaleTag */ 12, /* ModelTiepointTag */ 3, /* GeoKeyDirectoryTag */ 12, /* GeoDoubleParamsTag */ 2 /* GeoAsciiParamsTag */ }; int tiff_offset[] = { 0, /* NewSubfileType */ 0, /* ImageWidth */ 0, /* ImageLength */ 256, /* BitsPerSample */ 0, /* Compression */ 0, /* PhotometricInterpretation */ 0, /* StripOffsets */ 0, /* SamplesPerPixel */ 0, /* RowsPerStrip */ 0, /* StripByteCounts */ 264, /* XResolution */ 272, /* YResolution */ 0, /* ResolutionUnit */ 280, /* ModelPixelScaleTag */ 304, /* ModelTiepointTag */ 352, /* GeoKeyDirectoryTag */ 400, /* GeoDoubleParamsTag */ 448 /* GeoAsciiParamsTag */ }; /*--------------------------------------------------------------------*/ #define THIS_MODULE_NAME "mbgrdtiff" #define THIS_MODULE_LIB "core" #define THIS_MODULE_PURPOSE "Project grids or images and plot them on maps" #include "gmt_dev.h" //#define GMT_PROG_OPTIONS "->BJKOPRUVXYcfnptxy" GMT_OPT("S") #define GMT_PROG_OPTIONS "->JRVn" GMT_OPT("S") /* Control structure for mbgrdtiff */ struct MBGRDTIFF_CTRL { struct A { /* -A to write a GDAL file */ bool active; char *file; char *driver; } A; struct C { /* -C */ bool active; char *file; } C; struct D { /* -D to read GDAL file */ bool active; bool mode; /* Use info of -R option to reference image */ } D; struct E { /* -Ei| */ bool active; bool device_dpi; unsigned int dpi; } E; struct G { /* -G[f|b] */ bool active; double f_rgb[4]; double b_rgb[4]; } G; struct I { /* -I> */ bool active; bool do_rgb; char *file[3]; } I; struct Intensity { /* -I| */ bool active; bool constant; double value; char *file; } Intensity; struct M { /* -M */ bool active; } M; struct N { /* -N */ bool active; } N; struct mbO { /* -O */ bool active; char *file; } O; struct Q { /* -Q */ bool active; } Q; }; /*--------------------------------------------------------------------*/ void *New_mbgrdtiff_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */ struct MBGRDTIFF_CTRL *Ctrl; Ctrl = GMT_memory (GMT, NULL, 1, struct MBGRDTIFF_CTRL); /* Initialize values whose defaults are not 0/false/NULL */ Ctrl->G.b_rgb[0] = Ctrl->G.b_rgb[1] = Ctrl->G.b_rgb[2] = 1.0; return (Ctrl); } /*--------------------------------------------------------------------*/ void Free_mbgrdtiff_Ctrl (struct GMT_CTRL *GMT, struct MBGRDTIFF_CTRL *Ctrl) { /* Deallocate control structure */ int k; if (!Ctrl) return; for (k = 0; k < 3; k++) if (Ctrl->I.file[k]) free (Ctrl->I.file[k]); if (Ctrl->C.file) free (Ctrl->C.file); if (Ctrl->Intensity.file) free (Ctrl->Intensity.file); GMT_free (GMT, Ctrl); } /*--------------------------------------------------------------------*/ int GMT_mbgrdtiff_usage (struct GMTAPI_CTRL *API, int level) { GMT_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_NAME, THIS_MODULE_PURPOSE); if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR); GMT_Message (API, GMT_TIME_NONE, "usage: mbgrdtiff | %s [%s] [-C] [-Ei[|]]\n", GMT_J_OPT, GMT_B_OPT); GMT_Message (API, GMT_TIME_NONE, "\t[-G[f|b]] [-I|] [-K] [-M] [-N] [-O] [-P] [-Q]\n"); GMT_Message (API, GMT_TIME_NONE, "\t[%s] [-T] [%s] [%s]\n", GMT_Rgeo_OPT, GMT_U_OPT, GMT_V_OPT); GMT_Message (API, GMT_TIME_NONE, "\t[%s] [%s] [%s] [%s]\n\t[%s]\n\t[%s] [%s]\n\n", GMT_X_OPT, GMT_Y_OPT, GMT_c_OPT, GMT_f_OPT, GMT_n_OPT, GMT_p_OPT, GMT_t_OPT); if (level == GMT_SYNOPSIS) return (EXIT_FAILURE); GMT_Message (API, GMT_TIME_NONE, "\t is data set to be plotted. Its z-values are in user units and will be\n"); GMT_Message (API, GMT_TIME_NONE, "\t converted to rgb colors via the cpt file. Alternatively, give three separate\n"); GMT_Message (API, GMT_TIME_NONE, "\t grid files that contain the red, green, and blue components in the 0-255 range.\n"); GMT_Option (API, "J-"); GMT_Message (API, GMT_TIME_NONE, "\n\tOPTIONS:\n"); GMT_Option (API, "B-"); GMT_Message (API, GMT_TIME_NONE, "\t-C Color palette file to convert z to rgb. Optionally, instead give name of a master cpt\n"); GMT_Message (API, GMT_TIME_NONE, "\t to automatically assign 16 continuous colors over the data range [rainbow].\n"); GMT_Message (API, GMT_TIME_NONE, "\t-E Set dpi for the projected grid which must be constructed [100]\n"); GMT_Message (API, GMT_TIME_NONE, "\t if -Jx or -Jm is not selected [Default gives same size as input grid].\n"); GMT_Message (API, GMT_TIME_NONE, "\t Give i to do the interpolation in PostScript at device resolution.\n"); GMT_rgb_syntax (API->GMT, 'G', "Set transparency color for images that otherwise would result in 1-bit images.\n\t "); GMT_Message (API, GMT_TIME_NONE, "\t-I Use illumination. Append name of intensity grid file.\n"); GMT_Message (API, GMT_TIME_NONE, "\t For a constant intensity, just give the value instead.\n"); GMT_Option (API, "K"); GMT_Message (API, GMT_TIME_NONE, "\t-M Force monochrome image.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-N Do not clip image at the map boundary.\n"); GMT_Option (API, "O,P"); GMT_Message (API, GMT_TIME_NONE, "\t-Q Use PS Level 3 colormasking to make nodes with z = NaN transparent.\n"); GMT_Option (API, "R"); GMT_Option (API, "U,V,X,c,f,n,p,t,."); return (EXIT_FAILURE); } /*--------------------------------------------------------------------*/ int GMT_mbgrdtiff_parse (struct GMT_CTRL *GMT, struct MBGRDTIFF_CTRL *Ctrl, struct GMT_OPTION *options) { /* This parses the options provided to mbgrdtiff and sets parameters in Ctrl. * Note Ctrl has already been initialized and non-zero default values set. * Any GMT common options will override values set previously by other commands. * It also replaces any file names specified as input or output with the data ID * returned when registering these sources/destinations with the API. */ unsigned int n_errors = 0, n_files = 0; struct GMT_OPTION *opt = NULL; struct GMTAPI_CTRL *API = GMT->parent; for (opt = options; opt; opt = opt->next) { /* Process all the options given */ switch (opt->option) { case '<': /* Input file (only one or three is accepted) */ Ctrl->I.active = true; if (n_files >= 3) break; if (GMT_check_filearg (GMT, '<', opt->arg, GMT_IN)) Ctrl->I.file[n_files++] = strdup (opt->arg); else n_errors++; break; /* Processes program-specific parameters */ case 'C': /* CPT file */ Ctrl->C.active = true; if (Ctrl->C.file) free (Ctrl->C.file); Ctrl->C.file = strdup (opt->arg); break; case 'E': /* Sets dpi */ Ctrl->E.active = true; if (opt->arg[0] == 'i') /* Interpolate image to device resolution */ Ctrl->E.device_dpi = true; else if (opt->arg[0] == '\0') Ctrl->E.dpi = 100; /* Default grid dpi */ else Ctrl->E.dpi = atoi (opt->arg); break; case 'G': /* 1-bit fore or background color for transparent masks */ Ctrl->G.active = true; switch (opt->arg[0]) { case 'F': case 'f': if (GMT_getrgb (GMT, &opt->arg[1], Ctrl->G.f_rgb)) { GMT_rgb_syntax (GMT, 'G', " "); n_errors++; } else Ctrl->G.b_rgb[0] = -1; break; case 'B': case 'b': if (GMT_getrgb (GMT, &opt->arg[1], Ctrl->G.b_rgb)) { GMT_rgb_syntax (GMT, 'G', " "); n_errors++; } else Ctrl->G.f_rgb[0] = -1; break; default: /* Same as -Gf */ if (GMT_getrgb (GMT, opt->arg, Ctrl->G.f_rgb)) { GMT_rgb_syntax (GMT, 'G', " "); n_errors++; } else Ctrl->G.b_rgb[0] = -1; break; } break; case 'I': /* Input file (only one or three is accepted) */ /* if no grid file specified yet then first -Ifile sets * the primary grid file, and second -Ifile will be the * intensity file */ if (n_files == 0) { Ctrl->I.active = true; Ctrl->I.file[n_files++] = strdup (opt->arg); } else { Ctrl->Intensity.active = true; if (!GMT_access (GMT, opt->arg, R_OK)) /* Got a file */ Ctrl->Intensity.file = strdup (opt->arg); else if (opt->arg[0] && !GMT_not_numeric (GMT, opt->arg)) { /* Looks like a constant value */ Ctrl->Intensity.value = atof (opt->arg); Ctrl->Intensity.constant = true; } else { GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -I: Requires a valid grid file or a constant\n"); n_errors++; } } break; case 'M': /* Monochrome image */ Ctrl->M.active = true; break; case 'N': /* Do not clip at map boundary */ Ctrl->N.active = true; break; case 'O': /* Output file */ Ctrl->O.active = true; Ctrl->O.file = strdup (opt->arg); break; case 'Q': /* PS3 colormasking */ Ctrl->Q.active = true; break; default: /* Report bad options */ n_errors += GMT_default_error (GMT, opt->option); break; } } if (n_files == 3) Ctrl->I.do_rgb = true; // n_errors += GMT_check_condition (GMT, !GMT->common.J.active, // "Syntax error: Must specify a map projection with the -J option\n"); n_errors += GMT_check_condition (GMT, !Ctrl->C.file && !Ctrl->I.do_rgb, "Syntax error: Must specify color palette table\n"); n_errors += GMT_check_condition (GMT, !(n_files == 1 || n_files == 3), "Syntax error: Must specify one (or three) input file(s)\n"); n_errors += GMT_check_condition (GMT, Ctrl->Intensity.active && !Ctrl->Intensity.constant && !Ctrl->Intensity.file, "Syntax error -I option: Must specify intensity file or value\n"); n_errors += GMT_check_condition (GMT, Ctrl->E.active && !Ctrl->E.device_dpi && Ctrl->E.dpi <= 0, "Syntax error -E option: dpi must be positive\n"); n_errors += GMT_check_condition (GMT, Ctrl->G.f_rgb[0] < 0 && Ctrl->G.b_rgb[0] < 0, "Syntax error -G option: Only one of fore/back-ground can be transparent for 1-bit images\n"); n_errors += GMT_check_condition (GMT, Ctrl->M.active && Ctrl->Q.active, "Syntax error -Q option: Cannot use -M when doing colormasking\n"); return (n_errors ? GMT_PARSE_ERROR : GMT_OK); } /*--------------------------------------------------------------------*/ void GMT_mbgrdtiff_set_proj_limits (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *r, struct GMT_GRID_HEADER *g, bool projected) { /* Sets the projected extent of the grid given the map projection * The extreme x/y coordinates are returned in r, and dx/dy, and * nx/ny are set accordingly. Not that some of these may change * if GMT_project_init is called at a later stage */ unsigned int i, k; bool all_lats = false, all_lons = false; double x, y; r->nx = g->nx; r->ny = g->ny; r->registration = g->registration; r->n_bands = g->n_bands; /* By default, use entire plot region */ GMT_memcpy (r->wesn, GMT->current.proj.rect, 4, double); if (GMT->current.proj.projection == GMT_GENPER && GMT->current.proj.g_width != 0.0) return; if (GMT_is_geographic (GMT, GMT_IN)) { all_lats = GMT_180_RANGE (g->wesn[YHI], g->wesn[YLO]); all_lons = GMT_grd_is_global (GMT, g); if (all_lons && all_lats) return; /* Whole globe */ } /* Must search for extent along perimeter */ r->wesn[XLO] = r->wesn[YLO] = +DBL_MAX; r->wesn[XHI] = r->wesn[YHI] = -DBL_MAX; k = (g->registration == GMT_GRID_NODE_REG) ? 1 : 0; for (i = 0; i < g->nx - k; i++) { /* South and north sides */ GMT_geo_to_xy (GMT, g->wesn[XLO] + i * g->inc[GMT_X], g->wesn[YLO], &x, &y); r->wesn[XLO] = MIN (r->wesn[XLO], x), r->wesn[XHI] = MAX (r->wesn[XHI], x); r->wesn[YLO] = MIN (r->wesn[YLO], y), r->wesn[YHI] = MAX (r->wesn[YHI], y); GMT_geo_to_xy (GMT, g->wesn[XHI] - i * g->inc[GMT_X], g->wesn[YHI], &x, &y); r->wesn[XLO] = MIN (r->wesn[XLO], x), r->wesn[XHI] = MAX (r->wesn[XHI], x); r->wesn[YLO] = MIN (r->wesn[YLO], y), r->wesn[YHI] = MAX (r->wesn[YHI], y); } for (i = 0; i < g->ny - k; i++) { /* East and west sides */ GMT_geo_to_xy (GMT, g->wesn[XLO], g->wesn[YHI] - i * g->inc[GMT_Y], &x, &y); r->wesn[XLO] = MIN (r->wesn[XLO], x), r->wesn[XHI] = MAX (r->wesn[XHI], x); r->wesn[YLO] = MIN (r->wesn[YLO], y), r->wesn[YHI] = MAX (r->wesn[YHI], y); GMT_geo_to_xy (GMT, g->wesn[XHI], g->wesn[YLO] + i * g->inc[GMT_Y], &x, &y); r->wesn[XLO] = MIN (r->wesn[XLO], x), r->wesn[XHI] = MAX (r->wesn[XHI], x); r->wesn[YLO] = MIN (r->wesn[YLO], y), r->wesn[YHI] = MAX (r->wesn[YHI], y); } if (projected) { if (all_lons) { /* Full 360, use min/max for x */ r->wesn[XLO] = GMT->current.proj.rect[XLO]; r->wesn[XHI] = GMT->current.proj.rect[XHI]; } if (all_lats) { /* Full -90/+90, use min/max for y */ r->wesn[YLO] = GMT->current.proj.rect[YLO]; r->wesn[YHI] = GMT->current.proj.rect[YHI]; } } } /*--------------------------------------------------------------------*/ #define bailout(code) {GMT_Free_Options (mode); return (code);} #define Return(code) {Free_mbgrdtiff_Ctrl (GMT, Ctrl); GMT_end_module (GMT, GMT_cpy); bailout (code);} int GMT_mbgrdtiff (void *V_API, int mode, void *args) { char program_name[] = "mbswath"; /* TIFF arrays */ mb_path world_file; size_t image_size = 0; int modeltype; int projectionid; mb_path projectionname; char tiff_header[TIFF_HEADER_SIZE]; void *tiff_image = NULL; char tiff_comment[TIFF_COMMENT_MAXLINE]; char NorS; FILE *tfp; int nscan; int utmzone; int keyindex; char *projection = "-Jx1.0"; short value_short; int value_int; double value_double; size_t write_size; int i; bool done, need_to_project, normal_x, normal_y, resampled = false, gray_only = false; bool nothing_inside = false, use_intensity_grid; unsigned int k, nx = 0, ny = 0, grid_registration = GMT_GRID_NODE_REG, n_grids; unsigned int colormask_offset = 0, try, row, actual_row, col; uint64_t node_RGBA = 0; /* uint64_t for the RGB(A) image array. */ uint64_t node, kk, nm, byte; int index = 0, ks, error = 0; int nx8, shift, b_or_w, nx_pixels, k8; unsigned char *bitimage_1 = NULL, *bitimage_8 = NULL, *bitimage_24 = NULL, *rgb_used = NULL, i_rgb[3]; double dx, dy, x_side, y_side, x0 = 0.0, y0 = 0.0, rgb[4] = {0.0, 0.0, 0.0, 0.0}; double *NaN_rgb = NULL, red[4] = {1.0, 0.0, 0.0, 0.0}, wesn[4]; struct GMT_GRID *Grid_orig[3] = {NULL, NULL, NULL}, *Grid_proj[3] = {NULL, NULL, NULL}; struct GMT_GRID *Intens_orig = NULL, *Intens_proj = NULL; struct GMT_PALETTE *P = NULL; struct MBGRDTIFF_CTRL *Ctrl = NULL; struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL; /* General GMT interal parameters */ struct GMT_OPTION *options = NULL; struct PSL_CTRL *PSL = NULL; /* General PSL interal parameters */ struct GMT_GRID_HEADER *header_work = NULL; /* Pointer to a GMT header for the image or grid */ struct GMTAPI_CTRL *API = GMT_get_API_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */ /*----------------------- Standard module initialization and parsing ----------------------*/ if (API == NULL) return (GMT_NOT_A_SESSION); if (mode == GMT_MODULE_PURPOSE) return (GMT_mbgrdtiff_usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */ options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */ if (!options || options->option == GMT_OPT_USAGE) bailout (GMT_mbgrdtiff_usage (API, GMT_USAGE)); /* Return the usage message */ if (options->option == GMT_OPT_SYNOPSIS) bailout (GMT_mbgrdtiff_usage (API, GMT_SYNOPSIS)); /* Return the synopsis */ /* Parse the command-line arguments */ GMT = GMT_begin_module (API, THIS_MODULE_LIB, THIS_MODULE_NAME, &GMT_cpy); /* Save current state */ if (GMT_Parse_Common (API, GMT_PROG_OPTIONS, options)) Return (API->error); Ctrl = New_mbgrdtiff_Ctrl (GMT); /* Allocate and initialize a new control structure */ if ((error = GMT_mbgrdtiff_parse (GMT, Ctrl, options))) Return (error); /*---------------------------- This is the mbgrdtiff main code ----------------------------*/ n_grids = (Ctrl->I.do_rgb) ? 3 : 1; use_intensity_grid = (Ctrl->Intensity.active && !Ctrl->Intensity.constant); /* We want to use the intensity grid */ /* Read the illumination grid header right away so we can use its region to set that of an image (if requested) */ if (use_intensity_grid) { /* Illumination grid wanted */ GMT_Report (API, GMT_MSG_VERBOSE, "Allocates memory and read intensity file\n"); if ((Intens_orig = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_HEADER_ONLY, NULL, Ctrl->Intensity.file, NULL)) == NULL) { /* Get header only */ Return (API->error); } } GMT_Report (API, GMT_MSG_VERBOSE, "Allocates memory and read data file\n"); if (!Ctrl->D.active) { for (k = 0; k < n_grids; k++) { if ((Grid_orig[k] = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_HEADER_ONLY, NULL, Ctrl->I.file[k], NULL)) == NULL) { /* Get header only */ Return (API->error); } } if (!Ctrl->C.active) Ctrl->C.active = true; /* Use default CPT stuff */ } if (n_grids) header_work = Grid_orig[0]->header; /* OK, we are in GRID mode and this was not set further above. Do it now */ if (n_grids && Ctrl->I.do_rgb) { /* Must ensure all three grids are coregistered */ if (!GMT_grd_same_region (GMT, Grid_orig[0], Grid_orig[1])) error++; if (!GMT_grd_same_region (GMT, Grid_orig[0], Grid_orig[2])) error++; if (!(Grid_orig[0]->header->inc[GMT_X] == Grid_orig[1]->header->inc[GMT_X] && Grid_orig[0]->header->inc[GMT_X] == Grid_orig[2]->header->inc[GMT_X])) error++; if (!(Grid_orig[0]->header->nx == Grid_orig[1]->header->nx && Grid_orig[0]->header->nx == Grid_orig[2]->header->nx)) error++; if (!(Grid_orig[0]->header->ny == Grid_orig[1]->header->ny && Grid_orig[0]->header->ny == Grid_orig[2]->header->ny)) error++; if (!(Grid_orig[0]->header->registration == Grid_orig[1]->header->registration && Grid_orig[0]->header->registration == Grid_orig[2]->header->registration)) error++; if (error) { GMT_Report (API, GMT_MSG_NORMAL, "The r, g, and b grids are not congruent\n"); Return (EXIT_FAILURE); } } /* Determine what wesn to pass to map_setup */ if (!GMT->common.R.active && n_grids) GMT_memcpy (GMT->common.R.wesn, Grid_orig[0]->header->wesn, 4, double); GMT_err_fail (GMT, GMT_map_setup (GMT, GMT->common.R.wesn), ""); /* Determine if grid is to be projected */ need_to_project = (GMT_IS_NONLINEAR_GRATICULE (GMT) || Ctrl->E.dpi > 0); if (need_to_project) GMT_Report (API, GMT_MSG_DEBUG, "Projected grid is non-orthogonal, nonlinear, or dpi was changed\n"); /* Determine the wesn to be used to read the grid file; or bail if file is outside -R */ if (!GMT_grd_setregion (GMT, header_work, wesn, need_to_project * GMT->common.n.interpolant)) nothing_inside = true; else if (use_intensity_grid && !GMT_grd_setregion (GMT, Intens_orig->header, wesn, need_to_project * GMT->common.n.interpolant)) nothing_inside = true; if (nothing_inside) { /* No grid to plot; just do empty map and bail */ PSL = GMT_plotinit (GMT, options); GMT_plane_perspective (GMT, GMT->current.proj.z_project.view_plane, GMT->current.proj.z_level); GMT_plotcanvas (GMT); /* Fill canvas if requested */ GMT_map_basemap (GMT); GMT_plane_perspective (GMT, -1, 0.0); GMT_plotend (GMT); Return (EXIT_SUCCESS); } if (n_grids) { nx = GMT_get_n (GMT, wesn[XLO], wesn[XHI], Grid_orig[0]->header->inc[GMT_X], Grid_orig[0]->header->registration); ny = GMT_get_n (GMT, wesn[YLO], wesn[YHI], Grid_orig[0]->header->inc[GMT_Y], Grid_orig[0]->header->registration); } // if (!Ctrl->A.active) { /* Otherwise we are not writting any postscript */ // PSL = GMT_plotinit (GMT, options); // GMT_plane_perspective (GMT, GMT->current.proj.z_project.view_plane, GMT->current.proj.z_level); // GMT_plotcanvas (GMT); /* Fill canvas if requested */ // if (!Ctrl->N.active) GMT_map_clip_on (GMT, GMT->session.no_rgb, 3); // } /* Read data */ for (k = 0; k < n_grids; k++) { if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_DATA_ONLY, wesn, Ctrl->I.file[k], Grid_orig[k]) == NULL) { /* Get grid data */ Return (API->error); } } /* If given, get intensity file or compute intensities */ if (use_intensity_grid) { /* Illumination wanted */ GMT_Report (API, GMT_MSG_VERBOSE, "Allocates memory and read intensity file\n"); /* Remember, the illumination header was already read at the top */ if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_DATA_ONLY, wesn, Ctrl->Intensity.file, Intens_orig) == NULL) { Return (API->error); /* Get grid data */ } if (n_grids && (Intens_orig->header->nx != Grid_orig[0]->header->nx || Intens_orig->header->ny != Grid_orig[0]->header->ny)) { GMT_Report (API, GMT_MSG_NORMAL, "Intensity file has improper dimensions!\n"); Return (EXIT_FAILURE); } } if (need_to_project) { /* Need to resample the grd file */ int nx_proj = 0, ny_proj = 0; double inc[2] = {0.0, 0.0}; GMT_Report (API, GMT_MSG_VERBOSE, "project grid files\n"); if (Ctrl->E.dpi == 0) { /* Use input # of nodes as # of projected nodes */ nx_proj = nx; ny_proj = ny; } for (k = 0; k < n_grids; k++) { if (!Grid_proj[k] && (Grid_proj[k] = GMT_Duplicate_Data (API, GMT_IS_GRID, GMT_DUPLICATE_NONE, Grid_orig[k])) == NULL) Return (API->error); /* Just to get a header we can change */ GMT_mbgrdtiff_set_proj_limits (GMT, Grid_proj[k]->header, Grid_orig[k]->header, need_to_project); if (grid_registration == GMT_GRID_NODE_REG) /* Force pixel if dpi is set */ grid_registration = (Ctrl->E.dpi > 0) ? GMT_GRID_PIXEL_REG : Grid_orig[k]->header->registration; GMT_err_fail (GMT, GMT_project_init (GMT, Grid_proj[k]->header, inc, nx_proj, ny_proj, Ctrl->E.dpi, grid_registration), Ctrl->I.file[k]); GMT_set_grddim (GMT, Grid_proj[k]->header); if (GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_GRID_DATA_ONLY, NULL, NULL, NULL, 0, 0, Grid_proj[k]) == NULL) Return (API->error); GMT_grd_project (GMT, Grid_orig[k], Grid_proj[k], false); if (GMT_Destroy_Data (API, &Grid_orig[k]) != GMT_OK) { Return (API->error); } } if (use_intensity_grid) { if ((Intens_proj = GMT_Duplicate_Data (API, GMT_IS_GRID, GMT_DUPLICATE_NONE, Intens_orig)) == NULL) Return (API->error); /* Just to get a header we can change */ if (n_grids) GMT_memcpy (Intens_proj->header->wesn, Grid_proj[0]->header->wesn, 4, double); if (Ctrl->E.dpi == 0) { /* Use input # of nodes as # of projected nodes */ nx_proj = Intens_orig->header->nx; ny_proj = Intens_orig->header->ny; } GMT_err_fail (GMT, GMT_project_init (GMT, Intens_proj->header, inc, nx_proj, ny_proj, Ctrl->E.dpi, grid_registration), Ctrl->Intensity.file); GMT_set_grddim (GMT, Intens_proj->header); if (GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_GRID_DATA_ONLY, NULL, NULL, NULL, 0, 0, Intens_proj) == NULL) Return (API->error); GMT_grd_project (GMT, Intens_orig, Intens_proj, false); if (GMT_Destroy_Data (API, &Intens_orig) != GMT_OK) { Return (API->error); } } resampled = true; } else { /* Simply set Grid_proj[i]/Intens_proj to point to Grid_orig[i]/Intens_orig */ struct GMT_GRID_HEADER tmp_header; for (k = 0; k < n_grids; k++) { /* Must get a copy so we can change one without affecting the other */ GMT_memcpy (&tmp_header, Grid_orig[k]->header, 1, struct GMT_GRID_HEADER); Grid_proj[k] = Grid_orig[k]; //GMT_mbgrdtiff_set_proj_limits (GMT, Grid_proj[k]->header, &tmp_header, need_to_project); } if (use_intensity_grid) Intens_proj = Intens_orig; if (n_grids) grid_registration = Grid_orig[0]->header->registration; } if (n_grids) { Grid_proj[0]->header->n_bands = 1; header_work = Grid_proj[0]->header; /* Later when need to refer to the header, use this copy */ } nm = header_work->nm; nx = header_work->nx; ny = header_work->ny; /* Get/calculate a color palette file */ if (!Ctrl->I.do_rgb) { if (Ctrl->C.active) { /* Read palette file */ if ((P = GMT_Get_CPT (GMT, Ctrl->C.file, GMT_CPT_OPTIONAL, header_work->z_min, header_work->z_max)) == NULL) { Return (API->error); } gray_only = (P && P->is_gray); } } if (P && P->has_pattern) GMT_Report (API, GMT_MSG_VERBOSE, "Warning: Patterns in cpt file only apply to -T\n"); GMT_Report (API, GMT_MSG_VERBOSE, "Evaluate pixel colors\n"); NaN_rgb = (P) ? P->patch[GMT_NAN].rgb : GMT->current.setting.color_patch[GMT_NAN]; if (Ctrl->Q.active) { if (gray_only) { GMT_Report (API, GMT_MSG_VERBOSE, "Your image is grayscale only but -Q requires 24-bit; image will be converted to 24-bit.\n"); gray_only = false; NaN_rgb = red; /* Arbitrarily pick red as the NaN color since image is gray only */ GMT_memcpy (P->patch[GMT_NAN].rgb, red, 4, double); } rgb_used = GMT_memory (GMT, NULL, 256*256*256, unsigned char); } if (Ctrl->M.active || gray_only) { image_size = nm; bitimage_8 = GMT_memory (GMT, NULL, image_size, unsigned char); tiff_image = bitimage_8; } else { if (Ctrl->Q.active) colormask_offset = 3; image_size = 3 * nm + colormask_offset; bitimage_24 = GMT_memory (GMT, NULL, image_size, unsigned char); if (P && Ctrl->Q.active) { for (k = 0; k < 3; k++) bitimage_24[k] = GMT_u255 (P->patch[GMT_NAN].rgb[k]); } tiff_image = bitimage_24; } normal_x = !(GMT->current.proj.projection == GMT_LINEAR && !GMT->current.proj.xyz_pos[0] && !resampled); normal_y = !(GMT->current.proj.projection == GMT_LINEAR && !GMT->current.proj.xyz_pos[1] && !resampled); normal_x = true; normal_y = true; for (try = 0, done = false; !done && try < 2; try++) { /* Evaluate colors at least once, or twice if -Q and we need to select another NaN color */ for (row = 0, byte = colormask_offset; row < ny; row++) { actual_row = (normal_y) ? row : ny - row - 1; kk = GMT_IJPGI (header_work, actual_row, 0); if (Ctrl->D.active && row == 0) node_RGBA = kk; /* First time per row equals 'node', after grows alone */ for (col = 0; col < nx; col++) { /* Compute rgb for each pixel */ node = kk + (normal_x ? col : nx - col - 1); if (Ctrl->I.do_rgb) { for (k = 0; k < 3; k++) { if (GMT_is_fnan (Grid_proj[k]->data[node])) { /* If one is NaN they are all assumed to be NaN */ k = 3; /* To exit the k-loop */ GMT_rgb_copy (rgb, NaN_rgb); index = GMT_NAN - 3; /* Ensures no illumination done later */ } else { /* Set color, let index = 0 so illuminate test will work */ rgb[k] = GMT_is255 (Grid_proj[k]->data[node]); if (rgb[k] < 0.0) rgb[k] = 0.0; else if (rgb[k] > 1.0) rgb[k] = 1.0; /* Clip */ index = 0; } } } else index = GMT_get_rgb_from_z (GMT, P, Grid_proj[0]->data[node], rgb); if (Ctrl->I.active && index != GMT_NAN - 3) { if (!n_grids) { /* Here we are illuminating an image. Must recompute "node" with the GMT_IJP macro */ node = GMT_IJP (Intens_proj->header, actual_row, 0) + (normal_x ? col : nx - col - 1); } if (use_intensity_grid) GMT_illuminate (GMT, Intens_proj->data[node], rgb); else GMT_illuminate (GMT, Ctrl->Intensity.value, rgb); } if (P && gray_only) /* Color table only has grays, pick r */ bitimage_8[byte++] = GMT_u255 (rgb[0]); else if (Ctrl->M.active) /* Convert rgb to gray using the GMT_YIQ transformation */ bitimage_8[byte++] = GMT_u255 (GMT_YIQ (rgb)); else { for (k = 0; k < 3; k++) bitimage_24[byte++] = i_rgb[k] = GMT_u255 (rgb[k]); if (Ctrl->Q.active && index != GMT_NAN - 3) /* Keep track of all r/g/b combinations used except for NaN */ rgb_used[(i_rgb[0]*256 + i_rgb[1])*256+i_rgb[2]] = true; } } if (!n_grids) node_RGBA += header_work->n_bands * (header_work->pad[XLO] + header_work->pad[XHI]); } if (P && Ctrl->Q.active) { /* Check that we found an unused r/g/b value so colormasking will work OK */ index = (GMT_u255(P->patch[GMT_NAN].rgb[0])*256 + GMT_u255(P->patch[GMT_NAN].rgb[1]))*256 + GMT_u255(P->patch[GMT_NAN].rgb[2]); if (rgb_used[index]) { /* This r/g/b already appears in the image as a non-NaN color; we must find a replacement NaN color */ for (index = 0, ks = -1; ks == -1 && index < 256*256*256; index++) if (!rgb_used[index]) ks = index; if (ks == -1) { GMT_Report (API, GMT_MSG_NORMAL, "Warning: Colormasking will fail as there is no unused color that can represent transparency\n"); done = true; } else { /* Pick the first unused color (i.e., k) and let it play the role of the NaN color for transparency */ bitimage_24[0] = (unsigned char)(ks >> 16); bitimage_24[1] = (unsigned char)((ks >> 8) & 255); bitimage_24[2] = (unsigned char)(ks & 255); GMT_Report (API, GMT_MSG_VERBOSE, "Warning: transparency color reset from %s to color %d/%d/%d\n", GMT_putrgb (GMT, P->patch[GMT_NAN].rgb), (int)bitimage_24[0], (int)bitimage_24[1], (int)bitimage_24[2]); for (k = 0; k < 3; k++) P->patch[GMT_NAN].rgb[k] = GMT_is255 (bitimage_24[k]); /* Set new NaN color */ } } } else done = true; } if (Ctrl->Q.active) GMT_free (GMT, rgb_used); for (k = 1; k < n_grids; k++) { /* Not done with Grid_proj[0] yet, hence we start loop at k = 1 */ if (need_to_project && GMT_Destroy_Data (API, &Grid_proj[k]) != GMT_OK) { GMT_Report (API, GMT_MSG_NORMAL, "Failed to free Grid_proj[k]\n"); } } if (use_intensity_grid) { if (need_to_project || !n_grids) { if (GMT_Destroy_Data (API, &Intens_proj) != GMT_OK) { GMT_Report (API, GMT_MSG_NORMAL, "Failed to free Intens_proj\n"); } } } /* Get actual size of each pixel */ dx = GMT_get_inc (GMT, header_work->wesn[XLO], header_work->wesn[XHI], header_work->nx, header_work->registration); dy = GMT_get_inc (GMT, header_work->wesn[YLO], header_work->wesn[YHI], header_work->ny, header_work->registration); /* Set lower left position of image on map */ x0 = header_work->wesn[XLO]; y0 = header_work->wesn[YLO]; if (grid_registration == GMT_GRID_NODE_REG) { /* Grid registration, move 1/2 pixel down/left */ x0 -= 0.5 * dx; y0 -= 0.5 * dy; } x_side = dx * header_work->nx; y_side = dy * header_work->ny; if (P && gray_only) for (kk = 0, P->is_bw = true; P->is_bw && kk < nm; kk++) if (!(bitimage_8[kk] == 0 || bitimage_8[kk] == 255)) P->is_bw = false; if (P && P->is_bw) { /* Can get away with 1 bit image */ GMT_Report (API, GMT_MSG_VERBOSE, "Creating 1-bit B/W image\n"); nx8 = irint (ceil (nx / 8.0)); /* Image width must equal a multiple of 8 bits */ nx_pixels = nx8 * 8; image_size = nx8 * ny; bitimage_1 = GMT_memory (GMT, NULL, image_size, unsigned char); tiff_image = bitimage_1; for (row = k = k8 = 0; row < ny; row++) { shift = 0; byte = 0; for (col = 0; col < nx; col++, k++) { b_or_w = (bitimage_8[k] == 255); byte |= b_or_w; shift++; if (shift == 8) { /* Time to dump out byte */ bitimage_1[k8++] = (unsigned char) byte; byte = shift = 0; } else byte <<= 1; /* Move the bits we have so far 1 step to the left */ } if (shift) { /* Set the remaining bits in this bit to white */ byte |= 1; shift++; while (shift < 8) { byte <<= 1; byte |= 1; shift++; } bitimage_1[k8++] = (unsigned char) byte; } } x_side = nx_pixels * dx; //PSL_plotbitimage (PSL, x0, y0, x_side, y_side, PSL_BL, bit, nx_pixels, ny, Ctrl->G.f_rgb, Ctrl->G.b_rgb); } else if ((P && gray_only) || Ctrl->M.active) { GMT_Report (API, GMT_MSG_VERBOSE, "Creating 8-bit grayshade image\n"); //PSL_plotcolorimage (PSL, x0, y0, x_side, y_side, PSL_BL, bitimage_8, nx, ny, (Ctrl->E.device_dpi ? -8 : 8)); } else { GMT_Report (API, GMT_MSG_VERBOSE, "Creating 24-bit color image\n"); //PSL_plotcolorimage (PSL, x0, y0, x_side, y_side, PSL_BL, bitimage_24, (Ctrl->Q.active ? -1 : 1) * // nx, ny, (Ctrl->E.device_dpi ? -24 : 24)); } /*------------------------- Write out the GeoTiff and world files -------------------------*/ /* try to get projection from the grd file remark */ if (strncmp(&(header_work->remark[2]), "Projection: ", 12) == 0) { if ((nscan = sscanf(&(header_work->remark[2]), "Projection: UTM%d%c", &utmzone, &NorS)) == 2) { if (NorS == 'N') { projectionid = 32600 + utmzone; } else if (NorS == 'S') { projectionid = 32700 + utmzone; } modeltype = ModelTypeProjected; sprintf(projectionname, "UTM%2.2d%c", utmzone, NorS); } else if ((nscan = sscanf(&(header_work->remark[2]), "Projection: epsg%d", &projectionid)) == 1) { sprintf(projectionname, "epsg%d", projectionid); modeltype = ModelTypeProjected; } else if (strncmp(&(header_work->remark[2]), "Projection: SeismicProfile", 26) == 0) { sprintf(projectionname, "SeismicProfile"); modeltype = ModelTypeProjected; } else { strcpy(projectionname, "Geographic WGS84"); modeltype = ModelTypeGeographic; projectionid = GCS_WGS_84; } } else { strcpy(projectionname, "Geographic WGS84"); modeltype = ModelTypeGeographic; projectionid = GCS_WGS_84; } /* set the TIFF comment */ sprintf(tiff_comment, "Image generated by %s|", program_name); /* set the TIFF header */ memset(tiff_header,0,TIFF_HEADER_SIZE); index = 0; tiff_header[0] = 'M'; tiff_header[1] = 'M'; index += 2; value_short = 42; mb_put_binary_short(MB_NO, value_short, &tiff_header[index]); index += 2; value_int = 8; mb_put_binary_int(MB_NO, value_int, &tiff_header[index]); index += 4; /* number of entries in IFD */ value_short = NUMBER_TAGS; mb_put_binary_short(MB_NO, value_short, &tiff_header[index]); index += 2; /* loop over all tags */ for (i=0;iinc[0]; mb_put_binary_double(MB_NO, value_double, &tiff_header[tiff_offset[i]]); value_double = header_work->inc[1]; mb_put_binary_double(MB_NO, value_double, &tiff_header[tiff_offset[i]+8]); value_double = 0.0; mb_put_binary_double(MB_NO, value_double, &tiff_header[tiff_offset[i]+16]); break; case ModelTiepointTag: value_int = 6; mb_put_binary_int(MB_NO, value_int, &tiff_header[index]); index += 4; value_int = tiff_offset[i]; mb_put_binary_int(MB_NO, value_int, &tiff_header[index]); index += 4; value_double = 0; mb_put_binary_double(MB_NO, value_double, &tiff_header[tiff_offset[i]]); value_double = 0; mb_put_binary_double(MB_NO, value_double, &tiff_header[tiff_offset[i] + 1 * 8]); value_double = 0; mb_put_binary_double(MB_NO, value_double, &tiff_header[tiff_offset[i] + 2 * 8]); value_double = header_work->wesn[0] - 0.5 * header_work->inc[0]; mb_put_binary_double(MB_NO, value_double, &tiff_header[tiff_offset[i] + 3 * 8]); value_double = header_work->wesn[3] + 0.5 * header_work->inc[1]; mb_put_binary_double(MB_NO, value_double, &tiff_header[tiff_offset[i] + 4 * 8]); value_double = 0.0; mb_put_binary_double(MB_NO, value_double, &tiff_header[tiff_offset[i] + 5 * 8]); break; case GeoKeyDirectoryTag: value_int = 20; mb_put_binary_int(MB_NO, value_int, &tiff_header[index]); index += 4; value_int = tiff_offset[i]; mb_put_binary_int(MB_NO, value_int, &tiff_header[index]); index += 4; /* index to geotiff geokey directory */ keyindex = tiff_offset[i]; /* geokey directory header (KeyDirectoryVersion, KeyRevision, MinorRevision, NumberOfKeys) */ value_short = 1; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = 0; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = 2; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = 4; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; /* GTModelTypeGeoKey */ value_short = GTModelTypeGeoKey; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = 0; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = 1; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = modeltype; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; /* GTRasterTypeGeoKey */ value_short = GTRasterTypeGeoKey; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = 0; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = 1; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = RasterPixelIsPoint; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; /* GTCitationGeoKey */ value_short = GTCitationGeoKey; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = GeoAsciiParamsTag; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = strlen(tiff_comment); mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = 0; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; if (modeltype == ModelTypeGeographic) { /* GeographicTypeGeoKey */ value_short = GeographicTypeGeoKey; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = 0; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = 1; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = projectionid; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; } else if (modeltype == ModelTypeProjected) { /* ProjectedCSTypeGeoKey */ value_short = ProjectedCSTypeGeoKey; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = 0; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = 1; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; value_short = projectionid; mb_put_binary_short(MB_NO, value_short, &tiff_header[keyindex]); keyindex += 2; } break; case GeoDoubleParamsTag: value_int = 1; mb_put_binary_int(MB_NO, value_int, &tiff_header[index]); index += 4; value_int = tiff_offset[i]; mb_put_binary_int(MB_NO, value_int, &tiff_header[index]); index += 4; break; case GeoAsciiParamsTag: value_int = strlen(tiff_comment); mb_put_binary_int(MB_NO, value_int, &tiff_header[index]); index += 4; value_int = tiff_offset[i]; mb_put_binary_int(MB_NO, value_int, &tiff_header[index]); index += 4; /* put in the string */ strncpy(&tiff_header[tiff_offset[i]], tiff_comment, 64); break; } } /* open TIFF file */ if ((tfp = fopen(Ctrl->O.file,"w")) == NULL) { API->error++; return(API->error); } /* set the TIFF comment */ sprintf(tiff_comment, "Image generated by %s|", program_name); /* write the header */ if ((write_size = fwrite(tiff_header,1,TIFF_HEADER_SIZE,tfp)) != TIFF_HEADER_SIZE) { API->error++; return(API->error); } /* write the image */ if ((write_size = fwrite(tiff_image,1,image_size,tfp)) != image_size) { API->error++; return(API->error); } /* close the tiff file */ fclose(tfp); /* open world file */ strcpy(world_file, Ctrl->O.file); world_file[strlen(Ctrl->O.file)-4] = '\0'; strcat(world_file,".tfw"); if ((tfp = fopen(world_file,"w")) == NULL) { API->error++; return(API->error); } /* write out world file contents */ fprintf(tfp, "%f\r\n0.0\r\n0.0\r\n%f\r\n%f\r\n%f\r\n", dx, -dy, header_work->wesn[XLO] - 0.5 * dx, header_work->wesn[YHI] + 0.5 * dy); /* close the world file */ fclose(tfp); fprintf(stderr,"3 Grid header:\n\tnx:%d ny:%d registration:%d\n\tWESN: %f %f %f %f\n\tinc: %f %f\n", Grid_orig[0]->header->nx,Grid_orig[0]->header->ny,Grid_orig[0]->header->registration, Grid_orig[0]->header->wesn[XLO],Grid_orig[0]->header->wesn[XHI],Grid_orig[0]->header->wesn[YLO],Grid_orig[0]->header->wesn[YHI], Grid_orig[0]->header->inc[0],Grid_orig[0]->header->inc[1]); fprintf(stderr,"3 Work header:\n\tnx:%d ny:%d registration:%d\n\tWESN: %f %f %f %f\n\tinc: %f %f\n", header_work->nx,header_work->ny,header_work->registration, header_work->wesn[XLO],header_work->wesn[XHI],header_work->wesn[YLO],header_work->wesn[YHI], header_work->inc[0],header_work->inc[1]); // if (!Ctrl->A.active) { // if (!Ctrl->N.active) GMT_map_clip_off (GMT); // // GMT_map_basemap (GMT); // GMT_plane_perspective (GMT, -1, 0.0); // GMT_plotend (GMT); // } /* Free bitimage arrays. GMT_free will not complain if they have not been used (NULL) */ if (P && P->is_bw) GMT_free (GMT, bitimage_1); if (bitimage_8) GMT_free (GMT, bitimage_8); if (bitimage_24) GMT_free (GMT, bitimage_24); if (need_to_project && n_grids && GMT_Destroy_Data (API, &Grid_proj[0]) != GMT_OK) { GMT_Report (API, GMT_MSG_NORMAL, "Failed to free Grid_proj[0]\n"); } if (!Ctrl->C.active && GMT_Destroy_Data (API, &P) != GMT_OK) { Return (API->error); } Return (EXIT_SUCCESS); }