Add support for raster reprojection on import (#2843)

git-svn-id: http://svn.osgeo.org/postgis/trunk@12789 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/NEWS b/NEWS
index bc4680e752a3e5883b2401feab2d5ce8bf311b40..d9952d7f1645365105dcc67c0dbdfb61274fb7e3 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -17,6 +17,7 @@ PostGIS 2.2.0
 
  * New Features *
 
+  - #2843, Support reprojection on raster import
   - #2349, Support for encoded_polyline input/output (Kashif Rasul)
   - #2159, report libjson version from postgis_full_version()
   - #2770, ST_MemSize(raster)

diff --git a/raster/loader/raster2pgsql.c b/raster/loader/raster2pgsql.c
index d9fe9a167b04708dbe318939181263beb7f6fb95..ea30c0b2ff1122cc8b262a3b1e231e42507ccbd0 100644 (file)
--- a/raster/loader/raster2pgsql.c
+++ b/raster/loader/raster2pgsql.c
@@ -330,9 +330,10 @@ usage() {
                "OPTIONS:\n"
        ));
        printf(_(
-               "  -s <srid> Set the raster's SRID. Defaults to %d. If SRID not\n"
-               "      provided or is %d, raster's metadata will be checked to\n"
-               "      determine an appropriate SRID.\n"
+               "  -s [<from>:]<srid> Set the SRID field. Defaults to %d.\n"
+               "     Optionally reprojects from given SRID (cannot be used with -Y).\n"
+               "     Raster's metadata will be checked to determine an appropriate SRID.\n"
+               "     If a srid of %d is provided (either as from or as target).\n"
        ), SRID_UNKNOWN, SRID_UNKNOWN);
        printf(_(
                "  -b <band> Index (1-based) of band to extract from raster. For more\n"
@@ -693,7 +694,7 @@ init_config(RTLOADERCFG *config) {
        config->overview = NULL;
        config->overview_table = NULL;
        config->quoteident = 0;
-       config->srid = SRID_UNKNOWN;
+       config->srid = config->out_srid = SRID_UNKNOWN;
        config->nband = NULL;
        config->nband_count = 0;
        memset(config->tile_size, 0, sizeof(int) * 2);
@@ -822,7 +823,7 @@ static int
 insert_records(
        const char *schema, const char *table, const char *column,
        const char *filename, const char *file_column_name,
-       int copy_statements,
+       int copy_statements, int out_srid,
        STRINGBUFFER *tileset, STRINGBUFFER *buffer
 ) {
        char *fn = NULL;
@@ -865,7 +866,7 @@ insert_records(
        }
        /* INSERT statements */
        else {
-               len = strlen("INSERT INTO  () VALUES (''::raster);") + 1;
+               len = strlen("INSERT INTO  () VALUES (ST_Transform(''::raster,xxxxxxxxx));") + 1;
                if (schema != NULL)
                        len += strlen(schema);
                len += strlen(table);
@@ -878,6 +879,7 @@ insert_records(
                        fn = strreplace(filename, "'", "''", NULL);
 
                for (x = 0; x < tileset->length; x++) {
+                       char *ptr;
                        int sqllen = len;
 
                        sqllen += strlen(tileset->line[x]);
@@ -889,17 +891,27 @@ insert_records(
                                rterror(_("insert_records: Could not allocate memory for INSERT statement"));
                                return 0;
                        }
-                       sprintf(sql, "INSERT INTO %s%s (%s%s%s) VALUES ('%s'::raster%s%s%s);",
-                               (schema != NULL ? schema : ""),
-                               table,
-                               column,
-                               (filename != NULL ? "," : ""),
-                               (filename != NULL ? file_column_name : ""),
-                               tileset->line[x],
-                               (filename != NULL ? ",'" : ""),
-                               (filename != NULL ? fn : ""),
-                               (filename != NULL ? "'" : "")
-                       );
+                       ptr = sql;
+                       ptr += sprintf(sql, "INSERT INTO %s%s (%s%s%s) VALUES (",
+                                       (schema != NULL ? schema : ""),
+                                       table,
+                                       column,
+                                       (filename != NULL ? "," : ""),
+                                       (filename != NULL ? file_column_name : "")
+                               );
+                       if (out_srid != SRID_UNKNOWN) {
+                               ptr += sprintf(ptr, "ST_Transform(");
+                       }
+                       ptr += sprintf(ptr, "'%s'::raster",
+                                       tileset->line[x]
+                               );
+                       if (out_srid != SRID_UNKNOWN) {
+                               ptr += sprintf(ptr, ", %d)", out_srid);
+                       }
+                       if (filename != NULL) {
+                               ptr += sprintf(ptr, ",'%s'", fn);
+                       }
+                       ptr += sprintf(ptr, ");");
 
                        append_sql_to_buffer(buffer, sql);
                        sql = NULL;
@@ -1507,7 +1519,7 @@ build_overview(int idx, RTLOADERCFG *config, RASTERINFO *info, int ovx, STRINGBU
                                if (!insert_records(
                                        config->schema, ovtable, config->raster_column,
                                        (config->file_column ? config->rt_filename[idx] : NULL), config->file_column_name,
-                                       config->copy_statements,
+                                       config->copy_statements, config->out_srid,
                                        tileset, buffer
                                )) {
                                        rterror(_("build_overview: Could not convert raster tiles into INSERT or COPY statements"));
@@ -1597,6 +1609,12 @@ convert_raster(int idx, RTLOADERCFG *config, RASTERINFO *info, STRINGBUFFER *til
                }
        }
 
+       if ( info->srid == SRID_UNKNOWN && config->out_srid != SRID_UNKNOWN ) {
+                 rterror(_("convert_raster: could not determine source srid, cannot transform to target srid %d"), config->out_srid);
+                 GDALClose(hdsSrc);
+                 return 0;
+       }
+
        /* record geotransform matrix */
        if (GDALGetGeoTransform(hdsSrc, info->gt) != CE_None) {
                rtinfo(_("Using default geotransform matrix (0, 1, 0, 0, 0, -1) for raster: %s"), config->rt_file[idx]);
@@ -1821,7 +1839,7 @@ convert_raster(int idx, RTLOADERCFG *config, RASTERINFO *info, STRINGBUFFER *til
                                        if (!insert_records(
                                                config->schema, config->table, config->raster_column,
                                                (config->file_column ? config->rt_filename[idx] : NULL), config->file_column_name,
-                                               config->copy_statements,
+                                               config->copy_statements, config->out_srid,
                                                tileset, buffer
                                        )) {
                                                rterror(_("convert_raster: Could not convert raster tiles into INSERT or COPY statements"));
@@ -1940,7 +1958,7 @@ convert_raster(int idx, RTLOADERCFG *config, RASTERINFO *info, STRINGBUFFER *til
                                        if (!insert_records(
                                                config->schema, config->table, config->raster_column,
                                                (config->file_column ? config->rt_filename[idx] : NULL), config->file_column_name,
-                                               config->copy_statements,
+                                               config->copy_statements, config->out_srid,
                                                tileset, buffer
                                        )) {
                                                rterror(_("convert_raster: Could not convert raster tiles into INSERT or COPY statements"));
@@ -2056,7 +2074,7 @@ process_rasters(RTLOADERCFG *config, STRINGBUFFER *buffer) {
                        if (tileset.length && !insert_records(
                                config->schema, config->table, config->raster_column,
                                (config->file_column ? config->rt_filename[i] : NULL), config->file_column_name,
-                               config->copy_statements,
+                               config->copy_statements, config->out_srid,
                                &tileset, buffer
                        )) {
                                rterror(_("process_rasters: Could not convert raster tiles into INSERT or COPY statements"));
@@ -2103,7 +2121,7 @@ process_rasters(RTLOADERCFG *config, STRINGBUFFER *buffer) {
                                        if (tileset.length && !insert_records(
                                                config->schema, config->overview_table[j], config->raster_column,
                                                (config->file_column ? config->rt_filename[i] : NULL), config->file_column_name,
-                                               config->copy_statements,
+                                               config->copy_statements, config->out_srid,
                                                &tileset, buffer
                                        )) {
                                                rterror(_("process_rasters: Could not convert overview tiles into INSERT or COPY statements"));
@@ -2300,9 +2318,18 @@ main(int argc, char **argv) {
        ****************************************************************************/
 
        for (i = 1; i < argc; i++) {
+               char *optarg, *ptr;
                /* srid */
                if (CSEQUAL(argv[i], "-s") && i < argc - 1) {
-                       config->srid = atoi(argv[++i]);
+                       optarg = argv[++i];
+                       ptr = strchr(optarg, ':');
+                       if (ptr) {
+                               *ptr++ = '\0';
+                               sscanf(optarg, "%d", &config->srid);
+                               sscanf(ptr, "%d", &config->out_srid);
+                       } else {
+                               config->srid = config->out_srid = atoi(optarg);
+                       }
                }
                /* band index */
                else if (CSEQUAL(argv[i], "-b") && i < argc - 1) {
@@ -2632,6 +2659,17 @@ main(int argc, char **argv) {
                }
        }
 
+       if (config->srid != config->out_srid) {
+               if (config->copy_statements) {
+                       rterror(_("Invalid argument combination - cannot use -Y with -s FROM_SRID:TO_SRID"));
+                       exit(1);
+               }
+               if (config->out_srid == SRID_UNKNOWN) {
+                       rterror(_("Unknown target SRID is invalid when source SRID is given"));
+                       exit(1);
+               }
+       }
+
        /* register GDAL drivers */
        GDALAllRegister();
 

diff --git a/raster/loader/raster2pgsql.h b/raster/loader/raster2pgsql.h
index d0579f755d0411c7a11d127a0f65418dc2460a80..15d4a0c773caea4f38db6a55e9cad12f3b961d51 100644 (file)
--- a/raster/loader/raster2pgsql.h
+++ b/raster/loader/raster2pgsql.h
@@ -97,9 +97,12 @@ typedef struct raster_loader_config {
        /* case-sensitive of identifiers, 1 = yes, 0 = no (default) */
        int quoteident;
 
-       /* SRID of raster */
+       /* SRID of input raster */
        int srid;
 
+       /* SRID of output raster (reprojection) */
+       int out_srid;
+
        /* bands to extract */
        int *nband; /* 1-based */
        int nband_count;

diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in
index 457f1f7bc2cf1185726fd50df4e31d924069d3dc..cf21cb663e0c47ddfa55e84a3e6f88ca3efe18ac 100644 (file)
--- a/raster/test/regress/Makefile.in
+++ b/raster/test/regress/Makefile.in
@@ -143,6 +143,7 @@ TEST_BUGS = \
 
 TEST_LOADER = \
        loader/Basic \
+       loader/Projected \
        loader/BasicCopy \
        loader/BasicFilename \
        loader/BasicOutDB \

diff --git a/raster/test/regress/loader/Projected-post.sql b/raster/test/regress/loader/Projected-post.sql
new file mode 100644 (file)
index 0000000..a2a66b4
--- /dev/null
+++ b/raster/test/regress/loader/Projected-post.sql
@@ -0,0 +1,2 @@
+TRUNCATE spatial_ref_sys;
+

diff --git a/raster/test/regress/loader/Projected-pre.sql b/raster/test/regress/loader/Projected-pre.sql
new file mode 100644 (file)
index 0000000..61918bd
--- /dev/null
+++ b/raster/test/regress/loader/Projected-pre.sql
@@ -0,0 +1,8 @@
+-- NOTE: the need for truncation here reveals a leak bug in previous tests
+TRUNCATE spatial_ref_sys;
+INSERT INTO "spatial_ref_sys" ("srid","auth_name","auth_srid","srtext","proj4text")
+VALUES
+(4326,'EPSG',4326,'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]','+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ')
+,
+('3857','EPSG','3857','PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],AUTHORITY["EPSG","3857"],AXIS["X",EAST],AXIS["Y",NORTH]]','+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs')
+;

diff --git a/raster/test/regress/loader/Projected.opts b/raster/test/regress/loader/Projected.opts
new file mode 100644 (file)
index 0000000..98cdce4
--- /dev/null
+++ b/raster/test/regress/loader/Projected.opts
@@ -0,0 +1 @@
+-s :3857

diff --git a/raster/test/regress/loader/Projected.select.expected b/raster/test/regress/loader/Projected.select.expected
new file mode 100644 (file)
index 0000000..203172c
--- /dev/null
+++ b/raster/test/regress/loader/Projected.select.expected
@@ -0,0 +1,2 @@
+3857
+BOX(-20037508 25860122,-19887373 44927335)

diff --git a/raster/test/regress/loader/Projected.select.sql b/raster/test/regress/loader/Projected.select.sql
new file mode 100644 (file)
index 0000000..d910efc
--- /dev/null
+++ b/raster/test/regress/loader/Projected.select.sql
@@ -0,0 +1,2 @@
+SELECT st_srid(rast) from loadedrast limit 1;
+SELECT st_extent(ST_SnapToGrid(rast::geometry,1)) from loadedrast;

diff --git a/raster/test/regress/loader/Projected.tif b/raster/test/regress/loader/Projected.tif
new file mode 100644 (file)
index 0000000..312692f
Binary files /dev/null and b/raster/test/regress/loader/Projected.tif differ