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