From 9e569f9f3c690b313e3dfa41f23d450695ea7411 Mon Sep 17 00:00:00 2001 From: Sandro Santilli Date: Thu, 17 Jul 2014 08:59:53 +0000 Subject: [PATCH] Add support for raster reprojection on import (#2843) git-svn-id: http://svn.osgeo.org/postgis/trunk@12789 b70326c6-7e19-0410-871a-916f4a2858ee --- NEWS | 1 + raster/loader/raster2pgsql.c | 84 +++++++++++++----- raster/loader/raster2pgsql.h | 5 +- raster/test/regress/Makefile.in | 1 + raster/test/regress/loader/Projected-post.sql | 2 + raster/test/regress/loader/Projected-pre.sql | 8 ++ raster/test/regress/loader/Projected.opts | 1 + .../regress/loader/Projected.select.expected | 2 + .../test/regress/loader/Projected.select.sql | 2 + raster/test/regress/loader/Projected.tif | Bin 0 -> 24836 bytes 10 files changed, 82 insertions(+), 24 deletions(-) create mode 100644 raster/test/regress/loader/Projected-post.sql create mode 100644 raster/test/regress/loader/Projected-pre.sql create mode 100644 raster/test/regress/loader/Projected.opts create mode 100644 raster/test/regress/loader/Projected.select.expected create mode 100644 raster/test/regress/loader/Projected.select.sql create mode 100644 raster/test/regress/loader/Projected.tif diff --git a/NEWS b/NEWS index bc4680e75..d9952d7f1 100644 --- 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 d9fe9a167..ea30c0b2f 100644 --- a/raster/loader/raster2pgsql.c +++ b/raster/loader/raster2pgsql.c @@ -330,9 +330,10 @@ usage() { "OPTIONS:\n" )); printf(_( - " -s 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 [:] 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 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 d0579f755..15d4a0c77 100644 --- 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 457f1f7bc..cf21cb663 100644 --- 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 index 000000000..a2a66b48f --- /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 index 000000000..61918bd1d --- /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 index 000000000..98cdce4d5 --- /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 index 000000000..203172c75 --- /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 index 000000000..d910efc3a --- /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 index 0000000000000000000000000000000000000000..312692f9081dcc336dd52ba9919654a009078096 GIT binary patch literal 24836 zcmeI4O-~bX5XNU)#HvV;1T^?oC75=Jns^{;QcuPc zpMZE3?jGRF_%(2Lt?cM*hQnkJ(>x39KFoh+cAnpEN!z4dShztABAOsd>Xfj0&bLJ9 z8tXd8rg*X5k7s$qeO+UHHmT7x@9*%!DC=?7E3D5YLwW5@))W2wtk+~;UtxW|pVQ== zspLdH!)ewR`+A)9C2`zoaok1LS-%<18T~Yp>vya_ejcG=j(=nQqgSUU>wVU{y%Eah z&F=PkwcPnw9w2(aVNN;PH+Z!{J6GxGc7Zd0UL$%jmt|%@t<{b7c58XFo4#M)eA#W? zZ+Bm%8*f^jrEc2pw9#*r&H06gk7#@CVDk0B_T9hxY^YK>1`V5lT=j5{SS!729nq0zW=%jwJe?!fO7B`nbmW;?6Oa*4=ZLk^yVemMd1lrG zWW>{X@M>*GSOozP009sH0T2KI5C8!W0r@;!YrC_~myw9Ia(=BNI`Yh{3CM`2bHrNd zUF(RBJTq$oGUDkRu~vH5I-(=b%$k6VcsfU{mEN_E=*Tm(CLkl8&Jk;+cda8js%~bx zGgd(W1V8`;KmY_l00cn5LqIY?TCLkjb9aT58-5IMO00JNY z0w4eaAOHd&;2~gsyXZiEj^6x6ReR9{Enm@SZI&MMdG;!6ZNt`=!Bshx=2t|jJT7Py zbggsG|21q>MVk^;oXav*QXoQQ6-7)WmdGNnOqCRf5NHS(=DEkU@!w;^q$R0pU3Q%o fF`X|XMP8YbnOc`s6fvDIBY32Mh9n@heQo>(ch8Zw literal 0 HcmV?d00001 -- 2.40.0