From: Bborie Park <bkpark at ucdavis.edu>
Date: Sat, 1 Dec 2012 01:19:45 +0000 (+0000)
Subject: Removed requirements of SRID for calling GDAL Warp API
X-Git-Tag: 2.1.0beta2~337
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=58e564f82b186cdfb6ecf8657972353d5969c6b3;p=postgis

Removed requirements of SRID for calling GDAL Warp API

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

diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c
index 763089424..07f1be1cb 100644
--- a/raster/rt_core/rt_api.c
+++ b/raster/rt_core/rt_api.c
@@ -9240,8 +9240,8 @@ rt_raster_from_gdal_dataset(GDALDatasetH ds) {
  * @return the warped raster or NULL
  */
 rt_raster rt_raster_gdal_warp(
-	rt_raster raster, const char *src_srs,
-	const char *dst_srs,
+	rt_raster raster,
+	const char *src_srs, const char *dst_srs,
 	double *scale_x, double *scale_y,
 	int *width, int *height,
 	double *ul_xw, double *ul_yw,
@@ -9288,7 +9288,6 @@ rt_raster rt_raster_gdal_warp(
 	RASTER_DEBUG(3, "starting");
 
 	assert(NULL != raster);
-	assert(NULL != src_srs);
 
 	/*
 		max_err must be gte zero
@@ -9298,20 +9297,34 @@ rt_raster rt_raster_gdal_warp(
 	if (max_err < 0.) max_err = 0.125;
 	RASTER_DEBUGF(4, "max_err = %f", max_err);
 
-	_src_srs = rt_util_gdal_convert_sr(src_srs, 0);
-	/* dst_srs not provided, set to src_srs */
-	if (NULL == dst_srs)
-		_dst_srs = rt_util_gdal_convert_sr(src_srs, 0);
-	else
-		_dst_srs = rt_util_gdal_convert_sr(dst_srs, 0);
+	/* handle srs */
+	if (src_srs != NULL) {
+		/* reprojection taking place */
+		if (dst_srs != NULL && strcmp(src_srs, dst_srs) != 0) {
+			RASTER_DEBUG(4, "Warp operation does include a reprojection");
+			_src_srs = rt_util_gdal_convert_sr(src_srs, 0);
+			_dst_srs = rt_util_gdal_convert_sr(dst_srs, 0);
+		}
+		/* no reprojection, a stub just for clarity */
+		else {
+			RASTER_DEBUG(4, "Warp operation does NOT include reprojection");
+		}
+	}
+	else if (dst_srs != NULL) {
+		/* dst_srs provided but not src_srs */
+		rterror("rt_raster_gdal_warp: SRS required for input raster if SRS provided for warped raster");
+		return NULL;
+	}
 
 	/* load raster into a GDAL MEM dataset */
 	src_ds = rt_raster_to_gdal_mem(raster, _src_srs, NULL, NULL, 0, &src_drv);
 	if (NULL == src_ds) {
 		rterror("rt_raster_gdal_warp: Unable to convert raster to GDAL MEM format");
 
-		CPLFree(_src_srs);
-		CPLFree(_dst_srs);
+		if (_src_srs != NULL) {
+			CPLFree(_src_srs);
+			CPLFree(_dst_srs);
+		}
 
 		return NULL;
 	}
@@ -9324,18 +9337,26 @@ rt_raster rt_raster_gdal_warp(
 
 		GDALClose(src_ds);
 
-		CPLFree(_src_srs);
-		CPLFree(_dst_srs);
+		if (_src_srs != NULL) {
+			CPLFree(_src_srs);
+			CPLFree(_dst_srs);
+		}
 
 		return NULL;
 	}
 	for (i = 0; i < transform_opts_len; i++) {
 		switch (i) {
 			case 1:
-				transform_opts[i] = (char *) rtalloc(sizeof(char) * (strlen("DST_SRS=") + strlen(_dst_srs) + 1));
+				if (_dst_srs != NULL)
+					transform_opts[i] = (char *) rtalloc(sizeof(char) * (strlen("DST_SRS=") + strlen(_dst_srs) + 1));
+				else
+					transform_opts[i] = (char *) rtalloc(sizeof(char) * (strlen("DST_SRS=") + 1));
 				break;
 			case 0:
-				transform_opts[i] = (char *) rtalloc(sizeof(char) * (strlen("SRC_SRS=") + strlen(_src_srs) + 1));
+				if (_src_srs != NULL)
+					transform_opts[i] = (char *) rtalloc(sizeof(char) * (strlen("SRC_SRS=") + strlen(_src_srs) + 1));
+				else
+					transform_opts[i] = (char *) rtalloc(sizeof(char) * (strlen("SRC_SRS=") + 1));
 				break;
 		}
 		if (NULL == transform_opts[i]) {
@@ -9345,35 +9366,44 @@ rt_raster rt_raster_gdal_warp(
 			rtdealloc(transform_opts);
 			GDALClose(src_ds);
 
-			CPLFree(_src_srs);
-			CPLFree(_dst_srs);
+			if (_src_srs != NULL) {
+				CPLFree(_src_srs);
+				CPLFree(_dst_srs);
+			}
 
 			return NULL;
 		}
 
 		switch (i) {
 			case 1:
-				snprintf(
-					transform_opts[i],
-					sizeof(char) * (strlen("DST_SRS=") + strlen(_dst_srs) + 1),
-					"DST_SRS=%s",
-					_dst_srs
-				);
+				if (_dst_srs != NULL) {
+					snprintf(
+						transform_opts[i],
+						sizeof(char) * (strlen("DST_SRS=") + strlen(_dst_srs) + 1),
+						"DST_SRS=%s",
+						_dst_srs
+					);
+				}
+				else
+					sprintf(transform_opts[i], "%s", "DST_SRS=");
 				break;
 			case 0:
-				snprintf(
-					transform_opts[i],
-					sizeof(char) * (strlen("SRC_SRS=") + strlen(_src_srs) + 1),
-					"SRC_SRS=%s",
-					_src_srs
-				);
+				if (_src_srs != NULL) {
+					snprintf(
+						transform_opts[i],
+						sizeof(char) * (strlen("SRC_SRS=") + strlen(_src_srs) + 1),
+						"SRC_SRS=%s",
+						_src_srs
+					);
+				}
+				else
+					sprintf(transform_opts[i], "%s", "SRC_SRS=");
 				break;
 		}
 		RASTER_DEBUGF(4, "transform_opts[%d] = %s", i, transform_opts[i]);
 	}
 	transform_opts[transform_opts_len] = NULL;
-	CPLFree(_src_srs);
-	CPLFree(_dst_srs);
+	if (_src_srs != NULL) CPLFree(_src_srs);
 
 	/* transformation object for building dst dataset */
 	transform_arg = GDALCreateGenImgProjTransformer2(src_ds, NULL, transform_opts);
@@ -9384,6 +9414,7 @@ rt_raster rt_raster_gdal_warp(
 
 		for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
 		rtdealloc(transform_opts);
+		if (_dst_srs != NULL) CPLFree(_dst_srs);
 
 		return NULL;
 	}
@@ -9399,6 +9430,7 @@ rt_raster rt_raster_gdal_warp(
 
 		for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
 		rtdealloc(transform_opts);
+		if (_dst_srs != NULL) CPLFree(_dst_srs);
 
 		return NULL;
 	}
@@ -9436,6 +9468,7 @@ rt_raster rt_raster_gdal_warp(
 
 		for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
 		rtdealloc(transform_opts);
+		if (_dst_srs != NULL) CPLFree(_dst_srs);
 
 		return NULL;
 	}
@@ -9469,6 +9502,7 @@ rt_raster rt_raster_gdal_warp(
 
 		for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
 		rtdealloc(transform_opts);
+		if (_dst_srs != NULL) CPLFree(_dst_srs);
 
 		return NULL;
 	}
@@ -9535,6 +9569,7 @@ rt_raster rt_raster_gdal_warp(
 
 			for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
 			rtdealloc(transform_opts);
+			if (_dst_srs != NULL) CPLFree(_dst_srs);
 
 			return NULL;
 		}
@@ -9565,6 +9600,7 @@ rt_raster rt_raster_gdal_warp(
 
 		for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
 		rtdealloc(transform_opts);
+		if (_dst_srs != NULL) CPLFree(_dst_srs);
 
 		return NULL;
 	}
@@ -9605,6 +9641,7 @@ rt_raster rt_raster_gdal_warp(
 
 		for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
 		rtdealloc(transform_opts);
+		if (_dst_srs != NULL) CPLFree(_dst_srs);
 
 		return NULL;
 	}
@@ -9627,6 +9664,7 @@ rt_raster rt_raster_gdal_warp(
 
 			for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
 			rtdealloc(transform_opts);
+			if (_dst_srs != NULL) CPLFree(_dst_srs);
 
 			return NULL;
 		}
@@ -9661,6 +9699,7 @@ rt_raster rt_raster_gdal_warp(
 
 				for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
 				rtdealloc(transform_opts);
+				if (_dst_srs != NULL) CPLFree(_dst_srs);
 
 				return NULL;
 			}
@@ -9678,6 +9717,7 @@ rt_raster rt_raster_gdal_warp(
 
 				for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
 				rtdealloc(transform_opts);
+				if (_dst_srs != NULL) CPLFree(_dst_srs);
 
 				return NULL;
 			}
@@ -9705,6 +9745,7 @@ rt_raster rt_raster_gdal_warp(
 
 						for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
 						rtdealloc(transform_opts);
+						if (_dst_srs != NULL) CPLFree(_dst_srs);
 
 						return NULL;
 					}
@@ -9734,6 +9775,7 @@ rt_raster rt_raster_gdal_warp(
 
 						for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
 						rtdealloc(transform_opts);
+						if (_dst_srs != NULL) CPLFree(_dst_srs);
 
 						return NULL;
 					}
@@ -9784,6 +9826,7 @@ rt_raster rt_raster_gdal_warp(
 
 				for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
 				rtdealloc(transform_opts);
+				if (_dst_srs != NULL) CPLFree(_dst_srs);
 
 				return NULL;
 			}
@@ -9814,6 +9857,7 @@ rt_raster rt_raster_gdal_warp(
 
 				for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
 				rtdealloc(transform_opts);
+				if (_dst_srs != NULL) CPLFree(_dst_srs);
 
 				return NULL;
 			}
@@ -9842,6 +9886,7 @@ rt_raster rt_raster_gdal_warp(
 
 		for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
 		rtdealloc(transform_opts);
+		if (_dst_srs != NULL) CPLFree(_dst_srs);
 
 		return NULL;
 	}
@@ -9865,6 +9910,38 @@ rt_raster rt_raster_gdal_warp(
 
 		GDALClose(src_ds);
 
+		for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
+		rtdealloc(transform_opts);
+		if (_dst_srs != NULL) CPLFree(_dst_srs);
+
+		return NULL;
+	}
+
+	/* set dst srs */
+	if (_dst_srs != NULL) {
+		cplerr = GDALSetProjection(dst_ds, _dst_srs);
+		CPLFree(_dst_srs);
+		if (cplerr != CE_None) {
+			rterror("rt_raster_gdal_warp: Unable to set projection");
+
+			GDALClose(dst_ds);
+			GDALClose(src_ds);
+
+			for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
+			rtdealloc(transform_opts);
+
+			return NULL;
+		}
+	}
+
+	/* set dst geotransform */
+	cplerr = GDALSetGeoTransform(dst_ds, _gt);
+	if (cplerr != CE_None) {
+		rterror("rt_raster_gdal_warp: Unable to set geotransform");
+
+		GDALClose(dst_ds);
+		GDALClose(src_ds);
+
 		for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
 		rtdealloc(transform_opts);
 
@@ -9929,36 +10006,6 @@ rt_raster rt_raster_gdal_warp(
 		}
 	}
 
-	/* set dst srs */
-	_dst_srs = rt_util_gdal_convert_sr((NULL == dst_srs ? src_srs : dst_srs), 1);
-	cplerr = GDALSetProjection(dst_ds, _dst_srs);
-	CPLFree(_dst_srs);
-	if (cplerr != CE_None) {
-		rterror("rt_raster_gdal_warp: Unable to set projection");
-
-		GDALClose(dst_ds);
-		GDALClose(src_ds);
-
-		for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
-		rtdealloc(transform_opts);
-
-		return NULL;
-	}
-
-	/* set dst geotransform */
-	cplerr = GDALSetGeoTransform(dst_ds, _gt);
-	if (cplerr != CE_None) {
-		rterror("rt_raster_gdal_warp: Unable to set geotransform");
-
-		GDALClose(dst_ds);
-		GDALClose(src_ds);
-
-		for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
-		rtdealloc(transform_opts);
-
-		return NULL;
-	}
-
 	/* create transformation object */
 	transform_arg = imgproj_arg = GDALCreateGenImgProjTransformer2(src_ds, dst_ds, transform_opts);
 	if (NULL == transform_arg) {
diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h
index 71e209e87..7c9fc4c48 100644
--- a/raster/rt_core/rt_api.h
+++ b/raster/rt_core/rt_api.h
@@ -1594,8 +1594,9 @@ rt_raster rt_raster_from_gdal_dataset(GDALDatasetH ds);
  *
  * @return the warped raster or NULL
  */
-rt_raster rt_raster_gdal_warp(rt_raster raster, const char *src_srs,
-	const char *dst_srs,
+rt_raster rt_raster_gdal_warp(
+	rt_raster raster,
+	const char *src_srs, const char *dst_srs,
 	double *scale_x, double *scale_y,
 	int *width, int *height,
 	double *ul_xw, double *ul_yw,
diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c
index 1c58ee43f..f0da14e22 100644
--- a/raster/rt_pg/rt_pg.c
+++ b/raster/rt_pg/rt_pg.c
@@ -11259,17 +11259,21 @@ Datum RASTER_GDALWarp(PG_FUNCTION_ARGS)
 		dst_srid = src_srid;
 	POSTGIS_RT_DEBUGF(4, "destination SRID: %d", dst_srid);
 
-	/* source SRID = SRID_UNKNOWN */
-	if (src_srid == SRID_UNKNOWN) {
-		/* target SRID != src SRID, error */
-		if (dst_srid != src_srid) {
-			elog(ERROR, "RASTER_GDALWarp: Input raster has unknown (%d) SRID", src_srid);
-			rt_raster_destroy(raster);
-			PG_FREE_IF_COPY(pgraster, 0);
-			PG_RETURN_NULL();
+	/* target SRID != src SRID, error */
+	if (src_srid == SRID_UNKNOWN && dst_srid != src_srid) {
+		elog(ERROR, "RASTER_GDALWarp: Input raster has unknown (%d) SRID", src_srid);
+		rt_raster_destroy(raster);
+		PG_FREE_IF_COPY(pgraster, 0);
+		PG_RETURN_NULL();
+	}
+	/* target SRID == src SRID, no reprojection */
+	else if (dst_srid == src_srid) {
+		/* set geotransform */
+		if (src_srid == SRID_UNKNOWN) {
+			double gt[6] = {0, 10, 0, 0, 0, -10};
+			rt_raster_set_geotransform_matrix(raster, gt);
 		}
 
-		/* target SRID == src SRID, special */
 		no_srid = 1;
 	}
 
@@ -11368,12 +11372,8 @@ Datum RASTER_GDALWarp(PG_FUNCTION_ARGS)
 	}
 
 	/* get srses from srids */
-	/* source srs */
-	/* no SRID, use EPSG:4326 (WGS84) */
-	if (no_srid)
-		src_srs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs";
-	/* legitimate SRID */
-	else {
+	if (!no_srid) {
+		/* source srs */
 		src_srs = rtpg_getSR(src_srid);
 		if (NULL == src_srs) {
 			elog(ERROR, "RASTER_GDALWarp: Input raster has unknown SRID (%d)", src_srid);
@@ -11381,19 +11381,12 @@ Datum RASTER_GDALWarp(PG_FUNCTION_ARGS)
 			PG_FREE_IF_COPY(pgraster, 0);
 			PG_RETURN_NULL();
 		}
-	}
-	POSTGIS_RT_DEBUGF(4, "src srs: %s", src_srs);
+		POSTGIS_RT_DEBUGF(4, "src srs: %s", src_srs);
 
-	/* target srs */
-	/* no SRID, use src_srs */
-	if (no_srid)
-		dst_srs = src_srs;
-	/* legitimate SRID */
-	else if (dst_srid != SRID_UNKNOWN) {
 		dst_srs = rtpg_getSR(dst_srid);
 		if (NULL == dst_srs) {
 			elog(ERROR, "RASTER_GDALWarp: Target SRID (%d) is unknown", dst_srid);
-			if (!no_srid && NULL != src_srs) pfree(src_srs);
+			if (!no_srid) pfree(src_srs);
 			rt_raster_destroy(raster);
 			PG_FREE_IF_COPY(pgraster, 0);
 			PG_RETURN_NULL();
@@ -11413,17 +11406,21 @@ Datum RASTER_GDALWarp(PG_FUNCTION_ARGS)
 	rt_raster_destroy(raster);
 	PG_FREE_IF_COPY(pgraster, 0);
 	if (!no_srid) {
-		if (NULL != src_srs) pfree(src_srs);
-		if (NULL != dst_srs) pfree(dst_srs);
+		pfree(src_srs);
+		pfree(dst_srs);
 	}
 	if (!rast) {
 		elog(ERROR, "RASTER_band: Could not create transformed raster");
 		PG_RETURN_NULL();
 	}
 
-	/* add target SRID but only if we're not using a no SRID */
-	if (!no_srid)
-		rt_raster_set_srid(rast, dst_srid);
+	/* add target SRID */
+	rt_raster_set_srid(rast, dst_srid);
+
+	if (no_srid && src_srid == SRID_UNKNOWN) {
+		double gt[6] = {0, 1, 0, 0, 0, -1};
+		rt_raster_set_geotransform_matrix(rast, gt);
+	}
 
 	pgrast = rt_raster_serialize(rast);
 	rt_raster_destroy(rast);
diff --git a/raster/test/regress/rt_gdalwarp.sql b/raster/test/regress/rt_gdalwarp.sql
index fb05c0244..76593ee06 100644
--- a/raster/test/regress/rt_gdalwarp.sql
+++ b/raster/test/regress/rt_gdalwarp.sql
@@ -247,7 +247,30 @@ INSERT INTO raster_gdalwarp_dst (rid, rast) VALUES (
 		ST_SetSRID(rast, 0),
 		'NearestNeighbor', 0.125,
 		NULL,
-		500., 500.
+		NULL, NULL,
+		NULL, NULL,
+		NULL, NULL,
+		5, 5
+	) FROM raster_gdalwarp_src)
+), (
+	0.26, (SELECT _st_gdalwarp(
+		ST_SetSRID(rast, 0),
+		'NearestNeighbor', 0.125,
+		NULL,
+		NULL, NULL,
+		NULL, NULL,
+		NULL, NULL,
+		2, 2
+	) FROM raster_gdalwarp_src)
+), (
+	0.27, (SELECT _st_gdalwarp(
+		ST_SetSRID(rast, 0),
+		'NearestNeighbor', 0.125,
+		NULL,
+		NULL, NULL,
+		NULL, NULL,
+		NULL, NULL,
+		100, 100
 	) FROM raster_gdalwarp_src)
 );
 
diff --git a/raster/test/regress/rt_gdalwarp_expected b/raster/test/regress/rt_gdalwarp_expected
index ca376cc32..9ce521b09 100644
--- a/raster/test/regress/rt_gdalwarp_expected
+++ b/raster/test/regress/rt_gdalwarp_expected
@@ -22,7 +22,9 @@ NOTICE:  Values must be provided for both X and Y when specifying the scale.  Re
 0.22|993310|26|26|1|500.000|500.000|0.000|6.000|950452.000|1396632.000|t|t|t
 0.23|984269|12|8|1|0.012|-0.012|0.000|0.000|-107.029|50.206|t|t|t
 0.24|974269|12|8|1|0.012|-0.012|0.000|0.000|-107.029|50.206|t|t|t
-0.25|0|20|20|1|500.000|500.000|0.000|0.000|-500000.000|590000.000|t|t|t
+0.25|0|5|5|1|1.000|-1.000|0.000|0.000|0.000|0.000|t|t|t
+0.26|0|2|2|1|1.000|-1.000|0.000|0.000|0.000|0.000|t|t|t
+0.27|0|100|100|1|1.000|-1.000|0.000|0.000|0.000|0.000|t|t|t
 0.3|994269|12|8|1|0.012|-0.012|0.000|0.000|-107.029|50.206|t|t|t
 0.4|993310|24|24|1|500.000|500.000|0.000|0.000|950732.188|1397281.783|t|t|t
 0.5|992163|10|10|1|1000.000|-1000.000|0.000|0.000|-500000.000|600000.000|t|t|t