[gdal-dev] gdal_translate: a preliminary `-gcp_file' patch

[Ivan Shmakov]

	Since the current version of GDAL doesn't support reading GCPs
	from some of the data sets I use [1], I have to use
	`gdal_translate' to specify all the GCPs, like:

$ file=/where/is/AIRS-L2-RetStd.hdf
$ gdal_translate \
      $(paste <(hdfdump --text "$file" {Long,Lat}itude) \
            | gawk '$1 != -9999 && $2 != -9999 { print "-gcp", (NR - 1) % 30 + .5, int ((NR - 1) / 30) + .5, $1, $2; }') \
      "HDF4_EOS:EOS_SWATH:\"$file\":L2_Standard_atmospheric&surface_product:TSurfStd" \
      AIRS-L2-RetStd.TSurfStd.gcp.tiff
$ 

	There, the dimensions of the data set are 30x45.

	However, would I try to use the above on a larger data set (say,
	MODIS 1 km, which is 1354x2030 for a standard NASA granule),
	I'll quickly run into the command line size limit (on most of
	the operating systems I know.)

	Therefore, I suggest the new `-gcp_from FILE' option for
	`gdal_translate'.

	The patch is roughly as follows.  Beware that the getline ()
	function is a GNU extension, available either from GNU Libc [2],
	or Gnulib [3].

--- gdal_translate.cpp.~4~	2007-12-13 11:32:31.000000000 +0600
+++ gdal_translate.cpp	2007-12-13 13:57:30.000000000 +0600
@@ -57,6 +57,7 @@
             "       [-srcwin xoff yoff xsize ysize] [-projwin ulx uly lrx lry]\n"
             "       [-a_srs srs_def] [-a_ullr ulx uly lrx lry] [-a_nodata value]\n"
             "       [-gcp pixel line easting northing [elevation]]*\n" 
+            "       [-gcp_from gcp_file]*\n" 
             "       [-mo \"META-TAG=VALUE\"]* [-quiet] [-sds]\n"
             "       [-co \"NAME=VALUE\"]*\n"
             "       src_dataset dst_dataset\n\n" );
@@ -79,6 +80,79 @@
 }
 
 /************************************************************************/
+/*                           handleGCPfrom ()                           */
+/************************************************************************/
+static int
+handleGCPfrom (const char *filename, GDAL_GCP **gcpp, int *ngcpp)
+{
+    /* NB: `ngcpp' should be of type `size_t *' rather than `int *' */
+    FILE *fp;
+    int lineno;
+    char *linebuf = 0;
+    size_t linesz = 0;
+    size_t new_points;
+
+    /* open the file, return -1 on failure */
+    fp = fopen (filename, "r");
+    if (fp == 0) {
+        /* use strerror () here? */
+        /* follow generic Unix ``program: message'' convention? */
+        fprintf (stderr, "WARNING: handleGCPfrom (): %s: %s\n",
+                 filename, "fopen () failed");
+        return -1;
+    }
+
+    /* parse the file */
+    /* NB: for better performance, CPLRealloc () the array in bigger
+     *     steps */
+    for (lineno = 0, new_points = 0; ! feof (fp); lineno++) {
+        int rv;
+        float pix, lin, x, y, elev;
+
+        /* NB: fscanf () is hardly a good idea */
+        {
+            ssize_t read;
+            if ((read = getline (&linebuf, &linesz, fp)) < 0) {
+                /* NB: a warning if not EOF here? */
+                break;
+            }
+        }
+
+        rv = sscanf (linebuf, " %f %f %f %f %f",
+                     &pix, &lin, &x, &y, &elev);
+        if (rv < 4) {
+            fprintf (stderr,
+                     ("WARNING: handleGCPfrom (): %s:%d"
+                      ": cannot parse line\n"),
+                     filename, lineno);
+            continue;
+        }
+
+        /* append the GCP just read to the array */
+
+        new_points++;
+        *gcpp = (GDAL_GCP *)
+            CPLRealloc (*gcpp, sizeof (**gcpp) * (++(*ngcpp)));
+        GDAL_GCP *tailp = *gcpp + *ngcpp - 1;
+        GDALInitGCPs (1, tailp);
+
+        tailp->dfGCPPixel = pix;
+        tailp->dfGCPLine = lin;
+        tailp->dfGCPX = x;
+        tailp->dfGCPY = y;
+    }
+
+    /* deallocate the line buffer */
+    CPLFree (linebuf);
+
+    /* close the file */
+    fclose (fp);
+
+    /* . */
+    return new_points;
+}
+
+/************************************************************************/
 /*                             ProxyMain()                              */
 /************************************************************************/
 
@@ -221,6 +295,11 @@
             /* should set id and info? */
         }   
 
+        else if (EQUAL (argv[i], "-gcp_from") && i < argc - 1) {
+            handleGCPfrom (argv[++i], &pasGCPs, &nGCPCount);
+            ++i;
+        }
+
         else if( EQUAL(argv[i],"-a_nodata") && i < argc - 1 )
         {
             bSetNoData = TRUE;

[1] http://trac.osgeo.org/gdal/ticket/2079
[2] http://www.gnu.org/software/libc/
[3] http://www.gnu.org/software/gnulib/