1 : /******************************************************************************
2 : * File : wktrasterdataset.cpp
3 : * Project: WKT Raster driver
4 : * Purpose: GDAL Dataset code for WKTRaster driver
5 : * Author: Jorge Arevalo, jorgearevalo@gis4free.org
6 : *
7 : * Last changes:
8 : * $Id: wktrasterdataset.cpp 18325 2009-12-17 17:29:51Z jorgearevalo $
9 : *
10 : * NOTE:
11 : * There are many small functions in this class. Some of these functions
12 : * perform database queries. From the point of view of performance, this
13 : * may be a bad idea, but from point of view of code manintenance, I think
14 : * is better.
15 : * Anyway, the WKT Raster project is being developed at same time that this
16 : * driver,and I think that is better to allow fast changes in the code
17 : * (easy to understand and to change) than performance, at this moment
18 : * (August 2009)
19 : *
20 : * TODO:
21 : * - Allow non-regular blocking table arrangements (in general)
22 : * - Allow PQ connection string parameters in any order
23 : * - Update data blocks in SetProjection
24 : * - Update data blocks in SetGeoTransform
25 : * - Disable sequential scanning in OpenConnection in the database instance
26 : * has support for GEOMETRY type (is a good idea?)
27 : * - In SetRasterProperties, if we can't create OGRGeometry from database,
28 : * we should "attack" directly the WKT representation of geometry, not
29 : * abort the program. For example, when pszProjectionRef is NULL
30 : * (SRID = -1), when poSR or prGeom are NULL...
31 : *
32 : * Low priority:
33 : * - Add support for rotated images in SetRasterProperties
34 : * - Check if the tiles of a table REALLY have regular blocking arrangement.
35 : * The fact that "regular_blocking" field in RASTER_COLUMNS table is set
36 : * to TRUE doesn't proof it. Method TableHasRegularBlocking.
37 : *
38 : ******************************************************************************
39 : * Copyright (c) 2009, Jorge Arevalo, jorgearevalo@gis4free.org
40 : *
41 : * Permission is hereby granted, free of charge, to any person obtaining a
42 : * copy of this software and associated documentation files (the "Software"),
43 : * to deal in the Software without restriction, including without limitation
44 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
45 : * and/or sell copies of the Software, and to permit persons to whom the
46 : * Software is furnished to do so, subject to the following conditions:
47 : *
48 : * The above copyright notice and this permission notice shall be included
49 : * in all copies or substantial portions of the Software.
50 : *
51 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
52 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
53 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
54 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
55 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
56 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
57 : * DEALINGS IN THE SOFTWARE.
58 : ******************************************************************************/
59 :
60 : #include "wktraster.h"
61 : #include <stdlib.h>
62 : #include "ogr_api.h"
63 : #include "ogr_geometry.h"
64 : #include "gdal.h"
65 : #include "cpl_conv.h"
66 : #include "cpl_string.h"
67 : #include "gdal_priv.h"
68 : #include <math.h>
69 : #include <cpl_error.h>
70 :
71 : #ifdef _WIN32
72 : #define rint(x) floor((x) + 0.5)
73 : #endif
74 :
75 : CPL_CVSID("$Id: wktrasterdataset.cpp 18325 2009-12-17 17:29:51Z jorgearevalo $");
76 :
77 : CPL_C_START
78 : void GDALRegister_WKTRaster(void);
79 :
80 : CPL_C_END
81 :
82 : /**
83 : * Constructor. Init the class properties with default values
84 : */
85 0 : WKTRasterDataset::WKTRasterDataset() {
86 0 : hPGconn = NULL;
87 0 : bCloseConnection = FALSE;
88 0 : pszSchemaName = NULL;
89 0 : pszTableName = NULL;
90 0 : pszRasterColumnName = NULL;
91 0 : pszWhereClause = NULL;
92 0 : pszProjection = NULL;
93 0 : pszWorkingMode = NULL;
94 0 : bTableHasGISTIndex = FALSE;
95 0 : nVersion = 0;
96 0 : nBlockSizeX = 0;
97 0 : nBlockSizeY = 0;
98 0 : dfPixelSizeX = 0.0;
99 0 : dfPixelSizeY = 0.0;
100 0 : dfUpperLeftX = 0.0;
101 0 : dfUpperLeftY = 0.0;
102 0 : dfLowerRightX = 0.0;
103 0 : dfLowerRightY = 0.0;
104 0 : dfRotationX = 0.0;
105 0 : dfRotationY = 0.0;
106 0 : nSrid = -1;
107 0 : nOverviews = 0;
108 0 : papoWKTRasterOv = NULL;
109 0 : poOutdbRasterDS = NULL;
110 :
111 : // TO BE DELETED
112 0 : papoBlocks = NULL;
113 0 : nBlocks = 0;
114 :
115 0 : }
116 :
117 : /**
118 : * Destructor. Frees allocated memory.
119 : */
120 0 : WKTRasterDataset::~WKTRasterDataset() {
121 0 : CPLFree(pszSchemaName);
122 0 : CPLFree(pszTableName);
123 0 : CPLFree(pszWorkingMode);
124 0 : CPLFree(pszWhereClause);
125 0 : CPLFree(pszRasterColumnName);
126 :
127 0 : if (nOverviews > 0) {
128 0 : for (int i = 0; i < nOverviews; i++) {
129 0 : delete papoWKTRasterOv[i];
130 0 : papoWKTRasterOv[i] = NULL;
131 : }
132 :
133 0 : CPLFree(papoWKTRasterOv);
134 : }
135 :
136 :
137 0 : if (nBlocks > 0) {
138 0 : for(int i = 0; i < nBlocks; i++) {
139 0 : delete papoBlocks[i];
140 0 : papoBlocks[i] = NULL;
141 : }
142 :
143 0 : CPLFree(papoBlocks);
144 : }
145 :
146 : /**
147 : * bCloseConnection = TRUE for Dataset, FALSE for Overviews
148 : */
149 0 : if (bCloseConnection == TRUE && hPGconn != NULL) {
150 0 : PQfinish(hPGconn);
151 0 : hPGconn = NULL;
152 : }
153 :
154 0 : if (poOutdbRasterDS != NULL) {
155 0 : GDALClose((GDALDatasetH)poOutdbRasterDS);
156 : }
157 0 : }
158 :
159 : /**
160 : * Check if the table has an index
161 : * Parameters:
162 : * - PGconn *: pointer to connection
163 : * - const char *: table name
164 : * - const char *: table schema
165 : * Output:
166 : * - GBool: TRUE if the table has an index, FALSE otherwise
167 : */
168 : static
169 0 : GBool TableHasGISTIndex(PGconn * hPGconn, const char * pszTable,
170 : const char * pszSchema) {
171 :
172 : // Security checkings
173 0 : VALIDATE_POINTER1(hPGconn, "TableHasGISTIndex", FALSE);
174 0 : VALIDATE_POINTER1(pszTable, "TableHasGISTIndex", FALSE);
175 :
176 :
177 0 : CPLString osCommand;
178 0 : PGresult * hPGresult = NULL;
179 : char * pszResultString;
180 0 : GBool bTableHasIndex = FALSE;
181 :
182 : /*****************************************************************
183 : * Check if the table has an index in PostgreSQL system tables
184 : *****************************************************************/
185 : osCommand.Printf(
186 : "SELECT relhasindex \
187 : FROM pg_class, pg_attribute, pg_type, pg_namespace \
188 : WHERE \
189 : pg_namespace.nspname = '%s' and \
190 : pg_namespace.oid = pg_class.relnamespace and \
191 : pg_class.relname = '%s' and \
192 : pg_class.oid = pg_attribute.attrelid and \
193 : pg_attribute.atttypid = pg_type.oid and \
194 : pg_type.typname = 'raster'",
195 0 : pszTable, pszSchema);
196 :
197 0 : hPGresult = PQexec(hPGconn, osCommand.c_str());
198 0 : if (
199 : hPGresult != NULL &&
200 : PQresultStatus(hPGresult) == PGRES_TUPLES_OK &&
201 : PQntuples(hPGresult) > 0
202 : ) {
203 0 : pszResultString = PQgetvalue(hPGresult, 0, 0);
204 0 : if (pszResultString != NULL && EQUALN(pszResultString, "t",
205 : sizeof (char)) != 0) {
206 0 : bTableHasIndex = TRUE;
207 0 : PQclear(hPGresult);
208 : } else {
209 0 : bTableHasIndex = FALSE;
210 0 : PQclear(hPGresult);
211 : }
212 :
213 : } else {
214 0 : if (hPGresult != NULL)
215 0 : PQclear(hPGresult);
216 0 : bTableHasIndex = FALSE;
217 : }
218 :
219 0 : return bTableHasIndex;
220 :
221 : }
222 :
223 :
224 : /**
225 : * Check if the RASTER_COLUMNS table exists
226 : * Parameters:
227 : * - PGconn *: pointer to connection
228 : * Output:
229 : * - GBool: TRUE if the table RASTER_COLUMNS does exist, FALSE otherwise
230 : */
231 : static
232 0 : GBool ExistsRasterColumnsTable(PGconn * hPGconn) {
233 : // Security checkings
234 0 : VALIDATE_POINTER1(hPGconn, "ExistsRasterColumnsTable", FALSE);
235 :
236 0 : const char* szCommand = "SELECT relname FROM pg_class \
237 : WHERE relname = 'raster_columns'";
238 0 : PGresult * hPGresult = NULL;
239 : char * pszResultString;
240 0 : GBool bExistsRasterColumnsTable = FALSE;
241 :
242 : /************************************************************************
243 : * Look for RASTER_COLUMNS table
244 : ************************************************************************/
245 :
246 0 : hPGresult = PQexec(hPGconn, szCommand);
247 0 : if (
248 : hPGresult != NULL &&
249 : PQresultStatus(hPGresult) == PGRES_TUPLES_OK &&
250 : PQntuples(hPGresult) > 0
251 : ) {
252 0 : pszResultString = PQgetvalue(hPGresult, 0, 0);
253 0 : if (pszResultString != NULL &&
254 : EQUALN(pszResultString,
255 : "raster_columns", 14 * sizeof (char)) != 0) {
256 0 : bExistsRasterColumnsTable = TRUE;
257 0 : PQclear(hPGresult);
258 : } else {
259 0 : bExistsRasterColumnsTable = FALSE;
260 0 : PQclear(hPGresult);
261 : }
262 :
263 : } else {
264 0 : if (hPGresult != NULL)
265 0 : PQclear(hPGresult);
266 0 : bExistsRasterColumnsTable = FALSE;
267 : }
268 :
269 0 : return bExistsRasterColumnsTable;
270 : }
271 :
272 :
273 : /**
274 : * Check if exists an overviews table
275 : * Parameters:
276 : * - PGconn *: pointer to connection
277 : * Output:
278 : * - GBool: TRUE if the table RASTER_OVERVIEWS does exist, FALSE otherwise
279 : */
280 : static
281 0 : GBool ExistsOverviewsTable(PGconn * hPGconn) {
282 :
283 : // Security checkings
284 0 : VALIDATE_POINTER1(hPGconn, "ExistsOverviewsTable", FALSE);
285 :
286 : const char* szCommand =
287 0 : "SELECT relname FROM pg_class WHERE relname = 'raster_overviews'";
288 0 : PGresult * hPGresult = NULL;
289 : char * pszResultString;
290 0 : GBool bExistsOverviewTable = FALSE;
291 :
292 : /************************************************************************
293 : * Look for RASTER_OVERVIEWS table
294 : ************************************************************************/
295 :
296 0 : hPGresult = PQexec(hPGconn, szCommand);
297 0 : if (
298 : hPGresult != NULL &&
299 : PQresultStatus(hPGresult) == PGRES_TUPLES_OK &&
300 : PQntuples(hPGresult) > 0
301 : ) {
302 0 : pszResultString = PQgetvalue(hPGresult, 0, 0);
303 0 : if (pszResultString != NULL &&
304 : EQUALN(pszResultString,
305 : "raster_overviews", 16 * sizeof (char)) != 0) {
306 0 : bExistsOverviewTable = TRUE;
307 0 : PQclear(hPGresult);
308 : } else {
309 0 : bExistsOverviewTable = FALSE;
310 0 : PQclear(hPGresult);
311 : }
312 :
313 : } else {
314 0 : if (hPGresult != NULL)
315 0 : PQclear(hPGresult);
316 0 : bExistsOverviewTable = FALSE;
317 : }
318 :
319 0 : return bExistsOverviewTable;
320 :
321 : }
322 :
323 : /**
324 : * Check if the table has regular_blocking constraint.
325 : * NOTE: From now (July 2009), only tables with regular_blocking
326 : * are allowed
327 : * Parameters:
328 : * - PGconn *: pointer to connection
329 : * - const char *: table name
330 : * - const char *: raster column name
331 : * - const char *: schema name
332 : * Returns:
333 : * - GBool: TRUE if the table has regular_blocking constraint, FALSE otherwise
334 : */
335 : static
336 0 : GBool TableHasRegularBlocking(PGconn * hPGconn, const char * pszTable,
337 : const char * pszColumn, const char * pszSchema) {
338 :
339 : // Security checkings
340 0 : VALIDATE_POINTER1(hPGconn, "TableHasRegularBlocking", FALSE);
341 0 : VALIDATE_POINTER1(pszTable, "TableHasRegularBlocking", FALSE);
342 0 : VALIDATE_POINTER1(pszColumn, "TableHasRegularBlocking", FALSE);
343 :
344 0 : CPLString osCommand;
345 0 : PGresult * hPGresult = NULL;
346 : char * pszResultString;
347 0 : GBool bTableHasRegularBlocking = FALSE;
348 :
349 :
350 : /************************************************************************
351 : * Look for regular_blocking in RASTER_COLUMNS table
352 : ************************************************************************/
353 : osCommand.Printf(
354 : "SELECT \
355 : regular_blocking \
356 : FROM raster_columns \
357 : WHERE \
358 : r_table_name = '%s' and \
359 : r_column = '%s' and \
360 : r_table_schema = '%s'",
361 0 : pszTable, pszColumn, pszSchema);
362 :
363 :
364 :
365 0 : hPGresult = PQexec(hPGconn, osCommand.c_str());
366 0 : if (
367 : hPGresult != NULL &&
368 : PQresultStatus(hPGresult) == PGRES_TUPLES_OK &&
369 : PQntuples(hPGresult) > 0
370 : ) {
371 0 : pszResultString = PQgetvalue(hPGresult, 0, 0);
372 0 : if (pszResultString != NULL && EQUALN(pszResultString, "t", sizeof (char)) != 0) {
373 0 : bTableHasRegularBlocking = TRUE;
374 0 : PQclear(hPGresult);
375 : } else {
376 0 : bTableHasRegularBlocking = FALSE;
377 0 : PQclear(hPGresult);
378 : }
379 :
380 : } else {
381 0 : if (hPGresult != NULL)
382 0 : PQclear(hPGresult);
383 0 : bTableHasRegularBlocking = FALSE;
384 : }
385 :
386 0 : return bTableHasRegularBlocking;
387 :
388 : }
389 :
390 : /**
391 : * Get the name of the column of type raster
392 : * Parameters:
393 : * PGconn *: pointer to connection
394 : * const char *: schema name
395 : * const char *: table name
396 : * Output:
397 : * char *: The column name, or NULL if doesn't exist
398 : */
399 : static
400 0 : char * GetWKTRasterColumnName(PGconn * hPGconn, const char * pszSchemaName,
401 : const char * pszTable) {
402 :
403 : // Security checkings
404 0 : VALIDATE_POINTER1(hPGconn, "GetWKTRasterColumnName", NULL);
405 0 : VALIDATE_POINTER1(pszTable, "GetWKTRasterColumnName", NULL);
406 :
407 0 : CPLString osCommand;
408 0 : PGresult * hPGresult = NULL;
409 : char * pszResultString;
410 :
411 :
412 : /************************************************************
413 : * Get the attribute name of type 'raster' of the table
414 : ************************************************************/
415 : osCommand.Printf(
416 : "SELECT attname \
417 : FROM pg_class, pg_attribute, pg_type, pg_namespace \
418 : WHERE \
419 : pg_namespace.nspname = '%s' and \
420 : pg_namespace.oid = pg_class.relnamespace and \
421 : pg_class.relname = '%s' and \
422 : pg_class.oid = pg_attribute.attrelid and \
423 : pg_attribute.atttypid = pg_type.oid and \
424 : pg_type.typname = 'raster'",
425 0 : pszSchemaName, pszTable);
426 :
427 0 : hPGresult = PQexec(hPGconn, osCommand.c_str());
428 0 : if (
429 : hPGresult != NULL &&
430 : PQresultStatus(hPGresult) == PGRES_TUPLES_OK &&
431 : PQntuples(hPGresult) > 0
432 : ) {
433 0 : pszResultString = CPLStrdup(PQgetvalue(hPGresult, 0, 0));
434 0 : PQclear(hPGresult);
435 :
436 : } else {
437 0 : if (hPGresult != NULL)
438 0 : PQclear(hPGresult);
439 0 : pszResultString = NULL;
440 : }
441 :
442 :
443 0 : return pszResultString;
444 :
445 : }
446 :
447 : /**
448 : * Open a database connection and perform some security checkings
449 : * Parameters:
450 : * const char *: Connection string
451 : * Output:
452 : * PGconn *: pointer to object connection, or NULL it there was a problem
453 : */
454 : static
455 1 : PGconn * OpenConnection(const char *pszConnectionString) {
456 :
457 : // Security checkings
458 1 : VALIDATE_POINTER1(pszConnectionString, "OpenConnection", NULL);
459 :
460 1 : PGconn * hPGconn = NULL;
461 1 : PGresult * hPGresult = NULL;
462 :
463 :
464 : /********************************************************
465 : * Connect with database
466 : ********************************************************/
467 1 : hPGconn = PQconnectdb(pszConnectionString + 3);
468 1 : if (hPGconn == NULL || PQstatus(hPGconn) == CONNECTION_BAD) {
469 : CPLError(CE_Failure, CPLE_AppDefined, "PGconnectcb failed.\n%s",
470 1 : PQerrorMessage(hPGconn));
471 1 : PQfinish(hPGconn);
472 1 : return NULL;
473 : }
474 :
475 : /*****************************************************************
476 : * Test to see if this database instance has support for the
477 : * PostGIS Geometry type.
478 : * TODO: If so, disable sequential scanning so we will get the
479 : * value of the gist indexes. (is it a good idea?)
480 : *****************************************************************/
481 : hPGresult = PQexec(hPGconn,
482 0 : "SELECT oid FROM pg_type WHERE typname = 'geometry'");
483 0 : if (
484 : hPGresult == NULL ||
485 : PQresultStatus(hPGresult) != PGRES_TUPLES_OK ||
486 : PQntuples(hPGresult) <= 0
487 : ) {
488 : CPLError(CE_Failure, CPLE_AppDefined,
489 0 : "Can't find geometry type, is Postgis correctly installed ?\n");
490 0 : if (hPGresult)
491 0 : PQclear(hPGresult);
492 0 : PQfinish(hPGconn);
493 0 : hPGconn = NULL;
494 :
495 : }
496 :
497 0 : return hPGconn;
498 : }
499 :
500 : /**
501 : * Extract a field from the connection string:
502 : * PG:[host='<host>'] [user='<user>'] [password='<password>']
503 : * dbname='<dbname>' table='<raster_table>' [mode='<working_mode>'
504 : * where='<where_clause>']
505 : * The idea is to extract the fields that PQconnect function doesn't accept
506 : * (this is, 'table', 'mode' and 'where')
507 : *
508 : * Parameters:
509 : * char **: pointer to connection string, to be modified
510 : * char *: field init (ex.: "where=", "table=")
511 : * Output:
512 : * char *: The field (where clause or table name) of the connection string
513 : * if everything works fine, NULL otherwise
514 : */
515 : static
516 4 : char * ExtractField(char ** ppszConnectionString, const char * pszFieldInit) {
517 :
518 : // Security checkings
519 4 : VALIDATE_POINTER1(ppszConnectionString, "ExtractField", NULL);
520 4 : VALIDATE_POINTER1(*ppszConnectionString, "ExtractField", NULL);
521 4 : VALIDATE_POINTER1(pszFieldInit, "ExtractField", NULL);
522 :
523 4 : char * pszField = NULL;
524 4 : char * pszStart = NULL;
525 :
526 :
527 : /************************************************************************
528 : * Get the value of the field to extract, and delete the whole field
529 : * from the connection string
530 : ************************************************************************/
531 4 : int nFieldInitLen = strlen(pszFieldInit);
532 :
533 : /*
534 : * Get the initial position of the field in the string
535 : */
536 4 : pszStart = strstr(*ppszConnectionString, pszFieldInit);
537 4 : if (pszStart == NULL)
538 2 : return NULL;
539 :
540 2 : int bHasQuote = pszStart[nFieldInitLen] == '\'';
541 :
542 : // Copy the field of the connection string to another var
543 2 : pszField = CPLStrdup(pszStart + nFieldInitLen + bHasQuote);
544 2 : char* pszEndField = strchr(pszField, (bHasQuote) ? '\'' : ' ');
545 2 : if (pszEndField)
546 : {
547 2 : *pszEndField = '\0';
548 : char* pszEnd = pszStart + nFieldInitLen + bHasQuote +
549 2 : (pszEndField - pszField) + 1;
550 2 : memmove(pszStart, pszEnd, strlen(pszEnd) + 1);
551 : }
552 : else
553 : // Delete the field's part from the connection string
554 0 : *pszStart = '\0';
555 :
556 2 : return pszField;
557 : }
558 :
559 : /**
560 : * Populate the georeference information fields of dataset
561 : * Input:
562 : * Output:
563 : * CPLErr: CE_None if the fields are populated, CE_Failure otherwise
564 : */
565 0 : CPLErr WKTRasterDataset::SetRasterProperties() {
566 0 : CPLString osCommand;
567 0 : PGresult * hPGresult = NULL;
568 0 : char * pszWKTRaster = NULL;
569 0 : OGRSpatialReference * poSR = NULL;
570 0 : OGRGeometry * poGeom = NULL;
571 0 : OGRErr OgrErr = OGRERR_NONE;
572 0 : OGREnvelope * poE = NULL;
573 0 : char * pszProjectionRef = NULL;
574 :
575 :
576 : /********************************************************************
577 : * Check working mode. This is the regular_blocking mode
578 : ********************************************************************/
579 0 : if (EQUALN(pszWorkingMode,
580 : REGULARLY_TILED_MODE, strlen(REGULARLY_TILED_MODE))) {
581 : /********************************************
582 : * Populate georeference related fields
583 : * Note: From now (July 2009), only
584 : * north up images (unrotated) will be
585 : * supported.
586 : *
587 : * TODO: Support for rotated images?
588 : ********************************************/
589 0 : dfRotationX = 0.0;
590 0 : dfRotationY = 0.0;
591 :
592 : /**************************************************************
593 : * Try to fetch info from RASTER_COLUMNS table
594 : **************************************************************/
595 : osCommand.Printf(
596 : "SELECT srid, ST_AsText(extent), pixelsize_x, pixelsize_y, \
597 : blocksize_x, blocksize_y FROM raster_columns WHERE r_table_name\
598 : = '%s' AND r_column = '%s' and r_table_schema = '%s'",
599 0 : pszTableName, pszRasterColumnName, pszSchemaName);
600 :
601 0 : hPGresult = PQexec(hPGconn, osCommand.c_str());
602 0 : if (
603 : hPGresult != NULL &&
604 : PQresultStatus(hPGresult) == PGRES_TUPLES_OK &&
605 : PQntuples(hPGresult) > 0
606 : ) {
607 :
608 0 : nSrid = atoi(PQgetvalue(hPGresult, 0, 0));
609 :
610 0 : pszWKTRaster = CPLStrdup(PQgetvalue(hPGresult, 0, 1));
611 :
612 0 : dfPixelSizeX = atof(PQgetvalue(hPGresult, 0, 2));
613 0 : dfPixelSizeY = atof(PQgetvalue(hPGresult, 0, 3));
614 0 : nBlockSizeX = atoi(PQgetvalue(hPGresult, 0, 4));
615 0 : nBlockSizeY = atoi(PQgetvalue(hPGresult, 0, 5));
616 :
617 0 : PQclear(hPGresult);
618 :
619 : /******************************************************************
620 : * TODO If we can't create OGRGeometry from database, we should
621 : * "attack" directly the WKT representation of geometry, not abort
622 : * the program. For example, when pszProjectionRef is NULL
623 : * (SRID = -1), when poSR or prGeom are NULL...
624 : ******************************************************************/
625 0 : pszProjectionRef = (char *) GetProjectionRef();
626 0 : poSR = new OGRSpatialReference(pszProjectionRef);
627 : OgrErr = OGRGeometryFactory::createFromWkt(&pszWKTRaster, poSR,
628 0 : &poGeom);
629 0 : if (OgrErr != OGRERR_NONE) {
630 :
631 : CPLError(CE_Failure, CPLE_AppDefined,
632 0 : "Can't get WKT Raster extent from database\n");
633 0 : CPLFree(pszWKTRaster);
634 :
635 : /**
636 : * NOTE: Projection ref can be an empty string. Should it be
637 : * free in that case?
638 : */
639 0 : CPLFree(pszProjectionRef);
640 :
641 0 : return CE_Failure;
642 : }
643 :
644 :
645 : /****************************************************************
646 : * Get upper left and lower rights coordinates for this block.
647 : * The lower right coords are not really needed...
648 : ****************************************************************/
649 0 : poE = new OGREnvelope();
650 0 : poGeom->getEnvelope(poE);
651 0 : dfUpperLeftX = poE->MinX;
652 0 : dfUpperLeftY = poE->MaxY;
653 0 : dfLowerRightX = poE->MaxX;
654 0 : dfLowerRightY = poE->MinY;
655 :
656 : /****************************************************************
657 : * Get raster size from extent (Geometry)
658 : ****************************************************************/
659 : nRasterXSize = (int)
660 0 : fabs(rint((poE->MaxX - poE->MinX) / dfPixelSizeX));
661 : nRasterYSize = (int)
662 0 : fabs(rint((poE->MaxY - poE->MinY) / dfPixelSizeY));
663 :
664 :
665 0 : return CE_None;
666 : }
667 :
668 : /**
669 : * If the table isn't listed in RASTER_COLUMNS then it isn't
670 : * regularly_blocked for our purposes. This is only temporarily
671 : */
672 : else {
673 : CPLError(CE_Failure, CPLE_AppDefined,
674 0 : "Can't get georeference coordinates from database\n");
675 0 : return CE_Failure;
676 : }
677 :
678 : }
679 :
680 :
681 : else {
682 : // Here other working modes will be check
683 0 : return CE_Failure;
684 0 : }
685 : }
686 :
687 : /**
688 : * Explode a string representing an array into an array of strings. The input
689 : * string has this format: {element1,element2,element3,....,elementN}
690 : * Parameters:
691 : * - const char *: a string representing an array
692 : * - int *: pointer to an int that will contain the number of elements
693 : * Returns:
694 : * char **: An array of strings, one per element. Must be freed with CSLDestroy
695 : */
696 0 : char ** WKTRasterDataset::ExplodeArrayString(
697 : const char * pszPQarray, int * pnNumberOfElements) {
698 :
699 : // integrity checkings
700 0 : if (
701 : pszPQarray == NULL ||
702 0 : pszPQarray[0] != '{' ||
703 0 : pszPQarray[strlen(pszPQarray) - 1] != '}'
704 : ) {
705 0 : if (pnNumberOfElements)
706 0 : *pnNumberOfElements = 0;
707 0 : return NULL;
708 : }
709 :
710 0 : char* pszTemp = CPLStrdup(pszPQarray + 1);
711 0 : pszTemp[strlen(pszTemp) - 1] ='\0';
712 0 : char** papszRet = CSLTokenizeString2( pszTemp, ",", 0 );
713 0 : CPLFree(pszTemp);
714 :
715 0 : if (pnNumberOfElements)
716 0 : *pnNumberOfElements = CSLCount(papszRet);
717 0 : return papszRet;
718 : }
719 :
720 :
721 : /**
722 : * Implode an array of strings, to transform it into a PostgreSQL-style
723 : * array, with this format: {element1,element2,element3,....,elementN}
724 : * Input:
725 : * char **: An array of strings, one per element
726 : * int: An int that contains the number of elements
727 : * Output:
728 : * const char *: a string representing an array
729 : */
730 0 : char * WKTRasterDataset::ImplodeStrings(char ** papszElements, int nElements)
731 : {
732 0 : VALIDATE_POINTER1(papszElements, "ImplodeStrings", NULL);
733 :
734 :
735 0 : char * pszPQarray = NULL;
736 : char * ptrString;
737 : char szTemp[1024];
738 0 : unsigned int nCharsCopied = 0;
739 0 : unsigned int nNumberOfCommas = 0;
740 0 : unsigned int nNumAvailableBytes = 0;
741 0 : int nPosOpeningBracket = 0;
742 0 : int nPosClosingBracket = 0;
743 :
744 : /**************************************************************************
745 : * Array for the string. We could allocate memory for all the elements of
746 : * the array, plus commas and brackets, but 1024 bytes should be enough to
747 : * store small values, and the method is faster in this way
748 : **************************************************************************/
749 0 : pszPQarray = (char *)CPLCalloc(1024, sizeof(char));
750 0 : VALIDATE_POINTER1(pszPQarray, "ImplodeStrings", NULL);
751 :
752 :
753 : // empty string
754 0 : if (nElements <= 0) {
755 0 : nPosOpeningBracket = 0;
756 0 : nPosClosingBracket = 1;
757 : }
758 :
759 : // Without commas
760 0 : else if (nElements == 1) {
761 0 : nPosOpeningBracket = 0;
762 0 : nCharsCopied = MIN(1024, strlen(papszElements[0]));
763 : memcpy(pszPQarray + sizeof(char),
764 0 : papszElements[0], nCharsCopied*sizeof(char));
765 :
766 0 : nPosClosingBracket = nCharsCopied + sizeof(char);
767 : }
768 :
769 : // loop the array and create the string
770 : else {
771 0 : nPosOpeningBracket = 0;
772 :
773 0 : nNumberOfCommas = nElements - 1;
774 0 : nNumAvailableBytes = 1024 - (2 + nNumberOfCommas) * sizeof(char);
775 :
776 : // This should NEVER happen, it's really difficult...
777 0 : if (nNumAvailableBytes < 2) {
778 : CPLError(CE_Failure, CPLE_OutOfMemory,
779 0 : "Sorry, couldn't allocate enough space for PQ array");
780 0 : CPLFree(pszPQarray);
781 :
782 0 : return NULL;
783 : }
784 :
785 : // ptrString will move over the output string
786 0 : ptrString = pszPQarray + sizeof(char);
787 :
788 0 : for(int i = 0; i < nElements; i++) {
789 :
790 : // Check if element really exists
791 0 : if (papszElements[i] == NULL) {
792 : // one less comma, more space available
793 0 : nNumAvailableBytes += sizeof(char);
794 0 : continue;
795 : }
796 :
797 : // Copy the element to output string
798 : else {
799 0 : memset(szTemp, '\0', 1024*sizeof(char));
800 :
801 : // We can copy, at most, the number of available bytes
802 0 : nCharsCopied = MIN(strlen(papszElements[i]), nNumAvailableBytes);
803 0 : memcpy(szTemp, papszElements[i], nCharsCopied * sizeof(char));
804 :
805 : // Copy element to end string
806 0 : memcpy(ptrString, szTemp, strlen(szTemp) * sizeof(char));
807 0 : ptrString += strlen(szTemp) * sizeof(char);
808 :
809 : // No last element, add comma
810 0 : if (i < nElements - 1) {
811 : // If we don't have space for more elements, break
812 0 : if (nNumAvailableBytes == 1) {
813 0 : break;
814 : }
815 :
816 : // Add comma
817 0 : ptrString[0] = ',';
818 0 : ptrString += sizeof(char);
819 : }
820 : }
821 : }
822 :
823 0 : nPosClosingBracket = ptrString - pszPQarray;
824 : }
825 :
826 : // put brackets and return
827 0 : pszPQarray[nPosOpeningBracket] = '{';
828 0 : pszPQarray[nPosClosingBracket] = '}';
829 0 : return pszPQarray;
830 : }
831 :
832 : /**
833 : * Method: Open
834 : * Open a connection wiht PostgreSQL. The connection string will have this
835 : * format:
836 : * PG:[host='<host>]' user='<user>' [password='<password>]'
837 : * dbname='<dbname>' table='<raster_table>' [mode='working_mode'
838 : * where='<where_clause>'].
839 : * All the connection string, apart from table='table_name', where='sql_where'
840 : * and mode ='working_mode' is PQconnectdb-valid.
841 : *
842 : * NOTE: The table name can't include the schema in the form
843 : * <schema>.<table_name>,
844 : * and this is a TODO task.
845 : * Parameters:
846 : * - GDALOpenInfo *: pointer to a GDALOpenInfo structure, with useful
847 : * information for Dataset
848 : * Returns:
849 : * - GDALDataset *: pointer to a new GDALDataset instance on success,
850 : * NULL otherwise
851 : */
852 8417 : GDALDataset * WKTRasterDataset::Open(GDALOpenInfo * poOpenInfo) {
853 8417 : GBool bTableHasIndex = FALSE;
854 8417 : GBool bTableHasRegularBlocking = TRUE;
855 8417 : GBool bExistsOverviewsTable = FALSE;
856 8417 : char * pszRasterColumnName = NULL;
857 8417 : char * pszSchemaName = NULL;
858 8417 : char * pszTableName = NULL;
859 8417 : char * pszWhereClause = NULL;
860 8417 : WKTRasterDataset * poDS = NULL;
861 8417 : PGconn * hPGconn = NULL;
862 8417 : PGresult * hPGresult = NULL;
863 8417 : CPLString osCommand;
864 : char * pszWorkingMode;
865 8417 : char * pszArrayPixelTypes = NULL;
866 8417 : char * pszArrayNodataValues = NULL;
867 8417 : char ** papszPixelTypes = NULL;
868 8417 : char ** papszNodataValues = NULL;
869 8417 : GDALDataType hDataType = GDT_Byte;
870 8417 : double dfNoDataValue = 0.0;
871 8417 : int nCountNoDataValues = 0;
872 :
873 :
874 : /********************************************************
875 : * Security checkings
876 : ********************************************************/
877 8417 : if (
878 : poOpenInfo->pszFilename == NULL ||
879 : !EQUALN(poOpenInfo->pszFilename, "PG:", 3)
880 : ) {
881 : /**
882 : * Drivers must quietly return NULL if the passed file is not of their
883 : * format. They should only produce an error if the file does appear to
884 : * be of their supported format, but for some reason, unsupported or
885 : * corrupt
886 : */
887 8416 : return NULL;
888 : }
889 :
890 :
891 : /*************************************************
892 : * Check GEOS library existence
893 : *************************************************/
894 1 : if (OGRGeometryFactory::haveGEOS() == FALSE) {
895 : CPLError(CE_Failure, CPLE_AppDefined,
896 0 : "Couldn't find GEOS library installed");
897 0 : return NULL;
898 : }
899 :
900 :
901 : /*************************************************
902 : * Only support for Readonly access, from now
903 : * UPDATE July 09: GA_Update allowed
904 : *************************************************/
905 : /**
906 : if (poOpenInfo->eAccess != GA_ReadOnly) {
907 : CPLError(CE_Failure, CPLE_AppDefined,
908 : "UNDER DEVELOPMENT: Access levels supported: readonly");
909 : return NULL;
910 : }
911 : */
912 :
913 :
914 : /************************************************
915 : * Extract table name, mode and where clause
916 : ***********************************************/
917 1 : char* pszConnectionString = CPLStrdup(poOpenInfo->pszFilename);
918 :
919 1 : pszWhereClause = ExtractField(&pszConnectionString, "where=");
920 :
921 : // Working mode by default: REGULARLY_TILED_MODE
922 1 : pszWorkingMode = ExtractField(&pszConnectionString, "mode=");
923 1 : if (pszWorkingMode == NULL)
924 1 : pszWorkingMode = CPLStrdup(REGULARLY_TILED_MODE);
925 :
926 :
927 1 : pszSchemaName = ExtractField(&pszConnectionString, "schema=");
928 1 : if (pszSchemaName == NULL)
929 0 : pszSchemaName = CPLStrdup(DEFAULT_SCHEMA);
930 :
931 1 : pszTableName = ExtractField(&pszConnectionString, "table=");
932 :
933 1 : if (pszTableName == NULL) {
934 : CPLError(CE_Failure, CPLE_AppDefined,
935 : "Can't find a table name. Is connection string in the format \
936 : PG:[host=<host>] [user=<user>] [password=<password>] \
937 : dbname=<dbname> table=<raster_table> [schema=<schema>] \
938 0 : [mode=<working_mode>] [where=<where_clause>] ?\n");
939 0 : CPLFree(pszSchemaName);
940 0 : CPLFree(pszWorkingMode);
941 0 : CPLFree(pszWhereClause);
942 0 : CPLFree(pszConnectionString);
943 :
944 0 : return NULL;
945 :
946 : }
947 :
948 : /*********************************************************
949 : * Check working mode
950 : *********************************************************/
951 1 : if (!EQUALN(pszWorkingMode, REGULARLY_TILED_MODE,
952 : strlen(REGULARLY_TILED_MODE))) {
953 : CPLError(CE_Failure, CPLE_AppDefined,
954 : "Sorry, but the only working mode accepted from now is the \
955 0 : REGULARLY_TILED_MODE\n");
956 0 : CPLFree(pszTableName);
957 0 : CPLFree(pszSchemaName);
958 0 : CPLFree(pszWorkingMode);
959 0 : CPLFree(pszWhereClause);
960 0 : CPLFree(pszConnectionString);
961 :
962 0 : return NULL;
963 : }
964 :
965 : /**
966 : * After this, the connection string is prepared to be directly passed
967 : * to PGconnect
968 : */
969 :
970 :
971 : /********************************************************
972 : * Open a database connection
973 : ********************************************************/
974 1 : hPGconn = OpenConnection(pszConnectionString);
975 :
976 1 : CPLFree(pszConnectionString);
977 1 : pszConnectionString = NULL;
978 :
979 1 : if (hPGconn == NULL) {
980 1 : CPLFree(pszTableName);
981 1 : CPLFree(pszSchemaName);
982 1 : CPLFree(pszWorkingMode);
983 1 : CPLFree(pszWhereClause);
984 :
985 1 : return NULL;
986 : }
987 :
988 :
989 : /************************************************
990 : * Check if the RASTER_COLUMNS table does exist
991 : ************************************************/
992 0 : if (ExistsRasterColumnsTable(hPGconn) == FALSE) {
993 : CPLError(CE_Failure, CPLE_AppDefined,
994 : "Couldn't find RASTER_COLUMNS table. Please, check WKT Raster\
995 0 : extension is properly installed");
996 :
997 0 : PQfinish(hPGconn);
998 0 : CPLFree(pszTableName);
999 0 : CPLFree(pszSchemaName);
1000 0 : CPLFree(pszWorkingMode);
1001 0 : CPLFree(pszWhereClause);
1002 :
1003 0 : return NULL;
1004 : }
1005 :
1006 :
1007 :
1008 : /******************************************************
1009 : * Check if the table has a column of the type raster
1010 : ******************************************************/
1011 : pszRasterColumnName = GetWKTRasterColumnName(hPGconn, pszSchemaName,
1012 0 : pszTableName);
1013 0 : if (pszRasterColumnName == NULL) {
1014 : CPLError(CE_Failure, CPLE_AppDefined,
1015 : "Can't find a WKT Raster column in %s table\n",
1016 0 : pszTableName);
1017 0 : PQfinish(hPGconn);
1018 0 : CPLFree(pszTableName);
1019 0 : CPLFree(pszSchemaName);
1020 0 : CPLFree(pszWorkingMode);
1021 0 : CPLFree(pszWhereClause);
1022 :
1023 0 : return NULL;
1024 : }
1025 :
1026 :
1027 : /****************************************************
1028 : * Check if the table has regular_blocking
1029 : * constraint enabled. From now (July 2009)
1030 : * only regular_blocking tables are allowed.
1031 : *
1032 : * TODO: The fact that "regular_blocking"
1033 : * is set to TRUE doesn't mean that all the
1034 : * tiles read are regular. Need to check it
1035 : ****************************************************/
1036 : bTableHasRegularBlocking = TableHasRegularBlocking(hPGconn, pszTableName,
1037 0 : pszRasterColumnName, pszSchemaName);
1038 0 : if (bTableHasRegularBlocking == FALSE) {
1039 : CPLError(CE_Failure, CPLE_AppDefined,
1040 : "Sorry, but table %s doesn't seem to have regular blocking \
1041 : arrangement. Only tables with regular blocking arrangement can \
1042 : be read from now\n",
1043 0 : pszTableName);
1044 0 : PQfinish(hPGconn);
1045 0 : CPLFree(pszTableName);
1046 0 : CPLFree(pszSchemaName);
1047 0 : CPLFree(pszWorkingMode);
1048 0 : CPLFree(pszWhereClause);
1049 0 : CPLFree(pszRasterColumnName);
1050 0 : return NULL;
1051 :
1052 : }
1053 :
1054 : /********************************************
1055 : * Check if the table has an index
1056 : * Useful for spatial queries in raster band
1057 : ********************************************/
1058 0 : bTableHasIndex = TableHasGISTIndex(hPGconn, pszTableName, pszSchemaName);
1059 :
1060 :
1061 : /*********************************************************************
1062 : * At this point, we can start working with the given database and
1063 : * table. Of course, the working mode is essential from now.
1064 : *
1065 : * TODO: Only REGULARLY_TILED_WORKING_MODE (this is, regular_blocking) is
1066 : * allowed, but the code must be prepared for supporting the rest of the
1067 : * working modes
1068 : *********************************************************************/
1069 :
1070 :
1071 :
1072 : /********************************************************
1073 : * Instantiate a new Dataset object
1074 : ********************************************************/
1075 0 : poDS = new WKTRasterDataset();
1076 0 : poDS->pszRasterColumnName = pszRasterColumnName;
1077 0 : poDS->bTableHasGISTIndex = bTableHasIndex;
1078 0 : poDS->pszTableName = pszTableName;
1079 0 : poDS->pszSchemaName = pszSchemaName;
1080 0 : poDS->pszWorkingMode = pszWorkingMode;
1081 0 : poDS->pszWhereClause = pszWhereClause;
1082 0 : poDS->eAccess = poOpenInfo->eAccess;
1083 0 : poDS->hPGconn = hPGconn;
1084 0 : poDS->bCloseConnection = TRUE; // dataset must close connection
1085 :
1086 :
1087 : /********************************************
1088 : * Populate georeference related fields
1089 : ********************************************/
1090 0 : if (poDS->SetRasterProperties() == CE_Failure) {
1091 : CPLError(CE_Failure, CPLE_AppDefined,
1092 : "Sorry, but couldn't create the structure to read the \
1093 0 : raster on memory, aborting.");
1094 0 : delete poDS;
1095 0 : return NULL;
1096 : }
1097 :
1098 : /****************************************************************
1099 : * Get pixel_types and nodata string representations (PQ arrays)
1100 : ****************************************************************/
1101 : osCommand.Printf(
1102 : "select pixel_types, nodata_values from raster_columns \
1103 : where r_table_schema = '%s' and r_table_name = '%s' and \
1104 : r_column = '%s'",
1105 0 : poDS->pszSchemaName, poDS->pszTableName, poDS->pszRasterColumnName);
1106 0 : hPGresult = PQexec(poDS->hPGconn, osCommand.c_str());
1107 0 : if (
1108 : hPGresult == NULL ||
1109 : PQresultStatus(hPGresult) != PGRES_TUPLES_OK ||
1110 : PQntuples(hPGresult) <= 0
1111 : ) {
1112 0 : poDS->nBands = 1;
1113 : } else {
1114 0 : pszArrayPixelTypes = CPLStrdup(PQgetvalue(hPGresult, 0, 0));
1115 0 : pszArrayNodataValues = CPLStrdup(PQgetvalue(hPGresult, 0, 1));
1116 0 : int nBands = 1;
1117 :
1118 : // Explode PQ arrays into char **
1119 0 : if (pszArrayPixelTypes != NULL)
1120 : papszPixelTypes = poDS->ExplodeArrayString(pszArrayPixelTypes,
1121 0 : &nBands);
1122 0 : if (pszArrayNodataValues != NULL) {
1123 :
1124 : /**
1125 : * Has nBands been fetched by the previous call? If yes, we don't
1126 : * need to fetch it in the next call
1127 : */
1128 0 : if (nBands != 0)
1129 : papszNodataValues =
1130 : poDS->ExplodeArrayString(pszArrayNodataValues,
1131 0 : &nCountNoDataValues);
1132 :
1133 : // We don't have a value for nBands, try to fetch one
1134 : else
1135 : papszNodataValues =
1136 : poDS->ExplodeArrayString(pszArrayNodataValues,
1137 0 : &nBands);
1138 : }
1139 :
1140 0 : poDS->nBands = nBands;
1141 : }
1142 :
1143 :
1144 : /**************************************************************
1145 : * Create the raster bands
1146 : **************************************************************/
1147 0 : GBool bSignedByte = FALSE;
1148 0 : int nBitDepth = 8;
1149 0 : for (int iBand = 0; iBand < poDS->nBands; iBand++) {
1150 0 : if (pszArrayPixelTypes != NULL && papszPixelTypes != NULL) {
1151 0 : if (EQUALN(papszPixelTypes[iBand], "1BB", 3 * sizeof (char))) {
1152 0 : hDataType = GDT_Byte;
1153 0 : nBitDepth = 1;
1154 : }
1155 :
1156 0 : else if(EQUALN(papszPixelTypes[iBand], "2BUI", 4 * sizeof (char))) {
1157 0 : hDataType = GDT_Byte;
1158 0 : nBitDepth = 2;
1159 : }
1160 :
1161 0 : else if (EQUALN(papszPixelTypes[iBand], "4BUI",
1162 : 4 * sizeof (char))) {
1163 0 : hDataType = GDT_Byte;
1164 0 : nBitDepth = 4;
1165 : }
1166 :
1167 0 : else if (EQUALN(papszPixelTypes[iBand], "8BUI",
1168 : 4 * sizeof (char))) {
1169 0 : hDataType = GDT_Byte;
1170 0 : nBitDepth = 8;
1171 : }
1172 :
1173 0 : else if (EQUALN(papszPixelTypes[iBand], "8BSI",
1174 : 4 * sizeof (char))) {
1175 0 : hDataType = GDT_Byte;
1176 :
1177 : /**
1178 : * To indicate the unsigned byte values between 128 and 255
1179 : * should be interpreted as being values between -128 and -1 for
1180 : * applications that recognise the SIGNEDBYTE type.
1181 : */
1182 0 : bSignedByte = TRUE;
1183 :
1184 0 : nBitDepth = 8;
1185 : }
1186 :
1187 0 : else if (EQUALN(papszPixelTypes[iBand], "16BSI",
1188 : 5 * sizeof (char))) {
1189 0 : hDataType = GDT_Int16;
1190 0 : nBitDepth = 16;
1191 : }
1192 :
1193 0 : else if (EQUALN(papszPixelTypes[iBand], "16BUI",
1194 : 5 * sizeof (char))) {
1195 0 : hDataType = GDT_UInt16;
1196 0 : nBitDepth = 16;
1197 : }
1198 :
1199 0 : else if (EQUALN(papszPixelTypes[iBand], "32BSI",
1200 : 5 * sizeof (char))) {
1201 0 : hDataType = GDT_Int32;
1202 0 : nBitDepth = 32;
1203 : }
1204 :
1205 0 : else if (EQUALN(papszPixelTypes[iBand], "32BUI",
1206 : 5 * sizeof (char))) {
1207 0 : hDataType = GDT_UInt32;
1208 0 : nBitDepth = 32;
1209 : }
1210 :
1211 :
1212 0 : else if (EQUALN(papszPixelTypes[iBand], "32BF",
1213 : 4 * sizeof (char))) {
1214 0 : hDataType = GDT_Float32;
1215 0 : nBitDepth = 32;
1216 : }
1217 :
1218 0 : else if (EQUALN(papszPixelTypes[iBand], "64BF",
1219 : 4 * sizeof (char))) {
1220 0 : hDataType = GDT_Float64;
1221 0 : nBitDepth = 64;
1222 : }
1223 :
1224 : else {
1225 0 : hDataType = GDT_Byte;
1226 0 : nBitDepth = 8;
1227 : }
1228 :
1229 : }
1230 :
1231 : // array of data types doesn't exist or is an empty string
1232 : else {
1233 0 : hDataType = GDT_Byte;
1234 0 : nBitDepth = 8;
1235 : }
1236 :
1237 0 : if (pszArrayNodataValues != NULL && papszNodataValues != NULL &&
1238 : iBand < nCountNoDataValues) {
1239 0 : dfNoDataValue = atof(papszNodataValues[iBand]);
1240 : }
1241 :
1242 : // array of nodata values doesn't exist or is an empty string
1243 : else {
1244 0 : dfNoDataValue = 0.0;
1245 : }
1246 :
1247 :
1248 : // Create raster band objects
1249 : poDS->SetBand(iBand + 1, new WKTRasterRasterBand(poDS, iBand + 1,
1250 0 : hDataType, dfNoDataValue, bSignedByte, nBitDepth));
1251 : }
1252 :
1253 :
1254 : /******************************************************
1255 : * Create overviews datasets, if does exist
1256 : ******************************************************/
1257 0 : bExistsOverviewsTable = ExistsOverviewsTable(hPGconn);
1258 0 : if (bExistsOverviewsTable) {
1259 :
1260 :
1261 : // Count the number of overviews
1262 : osCommand.Printf(
1263 : "select o_table_name, overview_factor, o_column, \
1264 : o_table_schema from raster_overviews where \
1265 : r_table_schema = '%s' and r_table_name = '%s'",
1266 0 : poDS->pszSchemaName, poDS->pszTableName);
1267 0 : hPGresult = PQexec(poDS->hPGconn, osCommand.c_str());
1268 :
1269 : // no overviews
1270 0 : if (
1271 : hPGresult == NULL ||
1272 : PQresultStatus(hPGresult) != PGRES_TUPLES_OK ||
1273 : PQntuples(hPGresult) <= 0
1274 : ) {
1275 0 : poDS->nOverviews = 0;
1276 : }
1277 : // overviews
1278 : else {
1279 0 : poDS->nOverviews = PQntuples(hPGresult);
1280 : }
1281 :
1282 :
1283 : /**
1284 : * Create overviews's Datasets and metadata
1285 : */
1286 0 : if (poDS->nOverviews > 0) {
1287 :
1288 : poDS->papoWKTRasterOv =
1289 : (WKTRasterDataset **) VSICalloc(poDS->nOverviews,
1290 0 : sizeof (WKTRasterDataset *));
1291 0 : if (poDS->papoWKTRasterOv == NULL) {
1292 : CPLError(CE_Warning, CPLE_ObjectNull,
1293 : "Couldn't allocate memory for overviews. Number of \
1294 0 : overviews will be set to 0");
1295 0 : poDS->nOverviews = 0;
1296 : }
1297 : else {
1298 :
1299 0 : for (int i = 0; i < poDS->nOverviews; i++) {
1300 : /**
1301 : * CREATE DATASETS
1302 : */
1303 :
1304 : // AND RASTERXSIZE?? IT SHOULD BE REDUCED...
1305 0 : poDS->papoWKTRasterOv[i] = new WKTRasterDataset();
1306 0 : poDS->papoWKTRasterOv[i]->nBlockSizeX =
1307 0 : poDS->nBlockSizeX;
1308 0 : poDS->papoWKTRasterOv[i]->nBlockSizeY =
1309 0 : poDS->nBlockSizeY;
1310 :
1311 0 : int nOverviewFactor = atoi(PQgetvalue(hPGresult, i, 1));
1312 : char * pszOVTableName = CPLStrdup(PQgetvalue(hPGresult, i,
1313 0 : 0));
1314 : char * pszOVColumnName = CPLStrdup(PQgetvalue(hPGresult, i,
1315 0 : 2));
1316 : char * pszOVSchemaName = CPLStrdup(PQgetvalue(hPGresult, i,
1317 0 : 3));
1318 :
1319 : /**
1320 : * m/px is bigger in ov, because we are more far...
1321 : */
1322 0 : poDS->papoWKTRasterOv[i]->dfPixelSizeX =
1323 0 : poDS->dfPixelSizeX * nOverviewFactor;
1324 :
1325 0 : poDS->papoWKTRasterOv[i]->dfPixelSizeY =
1326 0 : poDS->dfPixelSizeY * nOverviewFactor;
1327 :
1328 : /**
1329 : * But raster size is smaller, for the same reason
1330 : * (we're more far covering the same area)
1331 : */
1332 0 : poDS->papoWKTRasterOv[i]->nRasterXSize =
1333 0 : poDS->nRasterXSize / nOverviewFactor;
1334 :
1335 0 : poDS->papoWKTRasterOv[i]->nRasterYSize =
1336 0 : poDS->nRasterYSize / nOverviewFactor;
1337 :
1338 : // Check these values (are correct??)
1339 0 : poDS->papoWKTRasterOv[i]->dfRotationX = poDS->dfRotationX;
1340 0 : poDS->papoWKTRasterOv[i]->dfRotationY = poDS->dfRotationY;
1341 0 : poDS->papoWKTRasterOv[i]->dfUpperLeftX = poDS->dfUpperLeftX;
1342 0 : poDS->papoWKTRasterOv[i]->dfUpperLeftY = poDS->dfUpperLeftY;
1343 :
1344 :
1345 0 : poDS->papoWKTRasterOv[i]->hPGconn = poDS->hPGconn;
1346 :
1347 : // children datasets don't close connection
1348 0 : poDS->papoWKTRasterOv[i]->bCloseConnection = FALSE;
1349 0 : poDS->papoWKTRasterOv[i]->pszTableName = pszOVTableName;
1350 0 : poDS->papoWKTRasterOv[i]->pszSchemaName = pszOVSchemaName;
1351 0 : poDS->papoWKTRasterOv[i]->pszRasterColumnName =
1352 0 : pszOVColumnName;
1353 :
1354 0 : poDS->papoWKTRasterOv[i]->pszWorkingMode =
1355 0 : CPLStrdup(poDS->pszWorkingMode);
1356 0 : poDS->papoWKTRasterOv[i]->pszWhereClause =
1357 : (poDS->pszWhereClause != NULL) ?
1358 : CPLStrdup(poDS->pszWhereClause):
1359 0 : NULL;
1360 :
1361 : GBool bOVTableHasIndex = TableHasGISTIndex(poDS->hPGconn,
1362 0 : pszOVTableName, pszOVSchemaName);
1363 0 : poDS->papoWKTRasterOv[i]->bTableHasGISTIndex =
1364 0 : bOVTableHasIndex;
1365 :
1366 0 : poDS->papoWKTRasterOv[i]->nSrid = poDS->nSrid;
1367 :
1368 :
1369 : /**
1370 : * CREATE OVERVIEWS RASTER BANDS
1371 : */
1372 0 : poDS->papoWKTRasterOv[i]->nBands = poDS->nBands;
1373 0 : for (int j = 0; j < poDS->nBands; j++) {
1374 : WKTRasterRasterBand * poWKTRB =
1375 : (WKTRasterRasterBand *)
1376 0 : (poDS->GetRasterBand(j + 1));
1377 0 : poDS->papoWKTRasterOv[i]->SetBand(j + 1,
1378 : new WKTRasterRasterBand(
1379 : poDS->papoWKTRasterOv[i],
1380 : j + 1, poWKTRB->GetRasterDataType(),
1381 0 : poWKTRB->GetNoDataValue(),
1382 : poWKTRB->IsSignedByteDataType(),
1383 0 : poWKTRB->GetNBitDepth()));
1384 : }
1385 :
1386 : } // end of creating overviews
1387 :
1388 : } // end else
1389 :
1390 : // free result
1391 0 : PQclear(hPGresult);
1392 :
1393 : } // end if nOverviews > 0
1394 : } // end if exists overview table
1395 :
1396 : /****************** Overviews created **********************/
1397 :
1398 :
1399 : /***********************************************************
1400 : * Free memory
1401 : ***********************************************************/
1402 0 : CSLDestroy(papszPixelTypes);
1403 0 : CSLDestroy(papszNodataValues);
1404 0 : CPLFree(pszArrayPixelTypes);
1405 0 : CPLFree(pszArrayNodataValues);
1406 :
1407 : // All WKT Raster bans are consecutives
1408 0 : poDS->SetMetadataItem("INTERLEAVE", "BAND", "IMAGE_STRUCTURE");
1409 :
1410 0 : return poDS;
1411 : }
1412 :
1413 : /**
1414 : * Method: GetGeoTransform
1415 : * Fetches the coefficients for transforming between pixel/line (P,L) raster
1416 : * space and projection coordinates (Xp, Yp) space
1417 : * The affine transformation performed is:
1418 : * Xp = padfTransform[0] + P*padfTransform[1] + L*padfTransform[2];
1419 : * Yp = padfTransform[3] + P*padfTransform[4] + L*padfTransform[5];
1420 : * Parameters:
1421 : * - double *: pointer to a double array, that will contains the affine
1422 : * transformation matrix coefficients
1423 : * Returns:
1424 : * - CPLErr: CE_None on success, or CE_Failure if no transform can be fetched.
1425 : */
1426 0 : CPLErr WKTRasterDataset::GetGeoTransform(double * padfTransform) {
1427 : // checking input parameters
1428 : // NOT NEEDED (Is illegal to call GetGeoTransform with a NULL
1429 : // argument. Thanks to Even Rouault)
1430 : /*
1431 : if (padfTransform == NULL)
1432 : {
1433 : // default matrix
1434 : padfTransform[0] = 0.0;
1435 : padfTransform[1] = 1.0;
1436 : padfTransform[2] = 0.0;
1437 : padfTransform[3] = 0.0;
1438 : padfTransform[4] = 0.0;
1439 : padfTransform[5] = 1.0;
1440 :
1441 : return CE_Failure;
1442 : }
1443 : */
1444 :
1445 : // copy necessary values in supplied buffer
1446 0 : padfTransform[0] = dfUpperLeftX;
1447 0 : padfTransform[1] = dfPixelSizeX;
1448 0 : padfTransform[2] = dfRotationX;
1449 0 : padfTransform[3] = dfUpperLeftY;
1450 0 : padfTransform[4] = dfRotationY;
1451 0 : padfTransform[5] = dfPixelSizeY;
1452 :
1453 0 : return CE_None;
1454 : }
1455 :
1456 : /**
1457 : * Fetch the projection definition string for this dataset in OpenGIS WKT
1458 : * format. It should be suitable for use with the OGRSpatialReference class
1459 : * Parameters: none
1460 : * Returns:
1461 : * - const char *: a pointer to an internal projection reference string.
1462 : * It should not be altered, freed or expected to last for
1463 : * long.
1464 : */
1465 0 : const char * WKTRasterDataset::GetProjectionRef() {
1466 0 : CPLString osCommand;
1467 : PGresult * hResult;
1468 :
1469 0 : if (nSrid == -1) {
1470 0 : return "";
1471 : }
1472 :
1473 : /********************************************************
1474 : * Reading proj from database
1475 : ********************************************************/
1476 0 : osCommand.Printf("SELECT srtext FROM spatial_ref_sys where SRID=%d", nSrid);
1477 :
1478 0 : hResult = PQexec(this->hPGconn, osCommand.c_str());
1479 :
1480 0 : if (hResult && PQresultStatus(hResult) == PGRES_TUPLES_OK
1481 : && PQntuples(hResult) > 0) {
1482 0 : CPLFree(pszProjection);
1483 0 : pszProjection = CPLStrdup(PQgetvalue(hResult, 0, 0));
1484 : }
1485 :
1486 0 : if (hResult)
1487 0 : PQclear(hResult);
1488 :
1489 0 : return pszProjection;
1490 : }
1491 :
1492 : /**
1493 : * Set the projection reference string for this dataset. The string should
1494 : * be in OGC WKT or PROJ.4 format
1495 : * Parameters:
1496 : * - const char *: projection reference string.
1497 : * Returns:
1498 : * - CE_Failure if an error occurs, otherwise CE_None.
1499 : */
1500 0 : CPLErr WKTRasterDataset::SetProjection(const char * pszProjectionRef)
1501 : {
1502 0 : VALIDATE_POINTER1(pszProjectionRef, "SetProjection", CE_Failure);
1503 :
1504 0 : CPLString osCommand;
1505 : PGresult * hResult;
1506 0 : int nFetchedSrid = -1;
1507 :
1508 :
1509 : /*****************************************************************
1510 : * Check if the dataset allows updating
1511 : *****************************************************************/
1512 0 : if (GetAccess() != GA_Update) {
1513 : CPLError(CE_Failure, CPLE_NoWriteAccess,
1514 0 : "This driver doesn't allow write access");
1515 0 : return CE_Failure;
1516 : }
1517 :
1518 : /*****************************************************************
1519 : * Look for projection with this text
1520 : *****************************************************************/
1521 :
1522 : // First, WKT text
1523 : osCommand.Printf("SELECT srid FROM spatial_ref_sys where srtext='%s'",
1524 0 : pszProjectionRef);
1525 0 : hResult = PQexec(hPGconn, osCommand.c_str());
1526 :
1527 0 : if (hResult && PQresultStatus(hResult) == PGRES_TUPLES_OK
1528 : && PQntuples(hResult) > 0) {
1529 :
1530 0 : nFetchedSrid = atoi(PQgetvalue(hResult, 0, 0));
1531 :
1532 : // update class attribute
1533 0 : nSrid = nFetchedSrid;
1534 :
1535 :
1536 : // update raster_columns table
1537 : osCommand.Printf("UPDATE raster_columns SET srid=%d WHERE \
1538 : r_table_name = '%s' AND r_column = '%s'",
1539 0 : nSrid, pszTableName, pszRasterColumnName);
1540 0 : hResult = PQexec(hPGconn, osCommand.c_str());
1541 0 : if (hResult == NULL || PQresultStatus(hResult) != PGRES_COMMAND_OK) {
1542 : CPLError(CE_Failure, CPLE_AppDefined,
1543 : "Couldn't update raster_columns table: %s",
1544 0 : PQerrorMessage(hPGconn));
1545 0 : return CE_Failure;
1546 : }
1547 :
1548 : // TODO: Update ALL blocks with the new srid...
1549 :
1550 0 : return CE_None;
1551 : }
1552 :
1553 : // If not, proj4 text
1554 : else {
1555 : osCommand.Printf(
1556 : "SELECT srid FROM spatial_ref_sys where proj4text='%s'",
1557 0 : pszProjectionRef);
1558 0 : hResult = PQexec(this->hPGconn, osCommand.c_str());
1559 :
1560 0 : if (hResult && PQresultStatus(hResult) == PGRES_TUPLES_OK
1561 : && PQntuples(hResult) > 0) {
1562 :
1563 0 : nFetchedSrid = atoi(PQgetvalue(hResult, 0, 0));
1564 :
1565 : // update class attribute
1566 0 : nSrid = nFetchedSrid;
1567 :
1568 : // update raster_columns table
1569 : osCommand.Printf("UPDATE raster_columns SET srid=%d WHERE \
1570 : r_table_name = '%s' AND r_column = '%s'",
1571 0 : nSrid, pszTableName, pszRasterColumnName);
1572 :
1573 0 : hResult = PQexec(hPGconn, osCommand.c_str());
1574 0 : if (hResult == NULL ||
1575 : PQresultStatus(hResult) != PGRES_COMMAND_OK) {
1576 : CPLError(CE_Failure, CPLE_AppDefined,
1577 : "Couldn't update raster_columns table: %s",
1578 0 : PQerrorMessage(hPGconn));
1579 0 : return CE_Failure;
1580 : }
1581 :
1582 : // TODO: Update ALL blocks with the new srid...
1583 :
1584 0 : return CE_None;
1585 : }
1586 :
1587 : else {
1588 : CPLError(CE_Failure, CPLE_WrongFormat,
1589 0 : "Couldn't find WKT neither proj4 definition");
1590 0 : return CE_Failure;
1591 : }
1592 0 : }
1593 : }
1594 :
1595 : /**
1596 : * Set the affine transformation coefficients.
1597 : * Parameters:
1598 : * - double *:a six double buffer containing the transformation
1599 : * coefficients to be written with the dataset.
1600 : * Returns:
1601 : * - CE_None on success, or CE_Failure if this transform cannot be written.
1602 : */
1603 0 : CPLErr WKTRasterDataset::SetGeoTransform(double * padfTransform)
1604 : {
1605 0 : VALIDATE_POINTER1(padfTransform, "SetGeoTransform", CE_Failure);
1606 :
1607 0 : dfUpperLeftX = padfTransform[0];
1608 0 : dfPixelSizeX = padfTransform[1];
1609 0 : dfRotationX = padfTransform[2];
1610 0 : dfUpperLeftY = padfTransform[3];
1611 0 : dfRotationY = padfTransform[4];
1612 0 : dfPixelSizeY = padfTransform[5];
1613 :
1614 : // TODO: Update all data blocks with these values
1615 :
1616 0 : return CE_None;
1617 : }
1618 :
1619 :
1620 : /************************************************************************/
1621 : /* GDALRegister_WKTRaster() */
1622 : /************************************************************************/
1623 338 : void GDALRegister_WKTRaster() {
1624 : GDALDriver *poDriver;
1625 :
1626 338 : if (GDALGetDriverByName("WKTRaster") == NULL) {
1627 336 : poDriver = new GDALDriver();
1628 :
1629 336 : poDriver->SetDescription("WKTRaster");
1630 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
1631 336 : "PostGIS WKT Raster driver");
1632 :
1633 336 : poDriver->pfnOpen = WKTRasterDataset::Open;
1634 :
1635 336 : GetGDALDriverManager()->RegisterDriver(poDriver);
1636 : }
1637 338 : }
1638 :
|
