/********************************************************************** * $Id: geotiff2pgraster.c 2007-07-22 XingLin $ * * PostGIS - Spatial Types for PostgreSQL * http://postgis.refractions.net * Copyright 2001-2003 Refractions Research Inc. * * This is free software; you can redistribute and/or modify it under * the terms of the GNU General Public Licence. See the COPYING file. * ********************************************************************** * * PGRaster - Raster/Coverage Model for PostgreSQL/PostGIS * ********************************************************************** * * File Name: geotiff2pgraster.c * * Author: Xing Lin * * Create Date: 2007-07-22 * * Description: GeoTIFF to PGRaster Converter * * Original Author: Xing Lin * * Maintainer: Xing Lin * * Last Update: 2007-07-22 * * Change Log: * (1) 2007-07-22 CREATED (Xing) * (2) 2007-07-22 Main part of functions (Xing) * (3) 2007-07-30 Finish GeoTIFF loader into GeoTIFF. (Xing) * **********************************************************************/ static char rcsid[] = "$Id: geotiff2pgraster.c 2007-07-22 Xing.Lin $"; #include "../config.h" #include #include #include #include #include #include #include #include #include /* Solaris9 does not provide stdint.h */ /* #include */ #include //GeoTIFF & TIFF library headers #include "geotiff.h" #include "xtiffio.h" #include "geo_normalize.h" #include "geovalues.h" #include "tiffio.h" //Libpq & PGRaster headers #include "libpq-fe.h" #include "getopt.h" #include "pgraster_const.h" #include "pgraster_types.h" #include "pgraster_utils.h" #include "pgraster_version.h" #include "wkb.h" #include "compat.h" /* Used to help create tempate tables names */ #include // for getpid() #include // for getpid() #ifdef HAVE_ICONV_H #include #endif #ifdef __CYGWIN__ #include #endif #ifndef TIFF_SWAB #define TIFF_SWAB 0x0080 /* byte swap file information */ #endif #ifdef USE_ICONV char *utf8(const char *fromcode, char *inputbuf); #endif /* Debug Control Flag */ //#define DEBUG 1 #define VERBOSE 1 /* Buffer length for SQL statement */ #define QUERY_BUF_LENGTH 1024 /* Global data */ PGconn *conn; PGresult *res; TIFF *tif=(TIFF*)0; /* TIFF-level descriptor */ GTIF *gtif=(GTIF*)0; /* GeoKey-level descriptor */ int rowbuflen; char *table, *geotiff_file, *schema, *geotiff_filename, *insertStmt; int keep_fieldname_case = 0; int big_endian = 0; int pgis_major_version; #define MODE_QUIET 0 #define MODE_NORMAL 1 #define MODE_VERBOSE 2 int outmode = 1; // 0=quiet mode; 1=normal mode, 2=verbose mode. bool beginTransaction = FALSE; char *PG_VERSION = "8.2.4\0"; PGRasterMetadata metadata; #ifdef USE_ICONV char *encoding=NULL; #endif /* TIFF Reader Global Data */ byte *pDataBuffer = NULL; short scanlineMode = 0; int nTIFFBytePerPixel = 0; int nTIFFBitsPerPixel = 0; int nRealBytePerPixel = 0; int nRealBitsPerPixel = 0; double dblNoDataValue = -1.0f; uint16 bps = 1; uint16 spp = 1; uint16 mode = PHOTOMETRIC_MINISBLACK; uint16 format = SAMPLEFORMAT_UINT; uint16 planar = PLANARCONFIG_CONTIG; uint16 preMultiplied = 0; uint16 fillorder = 0 ; byte tiff_bigendian = 0; // whether TIFF file uses big endian byte order. byte badIndians = 0; // whether TIFFBigEndian != big_endian byte hasAlpha = 0; // hasAlpha channel??? uint32 rows = 0; // number of TIFF image rows uint32 cols = 0; // number of TIFF image column uint32 bands = 0; // number of raster bands (image channel) uint32 nPyramidDepth = 0; // pyramide depth. /* // If the data buffer is smaller than whole image // we need such variables to navigate through the buffer. unsigned char* scanlineBuffer = NULL; unsigned char* decompressed = NULL; uint32 curRow = 0; uint32 curCol = 0; uint32 curMinRow = 0; uint32 curMaxRow = 0; uint32 curMinCol = 0; uint32 curMaxCol = 0; */ // variables for TIFF transformation CIE to RGB . static float refWhite[3]; static TIFFCIELabToRGB* cielab = 0; #ifdef DEBUG FILE *debug; #endif int getopt(); /* Prototypes */ char * getTableOID(char *schema, char *table); void usage(char* me, int status, FILE* out); int parse_commandline(int ARGC, char **ARGV); static void parse_table(char *spec); static void parse_geotiffname(char *spec); static void exit_nicely(PGconn *conn, int code); int get_pgraster_version(void); int get_postgis_version(void); /* Functions for Data Loading */ /* Step 0. Parse Paramaters and Create Connection*/ void ParseAndCreateConnection(int ARGC, char **ARGV); /* Step 1. Initialization*/ void Initialization(); void _ReadGeoTIFFHeader(); void _PreInsertMetadata(); void _CleanUpDemo(); /* Step 2. LoadData */ void LoadData(); int _BeginReadTIFF(); void _initCIELabConversion(TIFF* pTif); void _EndReadTIFF(); byte* _GetTIFFAt(int col, int row, int band, int *bytereturned); /* Step 3. CreateIndex() */ void CreateIndex(); /* Step 4. LoadSRS */ void LoadSRS(); /* Step 5. UpdateMetadata */ void UpdateMetadata(); /* Step 6. Cleanup */ void Cleanup(); /* Step 7. PrintResultMessage */ void PrintResultMessage(); #ifdef USE_ICONV char *utf8(const char *fromcode, char *inputbuf); #endif void *safe_malloc(size_t size) { void *ret = malloc(size); if ( ! ret ) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error**: Out of virtual memory\n"); exit_nicely(conn, 1); } return ret; } #define malloc(x) safe_malloc(x) /* Main Programe */ int main(int ARGC, char **ARGV) { /* Step 0. Parse Paramaters and Create Connection */ ParseAndCreateConnection(ARGC, ARGV); /* Step 1. Initialization & Create PreparedStatement for Data Insert*/ Initialization(); /* Step 2. Read GeoTIFF Data & Insert into PGRaster via PreparedStatement */ LoadData(); /* Step 3. Run Build GiST Index Statement at the Server-Side */ CreateIndex(); /* Step 4. LoadSRS */ LoadSRS(); /* Step 5. UpdateMetadata */ UpdateMetadata(); /* Step 5. Cleanup & Close Connection*/ Cleanup(); /* Step 6. Function Successfully Return and Send Message To Terminal or Log */ PrintResultMessage(1); return 0; } /******** Implementation of Functions ********/ /* * Return allocate memory. Free after use. */ char * getTableOID(char *schema, char *table) { PGresult *res3; char *query; char *ret; size_t size; size = strlen(table)+256; if ( schema ) size += strlen(schema)+1; query = (char *)malloc(size); if ( schema ) { sprintf(query, "SELECT oid FROM pg_class c, pg_namespace n WHERE c.relnamespace n.oid AND n.nspname = '%s' AND c.relname = '%s'", schema, table); } else { sprintf(query, "SELECT oid FROM pg_class WHERE relname = '%s'", table); } res3 = PQexec(conn, query); free(query); if ( ! res3 || PQresultStatus(res3) != PGRES_TUPLES_OK ) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error** when getting TableOID: %s", PQerrorMessage(conn)); exit_nicely(conn, 1); } if(PQntuples(res3) == 1 ){ ret = strdup(PQgetvalue(res3, 0, 0)); }else if(PQntuples(res3) == 0 ){ if(outmode != MODE_QUIET) fprintf(stdout, "Cannot find relation OID (does table exist?).\n"); PQclear(res3); return NULL; }else{ ret = strdup(PQgetvalue(res3, 0, 0)); if(outmode != MODE_QUIET) printf( "Warning: Multiple relations detected, the program will only dump the first relation.\n"); } PQclear(res3); return ret; } /* * Print the usage information to standard output destination (terminal, file or so on). */ void usage(char* me, int status, FILE* out) { fprintf(out,"PostgreSQL Version: %s \\ PostGIS Version: %s\n", PG_VERSION, POSTGIS_VERSION); fprintf(out,"RCSID: %s RELEASE: %s\n", rcsid, PGRASTER_VERSION); fprintf(out,"USAGE: %s [] geotiff-file [[.]]\n", me); fprintf(out,"\n"); fprintf(out,"OPTIONS:\n"); fprintf(out," -h \n"); fprintf(out," The database host to connect to. If ignored, the default address\n"); fprintf(out," of host will be \"localhost/127.0.0.1\". If no database running on\n"); fprintf(out," this host, the program will exit immediately with a error msg.\n"); fprintf(out," -p \n"); fprintf(out," The port to connect to on the database host. If ignored, the\n"); fprintf(out," standard port number of 5432 will be used.If no database running\n"); fprintf(out," on this port, the program will exit immediately with a error msg.\n"); fprintf(out," -d \n"); fprintf(out," The database to use when connecting. If ignored, the current OS\n"); fprintf(out," username will be used the name of database to login into. If this\n"); fprintf(out," database is refused to login, the program will exit immediately\n"); fprintf(out," with a error message.\n"); fprintf(out," -P \n"); fprintf(out," The password to use when connecting to the database. This option\n"); fprintf(out," can't be ignored. If this password is refused to login, the pro-\n"); fprintf(out," gram will exit immediately with a error message.\n"); fprintf(out," -U \n"); fprintf(out," The username to use when connecting to the database. If ignored,\n"); fprintf(out," the current OS username will be used to login into the database.\n"); fprintf(out," If this username is refused to login, the program will exit\n"); fprintf(out," immediately with a error message.\n"); fprintf(out," -k Keep postgresql identifiers case.\n"); fprintf(out," -?/--help Display this help screen.\n"); fprintf(out," -q Load data in quiet mode.\n"); fprintf(out," -v Load data in verbose mode.\n"); fprintf(out,"\n"); fprintf(out,"(Note that -a, -c, -d and -p are mutually exclusive.)\n"); fprintf(out,"####################### (NOT SUPPORTED YET) #########################\n"); fprintf(out,"(Note that -a, -c, -d and -p are mutually exclusive.)\n"); fprintf(out," -d Drops the database table before creating a new table with the\n"); fprintf(out," data in the geotiff file.\n"); fprintf(out," -a Appends data from the geotiff file into the existing database\n"); fprintf(out," table. Although it is not recommended to store raw data of sever-\n"); fprintf(out," ral raster dataset into one raster data table (sharing). But if\n"); fprintf(out," you insist, you can do it by this option\n"); fprintf(out," -c Creates a new table and populates it from the geotiff file. This\n"); fprintf(out," is the default mode.\n"); fprintf(out," -p Only produces the table creation SQL code, without adding any\n"); fprintf(out," actual data. This can be used if you need to completely seperate\n"); fprintf(out," the table creation and data loading steps.\n"); fprintf(out," (NOT SUPPORTED YET)\n"); fprintf(out," -D Use the PostgreSQL \"dump\" format for the output data. This can\n"); fprintf(out," be combined with -a, -c and -d. It is much faster to load than\n"); fprintf(out," the default \"insert\" SQL format. Use this for very large data\n"); fprintf(out," sets. If this option is set, then the following options could be\n"); fprintf(out," ignored.\n"); fprintf(out," (NOT SUPPORTED YET)\n"); fprintf(out,"####################### (NOT SUPPORTED YET) #########################\n"); fprintf(out,"\n"); exit (status); } /* Parse command line parameters */ int parse_commandline(int ARGC, char **ARGV) { if ( ARGC == 1 ) { usage(ARGV[0], 0, stdout); } int c, curindex; char buf[QUERY_BUF_LENGTH]; memset(buf,'\0', QUERY_BUF_LENGTH); /* Parse command line */ while ((c = getopt(ARGC, ARGV, "h:d:U:p:P:rk:vq")) != EOF){ switch (c) { case 'h': /*setenv("PGHOST", optarg, 1); */ snprintf(buf, 255, "PGHOST=%s", optarg); putenv(strdup(buf)); break; case 'd': /*setenv("PGHOST", optarg, 1); */ snprintf(buf, 255, "PGDATABASE=%s", optarg); putenv(strdup(buf)); break; case 'U': /*setenv("PGUSER", optarg, 1); */ snprintf(buf, 255, "PGUSER=%s", optarg); putenv(strdup(buf)); break; case 'p': /*setenv("PGPORT", optarg, 1); */ snprintf(buf, 255, "PGPORT=%s", optarg); putenv(strdup(buf)); break; case 'P': /*setenv("PGPASSWORD", optarg, 1); */ snprintf(buf, 255, "PGPASSWORD=%s", optarg); putenv(strdup(buf)); break; case 'v': // verbose mode outmode = MODE_VERBOSE; break; case 'q': // quiet mode outmode = MODE_QUIET; break; case 'k': keep_fieldname_case = 1; break; case '?': usage(ARGV[0], 0, stdout); default: return 0; } } curindex=0; for (; optinds * * 1. find # of tabs * 2. make new string * * we dont escape already escaped tabs */ char *result; char *ptr, *optr; int toescape = 0; size_t size; #ifdef USE_ICONV char *utf8str=NULL; if ( encoding ) { utf8str=utf8(encoding, str); if ( ! utf8str ) exit(1); str = utf8str; } #endif ptr = str; while (*ptr) { if ( *ptr == '\t' || *ptr == '\\' ) toescape++; ptr++; } if (toescape == 0) return str; size = ptr-str+toescape+1; result = calloc(1, size); optr=result; ptr=str; while (*ptr) { if ( *ptr == '\t' || *ptr == '\\' ) *optr++='\\'; *optr++=*ptr++; } *optr='\0'; #ifdef USE_ICONV if ( encoding ) free(str); #endif return result; } /* protect_quotes_string from shp_loader*/ char * protect_quotes_string(char *str) { /* * find all quotes and make them \quotes * find all '\' and make them '\\' * * 1. find # of characters * 2. make new string */ char *result; char *ptr, *optr; int toescape = 0; size_t size; #ifdef USE_ICONV char *utf8str=NULL; if ( encoding ) { utf8str=utf8(encoding, str); if ( ! utf8str ) exit(1); str = utf8str; } #endif ptr = str; while (*ptr) { if ( *ptr == '\'' || *ptr == '\\' ) toescape++; ptr++; } if (toescape == 0) return str; size = ptr-str+toescape+1; result = calloc(1, size); optr=result; ptr=str; while (*ptr) { if ( *ptr == '\'' || *ptr == '\\' ) *optr++='\\'; *optr++=*ptr++; } *optr='\0'; #ifdef USE_ICONV if ( encoding ) free(str); #endif return result; } /********************************************************************/ /******************** Functions for Data Loading ********************/ /********************************************************************/ void ParseAndCreateConnection(int ARGC, char **ARGV) { /* Definition and Initialization for local/global variables */ big_endian = is_bigendian(); geotiff_file = NULL; table = NULL; schema = NULL; geotiff_filename = NULL; rowbuflen=100; keep_fieldname_case = 0; conn = NULL; res = NULL; insertStmt = NULL; if ( getenv("ROWBUFLEN") ) rowbuflen=atoi(getenv("ROWBUFLEN")); /* * Make sure dates are returned in ISO * style (YYYY-MM-DD). * This is to allow goodDBFValue() function * to successfully extract YYYYMMDD format * expected in shapefile's dbf file. */ putenv("PGDATESTYLE=ISO"); if ( ! parse_commandline(ARGC, ARGV) ) { if(outmode != MODE_QUIET) { fprintf(stderr, "**ERROR** invalid option or command parameters.\n"); usage(ARGV[0], 2, stderr); } exit_nicely(conn,1); } #ifdef DEBUG if( getenv("PGHOST") != NULL && getenv("PGUSER") != NULL && getenv("PGPORT") != NULL && getenv("PGDATABASE") != NULL && getenv("PGPASSWORD") != NULL ) fprintf(stdout, "PGHOST=%s PGPORT=%s PGDATABASE=%s PGUSER=%s PGPASSWORD=%s\n", getenv("PGHOST"), getenv("PGPORT"), getenv("PGDATABASE"), getenv("PGUSER"), getenv("PGPASSWORD")); fprintf(stdout, "SCHEMA=%s TABLE=%s geotiff_file=%s geotiff_filename=%s\n",schema,table,geotiff_file,geotiff_filename); #endif if(outmode != MODE_QUIET) { fprintf(stdout, "****************************************************************\n"); fprintf(stdout, "****************** GeoTIFF2PGRaster Converter ******************\n"); fprintf(stdout, "****************************************************************\n"); fprintf(stdout, "To PostgreSQL:\n\t host=%s PGPORT=%s database=%s user=%s password=%s schema=%s table=%s \n", getenv("PGHOST"), getenv("PGPORT"), getenv("PGDATABASE"), getenv("PGUSER"), getenv("PGPASSWORD"),schema,table); fprintf(stdout, "To GeoTIFF:\n\tgeotiff_file=%s geotiff_filename=%s\n", geotiff_file,geotiff_filename); } } /* Step 1. Initialization*/ void Initialization() { /* Declare local variable */ char query[QUERY_BUF_LENGTH]; bool bResult = FALSE; if(outmode == MODE_VERBOSE) { fprintf(stdout, "Start initializing data load for GeoTIFF......\n"); } // Initialize PGRasterMetadata object initialize_metadata(&metadata); /* Make a connection to the specified database, and exit on failure */ // All the necessary information is stored in OS Env-Variables conn = PQconnectdb(""); if (PQstatus(conn) != CONNECTION_OK) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error**: when making connection to PostgreSQL: %s", PQerrorMessage(conn)); exit_nicely(conn, 1); } #ifdef DEBUG printf("Connection OK!"); debug = fopen("geotiff2pgraster.log", "wt"); PQtrace(conn, debug); #endif /* DEBUG */ #ifdef DEBUG _CleanUpDemo(); #endif /* * Check Input Schema & Table Valid for Insert??? */ sprintf(query, "SELECT CheckRasterValidity('%s','%s')", schema, table); #ifdef DEBUG printf("Query=%s\n", query); #endif res = PQexecParams(conn, query, 0, NULL, NULL, NULL, NULL, 1); /* ask for binary results */ if (PQresultStatus(res) != PGRES_TUPLES_OK) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error** when check raster name validity: %s", PQerrorMessage(conn)); PQclear(res); exit_nicely(conn,1); } else { bResult = ((byte*)PQgetvalue(res,0,0))[0]; PQclear(res); } if(bResult == FALSE) { if(outmode != MODE_QUIET) fprintf(stdout, "Current schema & table name exist. Please try again using a different pair."); exit_nicely(conn, 1); } else { if(outmode != MODE_QUIET) fprintf(stdout, "GeoTIFF Schema & Name OK. \n"); } /* * Read metadata from GeoTIFF file using libgeotiff & libtiff * & * PreInsert the metdata and generate a new object id */ _ReadGeoTIFFHeader(); /* * Insert metadata into metadata table. * It will be updated after finishing data insertion. */ _PreInsertMetadata(); if(outmode == MODE_VERBOSE) { fprintf(stdout, "Finish data loader initializaiton.\n"); } } void _ReadGeoTIFFHeader() { /* * Open the file, read the GeoTIFF information */ short Model = -1; short RasterModel = -1; // Area or Point int SRID = -1; double xGeoResolution = 0.0; double yGeoResolution = 0.0; double zGeoResolution = 0.0; double xCoordMax = 0.0; double yCoordMax = 0.0; double xCoordMin = 0.0; double yCoordMin = 0.0; PGRasterDataType datatype = DT_UNKNOWN; PGRasterBandType bandtype = BT_UNKNOWN; PGRasterValueType valuetype = VT_IMAGE; if(outmode == MODE_VERBOSE) { fprintf(stdout, "Start reading GeoTIFF & TIFF metadata......\n"); } //Use GeoTIFF & TIFF to Open and Read Header Information. tif=XTIFFOpen(geotiff_file,"r"); if (!tif) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error**: when openning TIFF files.\n"); exit_nicely(conn, 1); } gtif = GTIFNew(tif); if (!gtif) { if(outmode != MODE_QUIET) fprintf(stderr,"**Error**: failed to open GeoTIFF object.\n"); exit_nicely(conn, 1); } //Read Image Metdata Data { uint16 tif_magic = 0; FILE* ftifftemp = fopen(geotiff_file, "r"); // Judge TIFF Endian?? if(ftifftemp == NULL) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error** when open another instance of TIFF file\n"); exit_nicely(conn, 1); } fseek(ftifftemp,0,SEEK_SET); fread(&tif_magic,sizeof(tif_magic),1, ftifftemp); fclose(ftifftemp); ftifftemp = NULL; if( tif_magic == TIFF_BIGENDIAN) tiff_bigendian = 1; else tiff_bigendian = 0; #ifdef DEBUG fprintf(stderr, "TIFF Big Endian: %d, Big Endian: %d\n", tiff_bigendian, big_endian); #endif if(tiff_bigendian != big_endian) badIndians = 1; else badIndians = 0; //Read TIFF Header Data TIFFGetField(tif, TIFFTAG_BITSPERSAMPLE, &bps); TIFFGetField(tif, TIFFTAG_SAMPLESPERPIXEL, &spp); TIFFGetField(tif, TIFFTAG_PHOTOMETRIC, &mode); TIFFGetField(tif, TIFFTAG_SAMPLEFORMAT, &format); TIFFGetField(tif, TIFFTAG_PLANARCONFIG, &planar); TIFFGetField(tif, TIFFTAG_MATTEING, &preMultiplied); TIFFGetField(tif, TIFFTAG_FILLORDER, &fillorder); TIFFGetField(tif, TIFFTAG_IMAGEWIDTH, &cols); TIFFGetField(tif, TIFFTAG_IMAGELENGTH, &rows); //Calculate nTIFFBytePerPixel (Not so simple) nTIFFBytePerPixel = bps*spp/8; nTIFFBitsPerPixel = bps*spp; if((mode==PHOTOMETRIC_RGB) && (spp==4)) hasAlpha = 1; else hasAlpha = 0; if(bps == 1 && (mode == PHOTOMETRIC_MINISBLACK || mode == PHOTOMETRIC_MINISWHITE)) { //Black and White datatype = DT_1BIT; if(spp > 1) { bandtype = BT_MULTI; valuetype = VT_IMAGE; } else { bandtype = BT_SINGLE; valuetype = VT_NOMINAL; } bands = spp; nRealBytePerPixel = -8; nRealBitsPerPixel = bps*bands; } else if(mode == PHOTOMETRIC_RGB) { bands = 1; bandtype = BT_SINGLE; valuetype = VT_IMAGE; if(spp == 3) // RGB { bandtype = BT_RGB; datatype = DT_24BIT_RGB; } else { bandtype = BT_RGBA; datatype = DT_32BIT_RGBA; } nRealBytePerPixel = spp; nRealBitsPerPixel = spp * 8; } else if(mode == PHOTOMETRIC_CIELAB || mode == PHOTOMETRIC_PALETTE) { //GeoTIFF RGA or CYK mode bands = 1; bandtype = BT_SINGLE; valuetype = VT_IMAGE; bandtype = BT_RGB; datatype = DT_24BIT_RGB; nRealBytePerPixel = 3; nRealBytePerPixel = 8 * nRealBytePerPixel; } else { // Determind by format. bands = spp; if(spp > 1) bandtype = BT_MULTI; else bandtype = BT_SINGLE; nRealBitsPerPixel = bps*spp; if(bps < 8) { nRealBytePerPixel = (-1) * bps; } else nRealBytePerPixel = bps / 8; valuetype = VT_RATIO; switch (format) { case SAMPLEFORMAT_VOID: { if(bps == 8) datatype = DT_8BIT_U; else if(bps == 16) datatype = DT_16BIT_U; else if(bps == 32) datatype = DT_32BIT_REAL; else if(bps == 64) datatype = DT_64BIT_REAL; else { if(outmode != MODE_QUIET) fprintf(stderr , "**Error**: Unsupported TIFF sample format.\n"); exit_nicely(conn,1); } break; } case SAMPLEFORMAT_INT: { if(bps == 8) datatype = DT_8BIT_S; else if(bps == 16) datatype = DT_16BIT_S; else if(bps == 32) datatype = DT_32BIT_S; else {// Read spatial extent information from PGRASTER_METADATA if(outmode != MODE_QUIET) fprintf(stderr , "**Error**: Unsupported TIFF sample format.\n"); exit_nicely(conn,1); } break; } case SAMPLEFORMAT_UINT: { if(bps == 8) datatype = DT_8BIT_U; else if(bps == 16) datatype = DT_16BIT_U; else if(bps == 32) datatype = DT_32BIT_U; else { if(outmode != MODE_QUIET) fprintf(stderr, "**Error**: Unsupported TIFF sample format.\n"); exit_nicely(conn,1); } break; } case SAMPLEFORMAT_IEEEFP: { if(bps != 32) { if(outmode != MODE_QUIET) fprintf(stderr , "**Error**: Unsupported TIFF sample format.\n"); exit_nicely(conn,1); } else { datatype = DT_32BIT_REAL; } break; } case SAMPLEFORMAT_COMPLEXINT: case SAMPLEFORMAT_COMPLEXIEEEFP: default: { if(outmode != MODE_QUIET) fprintf(stderr , "**Error**: Unsupported TIFF sample format.\n"); exit_nicely(conn,1); } } } } //Start Reading Geo-Metedata Data { GTIFDefn defn; GTIFKeyGet( gtif, GTRasterTypeGeoKey, (void*)&RasterModel, 0, 1 ); #ifdef DEBUG if(RasterModel == RasterPixelIsArea) { printf("RasterModel : RasterPixelIsArea = %d\n", RasterModel); } else if(RasterModel == RasterPixelIsPoint) { printf("RasterModel: RasterPixelIsPoint = %d\n", RasterModel); } #endif if( GTIFGetDefn( gtif, &defn ) ) { Model = defn.Model; SRID = defn.PCS; //Get MBR in PCS / GCS xCoordMax = cols; xCoordMin = 0.0; yCoordMax = 0.0; yCoordMin = rows; if( !GTIFImageToPCS( gtif, &xCoordMax, &yCoordMax ) ) { xCoordMax = 0.0f; yCoordMax = 0.0f; } if( !GTIFImageToPCS( gtif, &xCoordMin, &yCoordMin ) ) { xCoordMax = 0.0f; yCoordMax = 0.0f; } if(yCoordMax < yCoordMin) { double temp = yCoordMax; yCoordMax = yCoordMin; yCoordMin = temp; } if(xCoordMax < xCoordMin) { double temp = xCoordMax; xCoordMax = xCoordMin; xCoordMin = temp; } //Calculate X & Y Spatial Resolution if(RasterModel == RasterPixelIsArea) { xGeoResolution = (xCoordMax - xCoordMin) / (cols); yGeoResolution = (yCoordMax - yCoordMin) / (rows); zGeoResolution = 0.0f; } else if (RasterModel == RasterPixelIsPoint) { xGeoResolution = (xCoordMax - xCoordMin) / (cols-1); yGeoResolution = (yCoordMax - yCoordMin) / (rows-1); zGeoResolution = 0.0f; } #ifdef DEBUG printf("xGeoResolution: %lf\n", xGeoResolution); printf("yGeoResolution: %lf\n", yGeoResolution); printf("zGeoResolution: %lf\n", zGeoResolution); #endif /* * Check if current SRID exist. If it is a user-defined project/geographic * coordinate system, please insert into spatial_ref_sys. */ { char query[QUERY_BUF_LENGTH]; uint32 nCount = 0; sprintf(query,"SELECT count(*)::INTEGER FROM \"%s\".\"spatial_ref_sys\" WHERE srid=%d",schema,SRID); res = PQexecParams( conn, query, 0, NULL, NULL, NULL, NULL, 1); if (PQresultStatus(res) != PGRES_TUPLES_OK) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error** when check whether SRID exist: %s. (Does table spatial_ref_sys existed?)", PQerrorMessage(conn)); PQclear(res); exit_nicely(conn,1); } else { nCount = ((uint32*)PQgetvalue(res,0,0))[0]; nCount = (uint32)(_ntohl(nCount)); PQclear(res); } if(nCount <= 0) { //XXX:NEW SRID, Need to Insert into spatial_ref_sys tables } //else, current srid existed. Good. } } else // NONE SRID { SRID = -1; } } //End Reading Geo-Metedata Data //Push Metata Into PGRasterMetadata object. { WKSEnvelope wksMBR; wksMBR.maxX = xCoordMax; wksMBR.minX = xCoordMin; wksMBR.maxY = yCoordMax; wksMBR.minY = yCoordMin; metadata.rasterObjectID = -1; metadata.rasterDimensions = 2; metadata.rasterBandType = BT_MULTI; metadata.rasterDataType = datatype; metadata.rasterValueType = VT_NOMINAL; metadata.rasterSchema = schema; metadata.rasterDataTable = table; metadata.rasterBandCount = bands; metadata.rasterRowCount = rows; metadata.rasterColumnCount = cols; metadata.rasterCellDepth = nRealBitsPerPixel; #ifdef DEBUG fprintf(stderr, "nRealBitsPerPixel: %d\n",nRealBitsPerPixel); #endif metadata.blockSizeBands = 1; metadata.blockCompression = CT_NONE; metadata.SRID = SRID; metadata.spatialResolutionY = yGeoResolution; metadata.spatialResolutionX = xGeoResolution; metadata.spatialResolutionZ = zGeoResolution; if(SRID != -1) metadata.geoReferenced = TRUE; metadata.spatialExtent = G_MBRToGeometry(&wksMBR); if(big_endian) metadata.spatialExtent->obj.byteOrder = wkbXDR; else metadata.spatialExtent->obj.byteOrder = wkbNDR; } if(outmode == MODE_VERBOSE) { fprintf(stdout, "Finish reading GeoTIFF & TIFF metadata.\n"); } } void _PreInsertMetadata() { char query[QUERY_BUF_LENGTH]; int fldLen = 0; char *pszValue = NULL; if(outmode == MODE_VERBOSE) { fprintf(stdout, "Updating PGRaster_Metadata ......\n"); } /* * Before data inserting: Begin a transaction * (a cursor can only be defined inside a transaction block) */ res=PQexec(conn, "BEGIN"); if ( ! res || PQresultStatus(res) != PGRES_COMMAND_OK ) { printf( "%s", PQerrorMessage(conn)); exit_nicely(conn, 1); } PQclear(res); beginTransaction = TRUE; /* Insert a new record in the PGRASTER_METADATA table * If error, such as there is no PGRASTER_METADATA table exist, * this programe will exit abnormally. * Using SQL functions. */ memset(query,'\0', QUERY_BUF_LENGTH); sprintf(query,"SELECT AddRasterTableFull('%s','%s',%f,%f,%d,%d,%d,%d,%d,%d,%d,%d,%d);", schema, table, metadata.spatialResolutionX, metadata.spatialResolutionY, metadata.rasterBandCount, metadata.rasterRowCount, metadata.rasterColumnCount, metadata.rasterDataType, metadata.rasterBandType, metadata.rasterValueType, metadata.blockBandInterleaving, metadata.blockCompression, metadata.SRID); res = PQexecParams(conn, query, 0, NULL, NULL, NULL, NULL, 0); if ( ! res || PQresultStatus(res) != PGRES_TUPLES_OK ) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error** when execute AddRasterTableFull at server side: %s\n.", PQerrorMessage(conn)); exit_nicely(conn, 1); } fldLen = PQgetlength(res,0,0); pszValue = PQgetvalue(res,0,0); metadata.name = malloc(fldLen * sizeof(char)+1); memset(metadata.name,'\0',sizeof(char)*fldLen+1); memcpy(metadata.name,(char*)pszValue,fldLen*sizeof(char)); PQclear(res); memset(query, '\0', QUERY_BUF_LENGTH); sprintf(query,"SELECT rasterObjectID FROM pgraster_metadata WHERE name='%s'", metadata.name); #ifdef DEBUG printf("Name=%s\n",metadata.name); #endif res = PQexecParams(conn, query, 0, NULL, NULL, NULL, NULL, 0); if ( ! res || PQresultStatus(res) != PGRES_TUPLES_OK ) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error** when trying to retrieval objectid from server side: %s\n.",PQerrorMessage(conn)); exit_nicely(conn, 1); } if(strcmp(PQgetvalue(res,0,0),"fail")==0) { if(outmode != MODE_QUIET) fprintf(stderr, "Current schema & table name exist. Please try again using a different pair. PQerrorMessage: %s\n", PQerrorMessage(conn)); exit_nicely(conn, 1); } else metadata.rasterObjectID = atoi((char*)PQgetvalue(res,0,0)); PQclear(res); { const char *paramValues[1]; int paramLengths[1]; int paramFormats[1]; paramValues[0] = (char*)(metadata.spatialExtent); paramLengths[0] = G_GetSize(metadata.spatialExtent); paramFormats[0] = 1; //Binary memset(query, '\0' , QUERY_BUF_LENGTH); sprintf(query,"UPDATE \"%s\".\"pgraster_metadata\" SET spatialextent=GeomFromWKB($1,%d) WHERE rasterObjectID=%Ld", schema, metadata.SRID, metadata.rasterObjectID); #ifdef DEBUG printf("Update Spatial Extent Query: %s\n",query); #endif res = PQexecParams(conn, query, 1, NULL, paramValues, paramLengths, paramFormats, 0); if ( PQresultStatus(res) != PGRES_COMMAND_OK ) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error** when trying to update spatial extent at server side: %s\n.",PQerrorMessage(conn)); exit_nicely(conn, 1); } PQclear(res); } { // Create Prepared Insert Statement for Raster Data Insertion insertStmt = (char*)malloc(strlen(schema) + strlen(table) + 1); memset(insertStmt, '\0' , strlen(schema) + strlen(table) + 1); memcpy(insertStmt, schema, strlen(schema)); memcpy(insertStmt+strlen(schema), table, strlen(table)); memset(query, '\0' , QUERY_BUF_LENGTH); sprintf(query,"INSERT INTO \"%s\".\"%s\"(rasterobjectid,pyramidlevel,bandblocknumber,rowblocknumber,columnblocknumber,blockbandsize,blockrowsize,blockcolumnsize,blockmbr,datablock) Values($1,$2,$3,$4,$5,$6,$7,$8,GeomFromWKB($9),$10)", schema, table); res = PQprepare( conn, insertStmt, query, 10, NULL); if ( PQresultStatus(res) != PGRES_COMMAND_OK ) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error** when trying to create prepared statement for data insertion at server side: %s\n.",PQerrorMessage(conn)); exit_nicely(conn, 1); } PQclear(res); } if(outmode == MODE_VERBOSE) { fprintf(stdout, "Finish writing PGRaster_Metadata first time.\n"); } } /* * Functions used to clean up temp tables by demo functions. * * Will only be invoked while testing/debuging mode. */ void _CleanUpDemo() { char query[QUERY_BUF_LENGTH]; memset(query, '\0' , QUERY_BUF_LENGTH); printf("Executing CleanUpDemo()!!!!\n"); sprintf(query,"SELECT DropRasterTableByName('%s','%s')", schema,table); printf("Query=%s\n",query); res = PQexecParams(conn, query, 0, NULL, NULL, NULL, NULL, 0); if ( ! res || PQresultStatus(res) != PGRES_TUPLES_OK ) { printf("Error when trying to clean up testing tables at server side: %s\n.",PQerrorMessage(conn)); } PQclear(res); } /* Step 2. LoadData */ void LoadData() { uint32 blockRowSize = metadata.blockSizeRows; uint32 blockColumnSize = metadata.blockSizeColumns; uint32 blockBandSize = metadata.blockSizeBands; byte *data = NULL; uint32 datalength = 0; uint32 objectid = metadata.rasterObjectID; WKBGeometry* blockMBR = NULL; WKSEnvelope wksMBR; const char *paramValues[10]; int paramLengths[10]; int paramFormats[10]; dblNoDataValue = metadata.nodataValue; int i; if(outmode == MODE_VERBOSE) { fprintf(stdout, "Start loading GeoTIFF data into PostgreSQL/PostGIS.\n"); } nPyramidDepth = 0; //Calculate Pyramid Depth nPyramidDepth = CalculatePyramidDepth(rows,cols,blockRowSize,blockColumnSize); //Begin TIFF Reader if(!_BeginReadTIFF()) { fprintf(stderr, "Error when initialize tiff reader.\n"); exit_nicely(conn,1); } //Read Data in Scanline Line Modes & Insert into RasterData Table. //By Default for(i = 0 ; i < bands ; i++) { uint32 nPyramidLevel = 1; while(1) { uint32 nBlockIndex_H = 1 ; uint32 nBlockIndex_V = 1 ; uint32 nBlockIndex_B = i + 1; long nColumns_P = 0; long nRows_P = 0; long nBlockHorizontal_P; long nBlockVertical_P; if(nPyramidLevel > nPyramidDepth) break; if(cols%(long)pow(2,nPyramidLevel-1) == 0) { nColumns_P = (long)(cols /(long)pow(2,nPyramidLevel-1)); } else { nColumns_P = (long)(cols /(long)pow(2,nPyramidLevel-1) + 1); } if(rows % (long)pow(2,nPyramidLevel-1) == 0) { nRows_P = (long)(rows/(long)pow(2,nPyramidLevel-1)); } else { nRows_P = (long)(rows/(long)pow(2,nPyramidLevel -1)+1); } nBlockHorizontal_P = (int)((nColumns_P%blockColumnSize) == 0? (nColumns_P/blockColumnSize):(nColumns_P/blockColumnSize + 1)); nBlockVertical_P = (int)((nRows_P%blockRowSize) == 0? (nRows_P/blockRowSize):(nRows_P/blockRowSize + 1)); while(1) { uint32 nBufferWidth = 0; uint32 nBufferHeight = 0; int row,col; int nCurByte = 0; int nCurBit = 0; // Use when bps < 8 long nMaxCol = 0; long nMinCol = 0; long nMaxRow = 0; long nMinRow = 0; if(nBlockIndex_H > nBlockHorizontal_P || nBlockIndex_V > nBlockVertical_P) break; if(nBlockIndex_H == nBlockHorizontal_P) nBufferWidth = (int)(nColumns_P - blockColumnSize*(nBlockIndex_H -1)); else nBufferWidth = blockColumnSize; if(nBlockIndex_V == nBlockVertical_P) nBufferHeight = (int)(nRows_P- blockRowSize*(nBlockIndex_V -1)); else nBufferHeight = blockRowSize; datalength = nBufferWidth * nBufferHeight * nRealBitsPerPixel; if(datalength % 8 != 0) datalength = (datalength / 8) +1; else datalength = datalength / 8; data = malloc(sizeof(char)* datalength); for( row = 0 ; row < nBufferHeight;row++) { for(col = 0 ; col < nBufferWidth;col++) { long nCol_Local = ((long)pow(2,nPyramidLevel-1))*((nBlockIndex_H-1)*blockColumnSize+col); long nRow_Local = ((long)pow(2,nPyramidLevel-1))*((nBlockIndex_V-1)*blockRowSize+row); int nLen = -1; //Read Data Here byte* pData = _GetTIFFAt(nCol_Local, nRow_Local,i, &nLen); //Store data to Byte buffer (should we consider the BSQ/BIL/BIP??) if(nLen < 8) { if(nLen == 1) { //black white data[nCurByte] = SetSingleBitAt(data[nCurByte],pData[0],nCurBit); nCurBit += 1; if(nCurBit == 8) { nCurByte += 1; nCurBit = 0; } } else if(nLen == 2) { //2 bits Grey Level data[nCurByte] = SetMultiBitAt(data[nCurByte],pData[0],2,nCurBit); nCurBit += 2; if(nCurBit == 8) { nCurByte += 1; nCurBit = 0; } } else if(nLen == 4) { //4 bits Grey Level data[nCurByte] = SetMultiBitAt(data[nCurByte],pData[0],4,nCurBit); nCurBit += 4; if(nCurBit == 8) { nCurByte += 1; nCurBit = 0; } } } else { memcpy(data+nCurByte,pData,nLen/8); nCurByte += (nLen / 8); } free(pData); } } // Calculate Max/Min Column & Row & Make block MBR nMaxCol = ((long)pow(2,nPyramidLevel-1))*((nBlockIndex_H-1)*blockColumnSize+nBufferWidth) - 1; if(nMaxCol >= cols) nMaxCol = cols - 1; nMinCol = ((long)pow(2,nPyramidLevel-1))*((nBlockIndex_H-1)*blockColumnSize); nMaxRow = ((long)pow(2,nPyramidLevel-1))*((nBlockIndex_V-1)*blockRowSize+nBufferHeight) - 1; if(nMaxRow >= rows) nMaxRow = rows - 1; nMinRow = ((long)pow(2,nPyramidLevel-1))*((nBlockIndex_V-1)*blockRowSize); if(outmode == MODE_VERBOSE) fprintf(stdout, "Pyramid:%d, Block:(%d,%d,%d), MBR(%ld %ld,%ld %ld)\n", (int)nPyramidLevel, (int)nBlockIndex_H, (int)nBlockIndex_V, (int)nBlockIndex_B, nMinCol, nMinRow, nMaxCol, nMaxRow); wksMBR.maxX = nMaxCol; wksMBR.minX = nMinCol; wksMBR.maxY = nMaxRow; wksMBR.minY = nMinRow; blockMBR = G_MBRToGeometry(&wksMBR); // Insert Data Record in Raster Data Table // 1. rasterobjectid { //Number in network byte order uint32 n_objectid64, n_pyramidlevel, n_bandblocknumber, n_rowblocknumber, n_columnblocknumber, n_blockbandsize, n_blockrowsize, n_blockcolumnsize; //n_objectid = _htonl(objectid); n_objectid64 = _htonl((_uint64)objectid); paramValues[0] = (char*)(&(n_objectid64)); paramLengths[0] = sizeof(n_objectid64); paramFormats[0] = 1; // Binary / Endian?? // 2. pyramidlevel n_pyramidlevel = _htonl(nPyramidLevel); paramValues[1] = (char*)(&(n_pyramidlevel)); paramLengths[1] = sizeof(n_pyramidlevel); paramFormats[1] = 1; // Binary / Endian?? // 3. bandblocknumber n_bandblocknumber = _htonl(nBlockIndex_B); paramValues[2] = (char*)(&(n_bandblocknumber)); paramLengths[2] = sizeof(n_bandblocknumber); paramFormats[2] = 1; // Binary / Endian?? // 4. rowblocknumber n_rowblocknumber = _htonl(nBlockIndex_V); paramValues[3] = (char*)(&(n_rowblocknumber)); paramLengths[3] = sizeof(n_rowblocknumber); paramFormats[3] = 1; // Binary / Endian?? // 5. columnblocknumber n_columnblocknumber = _htonl(nBlockIndex_H); paramValues[4] = (char*)(&(n_columnblocknumber)); paramLengths[4] = sizeof(n_columnblocknumber); paramFormats[4] = 1; // Binary / Endian?? // 6. blockbandsize n_blockbandsize = _htonl(blockBandSize); paramValues[5] = (char*)(&(n_blockbandsize)); paramLengths[5] = sizeof(n_blockbandsize); paramFormats[5] = 1; // Binary / Endian?? // 7. blockrowsize n_blockrowsize = _htonl(nBufferHeight); paramValues[6] = (char*)(&(n_blockrowsize)); paramLengths[6] = sizeof(n_blockrowsize); paramFormats[6] = 1; // Binary / Endian?? // 8. blockcolumnsize n_blockcolumnsize = _htonl(nBufferWidth); paramValues[7] = (char*)(&(n_blockcolumnsize)); paramLengths[7] = sizeof(n_blockcolumnsize); paramFormats[7] = 1; // Binary / Endian?? // 9. blockmbr paramValues[8] = (char*)blockMBR; paramLengths[8] = G_GetSize(blockMBR); paramFormats[8] = 1; // Binary?? //Should we compressed the data here before insertion. // 10. datablock paramValues[9] = (char*)data; paramLengths[9] = datalength; paramFormats[9] = 1; // Binary?? // Execute Prepared Statement res = PQexecPrepared( conn, insertStmt, 10, paramValues, paramLengths, paramFormats, 0 ); if ( PQresultStatus(res) != PGRES_COMMAND_OK ) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error** when trying to insert raster data at server side: %s\n.",PQerrorMessage(conn)); exit_nicely(conn, 1); } PQclear(res); } // Free Memory free(data); G_DestroyGeometry(blockMBR); blockMBR = NULL; // Update Old Trailers if(nBlockIndex_H == nBlockHorizontal_P && nBlockIndex_V < nBlockVertical_P) { nBlockIndex_H = 1; nBlockIndex_V ++ ; } else { nBlockIndex_H ++; } } // Update nPyramidLevel nPyramidLevel++ ; } } if(outmode == MODE_VERBOSE) { fprintf(stdout, "Finish GeoTIFF data loading.\n"); } } int _BeginReadTIFF() { int result = 0; uint32 bufferSize; if(!TIFFIsTiled(tif)) { // Process in Strips. unsigned long imageOffset; uint32 stripMax, stripCount; tsize_t stripSize; stripSize = TIFFStripSize (tif); stripMax = TIFFNumberOfStrips (tif); imageOffset = 0; bufferSize = stripMax * stripSize; if((pDataBuffer = (byte*) malloc(bufferSize)) == NULL) { if(outmode != MODE_QUIET) fprintf(stderr, "Could not allocate enough memory for the uncompressed TIFF image.\n"); return 0; } if( mode == PHOTOMETRIC_CIELAB ) _initCIELabConversion(tif); TIFFSetField(tif, TIFFTAG_ORIENTATION, ORIENTATION_TOPLEFT); for (stripCount = 0; stripCount < stripMax; stripCount++) { if((result = TIFFReadEncodedStrip (tif, stripCount, pDataBuffer + imageOffset, stripSize)) == -1) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error** read TIFF data in strip of number %d\n", (unsigned int)stripCount); return 0; } imageOffset += result; } return 1; } else { // Process in Tiles uint32 width = cols, height = rows; /* image width & height */ uint32 tile_width, tile_height; byte *raster = NULL; uint32 row, col; uint32 rows_read, cols_read; TIFFGetField(tif, TIFFTAG_TILEWIDTH, &tile_width); TIFFGetField(tif, TIFFTAG_TILELENGTH, &tile_height); bufferSize = width * height * nTIFFBytePerPixel; if((pDataBuffer = (byte*) malloc(bufferSize)) == NULL) { if(outmode != MODE_QUIET) fprintf(stderr, "Could not allocate enough memory for the uncompressed TIFF image.\n"); return 0; } raster = _TIFFmalloc(tile_width * tile_height * nTIFFBytePerPixel); if (raster == 0) { if(outmode != MODE_QUIET) fprintf(stderr, "Could not allocate enough memory for the TILE buffer.\n"); return 0; } for (row = 0; row < height; row += tile_height) { if((row+tile_height) < height) rows_read = tile_height; else rows_read = height - row; for (col = 0; col < width; col += tile_width) { uint32 i; // Read TIFFReadTile(tif, raster, col, row, 0, 0); if((col+tile_width) < width) cols_read = tile_width; else cols_read = width - col; // Put Data Into Data Buffer. for( i = 0 ; i luminance matrix */ { 3.2410F, -1.5374F, -0.4986F }, { -0.9692F, 1.8760F, 0.0416F }, { 0.0556F, -0.2040F, 1.0570F } }, 100.0F, 100.0F, 100.0F, /* Light o/p for reference white */ 255, 255, 255, /* Pixel values for ref. white */ 1.0F, 1.0F, 1.0F, /* Residual light o/p for black pixel */ 2.4F, 2.4F, 2.4F, /* Gamma values for the three guns */ }; void _initCIELabConversion(TIFF* pTif) { float *whitePoint; if( !cielab) { cielab = (TIFFCIELabToRGB *)_TIFFmalloc(sizeof(TIFFCIELabToRGB)); } TIFFGetFieldDefaulted(pTif, TIFFTAG_WHITEPOINT, &whitePoint); refWhite[1] = 100.0F; refWhite[0] = whitePoint[0] / whitePoint[1] * refWhite[1]; refWhite[2] = (1.0F - whitePoint[0] - whitePoint[1]) / whitePoint[1] * refWhite[1]; if (TIFFCIELabToRGBInit(cielab, &display_sRGB, refWhite) < 0) { if(outmode != MODE_QUIET) fprintf(stderr, "Failed to initialize CIE L*a*b*->RGB conversion state."); } } byte* _GetTIFFAt(int col, int row, int band, int *bitsreturned) { *bitsreturned = 0; byte* buffer = NULL; // By Default, we assume that TIFF uses a BSQ packing schema. (Only work for bps=n*8) int nIndex = (cols * row + col ) * nTIFFBytePerPixel * spp; // if there is no data buffer available or we need to move on the cursor if( (spp==3) || (spp==4) ) // Multi-Channel { if(bps == 8 && mode==PHOTOMETRIC_RGB) { *bitsreturned = spp * 8; buffer = (byte*)malloc(spp*sizeof(char)); memcpy(buffer, (byte*)(&(pDataBuffer[nIndex])), spp); } else if (bps==16 && mode==PHOTOMETRIC_RGB) { *bitsreturned = spp*16; buffer = (byte*)malloc(spp*sizeof(uint16)); if(spp == 3 ) // RGB { ((uint16*)buffer)[0] = ((uint16*)(&(pDataBuffer[nIndex])))[0]; // Red ((uint16*)buffer)[1] = ((uint16*)(&(pDataBuffer[nIndex])))[1]; // Green ((uint16*)buffer)[2] = ((uint16*)(&(pDataBuffer[nIndex])))[2]; // Blue if(tiff_bigendian) { // inverse byte order ((uint16*)buffer)[0] = _ntohs(((uint16*)buffer)[0]); ((uint16*)buffer)[1] = _ntohs(((uint16*)buffer)[1]); ((uint16*)buffer)[2] = _ntohs(((uint16*)buffer)[2]); } } else if(spp == 4) // RGBA { ((uint16*)buffer)[0] = ((uint16*)(&(pDataBuffer[nIndex])))[0]; // Red ((uint16*)buffer)[1] = ((uint16*)(&(pDataBuffer[nIndex])))[1]; // Green ((uint16*)buffer)[2] = ((uint16*)(&(pDataBuffer[nIndex])))[2]; // Blue ((uint16*)buffer)[3] = ((uint16*)(&(pDataBuffer[nIndex])))[3]; // Alpha if(tiff_bigendian) { // inverse byte order ((uint16*)buffer)[0] = _ntohs(((uint16*)buffer)[0]); ((uint16*)buffer)[1] = _ntohs(((uint16*)buffer)[1]); ((uint16*)buffer)[2] = _ntohs(((uint16*)buffer)[2]); ((uint16*)buffer)[3] = _ntohs(((uint16*)buffer)[3]); } } } else { if( format != SAMPLEFORMAT_UINT && format != SAMPLEFORMAT_INT && format != SAMPLEFORMAT_IEEEFP ) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error**: Unsupported TIFF sample format.\n"); exit_nicely(conn,1); } if( format == SAMPLEFORMAT_IEEEFP && bps != 32 ) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error**: Unsupported TIFF sample format.\n"); exit_nicely(conn,1); } if( planar != PLANARCONFIG_CONTIG ) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error**: Unsupported PLANARCONFIG_CONTIG sample format.\n"); exit_nicely(conn,1); } switch(mode) { case PHOTOMETRIC_RGB: { // Return RGBA data. // Can't happen here. It should happened at previous IF...THEN section. if(outmode != MODE_QUIET) fprintf(stdout, "**Error**: Can't have such TIFF sample format.\n"); break; } case PHOTOMETRIC_CIELAB: { // Coverted into RGBA mode and return char* p = (char*) (&(pDataBuffer[nIndex])); // a and b are signed, L is not... float X=0, Y=0, Z=0; uint32 r=0, g=0, b=0; TIFFCIELabToXYZ(cielab, ((unsigned char* )p)[0], // *NO* sign extension for L ! p[1], p[2], &X, &Y, &Z); TIFFXYZToRGB(cielab, X, Y, Z, &r, &g, &b); if(bps == 8) { *bitsreturned = 3 * 8; buffer = (byte*)malloc(3*sizeof(char)); buffer[0] = r; buffer[1] = g; buffer[2] = b; } else if(bps == 16) { *bitsreturned = 3 * 16; buffer = (byte*)malloc(6); ((uint16*)buffer)[0] = (uint16)r; // Red ((uint16*)buffer)[1] = (uint16)g; // Green ((uint16*)buffer)[2] = (uint16)b; // Blue } else { if(outmode != MODE_QUIET) fprintf(stderr, "**Error**: Unsupported sample format sample format.\n"); exit_nicely(conn,1); } break; } default: { *bitsreturned = bps; if(band >= spp) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error**: Unsupported PLANARCONFIG_CONTIG sample format.\n"); exit_nicely(conn,1); } // Directly return data if(bps < 8) { //less than 1 bytes int nIndexInByte = (nTIFFBitsPerPixel* (cols * row + col )) / 8; int nOffsetInBits = (nTIFFBitsPerPixel* (cols * row + col )) % 8; buffer = (byte*)malloc(sizeof(char)); byte value = ((byte*)pDataBuffer)[nIndexInByte]; if(bps == 1) { buffer[0] = GetSingleBitAt(value, nOffsetInBits); } else if(bps == 2 || bps == 4) { byte value = ((byte*)pDataBuffer)[nIndexInByte]; buffer[0] = GetMultiBitAt(value, nOffsetInBits,bps); } else { if(outmode != MODE_QUIET) fprintf(stderr, "**Error**: Unsupported PLANARCONFIG_CONTIG sample format.\n"); exit_nicely(conn,1); } } else { if(bps == 8) { buffer = (byte*)malloc(sizeof(byte)); buffer[0] = ((byte*)pDataBuffer)[nIndex+band]; // TIFF By Default BSQ. } else if(bps == 16) { buffer = (byte*)malloc(sizeof(_uint16)); *((_uint16*)buffer) = ((_uint16*)pDataBuffer)[nIndex+band]; // TIFF By Default BSQ. if(tiff_bigendian) { *((_uint16*)buffer) = _ntohs(*((_uint16*)buffer)); } } else if(bps == 32) { buffer = (byte*)malloc(sizeof(_uint32)); *((_uint32*)buffer) = ((_uint32*)pDataBuffer)[nIndex+band]; // TIFF By Default BSQ. if(tiff_bigendian) { *((_uint32*)buffer) = _ntohl(*((_uint32*)buffer)); } } else if(bps == 64) { buffer = (byte*)malloc(sizeof(_int64)); *((_int64*)buffer) = ((_int64*)pDataBuffer)[nIndex+band]; // TIFF By Default BSQ. if(tiff_bigendian) { *((_int64*)buffer) = _ntohll(*((_int64*)buffer)); } } else { if(outmode != MODE_QUIET) fprintf(stderr, "**Error**: Unsupported PLANARCONFIG_CONTIG sample format.\n"); exit_nicely(conn,1); } } break; } } } } else if( spp==1 ) // Single Channel { int colors = (bps <8) ? 1<>8; buffer[1] = green[c]>>8; buffer[2] = blue[c]>>8; } } else if( bps == 1 && ((mode==PHOTOMETRIC_MINISWHITE) || (mode==PHOTOMETRIC_MINISBLACK)) ) // Black White Pictures { int nIndexInByte = (nTIFFBitsPerPixel* (cols * row + col )) / 8; int nOffsetInBits = (nTIFFBitsPerPixel* (cols * row + col )) % 8; buffer = (byte*)malloc(sizeof(char)); byte c = ((byte*)pDataBuffer)[nIndexInByte]; if(mode == PHOTOMETRIC_MINISWHITE) c = ~c; if(fillorder != FILLORDER_MSB2LSB) { // We need to SWAP bits -- ABCDEFGH becomes HGFEDCBA byte tempbyte = 0; if(c & 128) tempbyte += 1; if(c & 64) tempbyte += 2; if(c & 32) tempbyte += 4; if(c & 16) tempbyte += 8; if(c & 8) tempbyte += 16; if(c & 4) tempbyte += 32; if(c & 2) tempbyte += 64; if(c & 1) tempbyte += 128; c = tempbyte; } // Pick Up the right bits and return it as bitslength = 1 buffer[0] = GetSingleBitAt(c,nOffsetInBits); *bitsreturned = 1; } else // Ordinary Grey Scale Images or special value image (DEM) { // Should I consider the Indians?? // All the data will be stored in BigEndian format. // Directly return data int nIndexInByte = (nTIFFBitsPerPixel* (cols * row + col )) / 8; int nOffsetInBits = (nTIFFBitsPerPixel* (cols * row + col )) % 8; if(bps < 8) { byte c = ((byte*)pDataBuffer)[nIndexInByte]; if(fillorder != FILLORDER_MSB2LSB) { // We need to SWAP bits -- ABCDEFGH becomes HGFEDCBA byte tempbyte = 0; if(c & 128) tempbyte += 1; if(c & 64) tempbyte += 2; if(c & 32) tempbyte += 4; if(c & 16) tempbyte += 8; if(c & 8) tempbyte += 16; if(c & 4) tempbyte += 32; if(c & 2) tempbyte += 64; if(c & 1) tempbyte += 128; c = tempbyte; } //less than 1 bytes buffer = (byte*)malloc(sizeof(char)); if(bps == 1) { buffer[0] = GetSingleBitAt(c, nOffsetInBits); } else if(bps == 2 || bps == 4) { buffer[0] = GetMultiBitAt(c, nOffsetInBits,bps); } else { if(outmode != MODE_QUIET) fprintf(stderr, "**Error**: Unsupported PLANARCONFIG_CONTIG sample format.\n"); exit_nicely(conn,1); } } else if(bps == 8 && format == SAMPLEFORMAT_INT) //char { *bitsreturned = 8; buffer = (byte*)malloc(1*sizeof(char)); buffer[0]= ((char*)pDataBuffer)[nIndex]; } else if(bps == 8) //unsigned char / byte { *bitsreturned = 8; buffer = (byte*)malloc(1*sizeof(char)); buffer[0]= ((byte*)pDataBuffer)[nIndex]; } else if(bps == 16 && format == SAMPLEFORMAT_INT) // short { int16 c; *bitsreturned = 2*8; buffer = (byte*)malloc(2*sizeof(char)); c= ((int16*)pDataBuffer)[nIndex]; memcpy(buffer,(byte*)&c,2); if(tiff_bigendian) *((int16*)buffer) = _ntohs(*((int16*)buffer)); } else if(bps == 16) // unsigned short { uint16 c; *bitsreturned = 2*8; buffer = (byte*)malloc(2*sizeof(char)); c= ((uint16*)(pDataBuffer+nIndex))[band]; memcpy(buffer,(byte*)&c,2); if(tiff_bigendian) *((uint16*)buffer) = _ntohs(*((uint16*)buffer)); } else if(bps == 32 && format == SAMPLEFORMAT_INT) // long { int32 c; *bitsreturned = 4*8; buffer = (byte*)malloc(4*sizeof(char)); c= ((int32*)(pDataBuffer+nIndex))[band]; memcpy(buffer,(byte*)&c,4); if(tiff_bigendian) *((int32*)buffer) = _ntohl(*((int32*)buffer)); } else if(bps == 32 && format == SAMPLEFORMAT_IEEEFP) // IEEE FP { float c; *bitsreturned = 4*8; buffer = (byte*)malloc(sizeof(float)); c= ((float*)(pDataBuffer+nIndex))[band]; memcpy(buffer,(byte*)&c,4); if(tiff_bigendian) *((float*)buffer) = _ntohf(*((float*)buffer)); } else if(bps == 32) // unsigned long { uint32 c; *bitsreturned = 4*8; buffer = (byte*)malloc(4*sizeof(char)); c= ((uint32*)(pDataBuffer+nIndex))[band]; memcpy(buffer,(byte*)&c,4); if(tiff_bigendian) *((uint32*)buffer) = _ntohl(*((uint32*)buffer)); } else { if(outmode != MODE_QUIET) fprintf(stderr, "Invalid PHOTOMETRIC Mode or Sample Format in single channel TIFF image (mode=%d)", mode); exit_nicely(conn,1); } } } return buffer; } void _EndReadTIFF() { // Free memory usage if(pDataBuffer != NULL) free(pDataBuffer); pDataBuffer = NULL; } /* Step 3. CreateIndex GiST Index on MBR of each block */ void CreateIndex() { char query[QUERY_BUF_LENGTH]; memset(query, '\0', QUERY_BUF_LENGTH); if(outmode == MODE_VERBOSE) fprintf(stdout, "Start building GiST on block mbr.\n"); sprintf(query,"CREATE INDEX \"index_%s_%s\" ON \"%s\".\"%s\" using GIST ( blockMBR GIST_GEOMETRY_OPS)", schema, table, schema, table); res = PQexecParams(conn, query, 0, NULL, NULL, NULL, NULL, 0); if ( ! res || PQresultStatus(res) != PGRES_COMMAND_OK ) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error** when create spatial index on block MBR at server side: %s\n.",PQerrorMessage(conn)); exit_nicely(conn, 1); } if(outmode == MODE_VERBOSE) fprintf(stdout, "Finish GiST on block mbr.\n"); } /* Step 4. LoadSRS */ void LoadSRS() { // Update record in PGRASTER_SRS table double adfCoeff[6], x, y; char query[QUERY_BUF_LENGTH]; memset(query,'\0',QUERY_BUF_LENGTH); if(outmode == MODE_VERBOSE) fprintf(stdout, "Start loading SRS information.\n"); /* * Compute the coefficients. */ x = 0.5; y = 0.5; if( !GTIFImageToPCS( gtif, &x, &y ) ) return; adfCoeff[4] = x; adfCoeff[5] = y; x = 1.5; y = 0.5; if( !GTIFImageToPCS( gtif, &x, &y ) ) return; adfCoeff[0] = x - adfCoeff[4]; adfCoeff[1] = y - adfCoeff[5]; x = 0.5; y = 1.5; if( !GTIFImageToPCS( gtif, &x, &y ) ) return; adfCoeff[2] = x - adfCoeff[4]; adfCoeff[3] = y - adfCoeff[5]; //Insert into PostgreSQL // adfCoeff[0-5] = [A-F] // A(0), B(1), C(2), D(3), E(4), F(5) // Xgeo = E + Xpixel * A + Ypixel * C // Ygeo = F + Ypixel * D + Xpixel * B sprintf(query, "UPDATE \"%s\".\"pgraster_srs\" SET rowdenominator='{%lf,%lf,%lf}', columndenominator='{%lf,%lf,%lf}' WHERE rasterObjectID=%Ld",schema, adfCoeff[5], adfCoeff[3],adfCoeff[1],adfCoeff[4], adfCoeff[0],adfCoeff[2],metadata.rasterObjectID); res = PQexec(conn,query); if ( ! res || PQresultStatus(res) != PGRES_COMMAND_OK ) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error** when update SRS information at server side: %s\n.",PQerrorMessage(conn)); exit_nicely(conn, 1); } PQclear(res); if(outmode == MODE_VERBOSE) fprintf(stdout, "Finish loading SRS information.\n"); } /* Step 5. UpdateMetadata */ void UpdateMetadata() { // Some metadata need to be updated after the data loading. char query[QUERY_BUF_LENGTH]; memset(query, '\0' , QUERY_BUF_LENGTH); if(outmode == MODE_VERBOSE) fprintf(stdout, "Start updating metadata second time.\n"); sprintf(query,"UPDATE \"%s\".\"pgraster_metadata\" SET rasterpyramiddepth=%d WHERE rasterObjectID=%Ld", schema, (int)nPyramidDepth, metadata.rasterObjectID); res = PQexecParams(conn, query, 0, NULL, NULL, NULL, NULL, 0); if ( ! res || PQresultStatus(res) != PGRES_COMMAND_OK ) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error** when updating spatial metadata at server side: %s\n.",PQerrorMessage(conn)); exit_nicely(conn, 1); } memset(query, '\0' , QUERY_BUF_LENGTH); sprintf(query,"UPDATE \"%s\".\"pgraster_metadata\" SET rastercelldepth=%d WHERE rasterObjectID=%Ld", schema, (int)metadata.rasterCellDepth, metadata.rasterObjectID); res = PQexecParams(conn, query, 0, NULL, NULL, NULL, NULL, 0); if ( ! res || PQresultStatus(res) != PGRES_COMMAND_OK ) { if(outmode != MODE_QUIET) fprintf(stderr, "**Error** when updating spatial metadata at server side: %s\n.",PQerrorMessage(conn)); exit_nicely(conn, 1); } if(outmode == MODE_VERBOSE) fprintf(stdout, "Finish updateing metadata second time.\n"); } /* Step 6. Cleanup */ void Cleanup() { /* end the transaction */ res = PQexec(conn, "END"); PQclear(res); beginTransaction = FALSE; /* Close Connection*/ PQfinish(conn); /* Close GeoTIFF File */ if (gtif) GTIFFree(gtif); if (tif) XTIFFClose(tif); #ifdef DEBUG if(debug != NULL) fclose(debug); #endif } /* Step 7. PrintResultMessage */ void PrintResultMessage(int sucessful) { if(outmode == MODE_QUIET) return; if(sucessful == 1) { //sucessful fprintf(stdout, "GeoTIFF file '%s' succesfully imported.\n", geotiff_file); } else { //fail fprintf(stdout, "GeoTIFF file '%s' importing failed.\n", geotiff_file); } } #ifdef USE_ICONV char * utf8 (const char *fromcode, char *inputbuf) { iconv_t cd; char *outputptr; char *outputbuf; size_t outbytesleft; size_t inbytesleft; inbytesleft = strlen (inputbuf); cd = iconv_open ("UTF-8", fromcode); if (cd == (iconv_t) - 1) { fprintf(stderr, "utf8: iconv_open: %s\n", strerror (errno)); return NULL; } outbytesleft = inbytesleft*3+1; /* UTF8 string can be 3 times larger */ /* then local string */ outputbuf = (char *) malloc (outbytesleft); if (!outputbuf) { fprintf(stderr, "utf8: malloc: %s\n", strerror (errno)); return NULL; } memset (outputbuf, 0, outbytesleft); outputptr = outputbuf; if (-1==iconv(cd, &inputbuf, &inbytesleft, &outputptr, &outbytesleft)) { fprintf(stderr, "utf8: %s", strerror (errno)); return NULL; } iconv_close (cd); return outputbuf; } #endif /* defined USE_ICONV */