Jan,
I concur with Jukka's analysis. Your need is rather uncommon, and I also
thinks it shouldn't be a default behaviour. However, adding that capability
to gdal_translate is quite easy. I've attached a patch for gdal_translate.cpp
that adds the capability of adding the 4 corners as 4 additional GCPs with
the "-add_corners_as_gcps" option. You might give it at try (it applies
cleanly on GDAL 1.6.0 and trunk). I'm not sure yet if it is valuable enough
to include it in baseline gdal_translate.cpp.
Best regards
Le Wednesday 11 February 2009 15:59:48 Jukka Rahkonen, vous avez écrit :
> Jan Hartmann <j.l.h.hartmann <at> uva.nl> writes:
> > Sorry to keep moaning about this, but I need an indication what's going
> > on here. Mind, I don't need an immediate solution: for the time being I
> > have a workaround. Just an idea whether this a real problem, a dumb
> > question, something that can be handled in the foreseeable future (or
> > not), perhaps with adequate funding. Everything is better than talking
> > to a blind wall.
> >
> > Sorry again Frank,
> >
> > Jan
>
> Hi,
>
> I think that your need to keep the original extents but add more ground
> control points inside the image frame is rather uncommon. Doesn't it mean
> that you trust that the image corners are correctly warped to a new
> projection, but there are local distortions in the middle of image which
> should be corrected with a few extra ground control points? For my mind it
> shoudn't be the default behaviour of gdal but it might be usable as an user
> selectable option sometimes.
>
> I know that missing gcp's at the image corners often leads to very odd
> result with polynomial warping because the formula shoots over. Even
> unaccurately placed gcp's could help a lot in preserving the original shape
> of the map. I guess that you are perhaps playing with scanned historical
> maps? I have a few old scanned parcel maps which are covering just the
> area of the farm, and two of the map corners are just white background. It
> is impossible to measure any real ground control points from the corners
> because there is nothing on the map to compare with, and warping with all
> gcp's on the middle area makes really funny looking curves into the hand
> drawn rectangle framing the original map.
>
> -Jukka-
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
Index: apps/gdal_translate.cpp
===================================================================
--- apps/gdal_translate.cpp (revision 16283)
+++ apps/gdal_translate.cpp (working copy)
@@ -57,7 +57,7 @@
" [-scale [src_min src_max [dst_min dst_max]]]\n"
" [-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 pixel line easting northing [elevation]]* [-add_corners_as_gcps]\n"
" [-mo \"META-TAG=VALUE\"]* [-quiet] [-sds]\n"
" [-co \"NAME=VALUE\"]*\n"
" src_dataset dst_dataset\n\n" );
@@ -115,6 +115,7 @@
int bSetNoData = FALSE;
double dfNoDataReal = 0.0;
int nRGBExpand = 0;
+ int bAddCornersAsGCPs = FALSE;
anSrcWin[0] = 0;
@@ -206,6 +207,9 @@
else if( EQUAL(argv[i],"-sds") )
bCopySubDatasets = TRUE;
+
+ else if( EQUAL(argv[i],"-add_corners_as_gcps") )
+ bAddCornersAsGCPs = TRUE;
else if( EQUAL(argv[i],"-gcp") && i < argc - 4 )
{
@@ -489,6 +493,59 @@
bDefBands = FALSE;
}
+ if (bAddCornersAsGCPs)
+ {
+ double adfGeoTransform[6];
+ GDALGetGeoTransform( hDataset, adfGeoTransform );
+
+ if (pszOutputSRS == NULL && GDALGetProjectionRef(hDataset) != NULL)
+ pszOutputSRS = CPLStrdup(GDALGetProjectionRef( hDataset ));
+
+ pasGCPs = (GDAL_GCP *)
+ CPLRealloc( pasGCPs, sizeof(GDAL_GCP) * (nGCPCount + 4) );
+ GDALInitGCPs( 4, pasGCPs + nGCPCount );
+
+ nGCPCount ++;
+ pasGCPs[nGCPCount-1].dfGCPPixel = 0;
+ pasGCPs[nGCPCount-1].dfGCPLine = 0;
+ pasGCPs[nGCPCount-1].dfGCPX = adfGeoTransform[0] +
+ pasGCPs[nGCPCount-1].dfGCPPixel * adfGeoTransform[1] +
+ pasGCPs[nGCPCount-1].dfGCPLine * adfGeoTransform[2];
+ pasGCPs[nGCPCount-1].dfGCPY = adfGeoTransform[3] +
+ pasGCPs[nGCPCount-1].dfGCPPixel * adfGeoTransform[4] +
+ pasGCPs[nGCPCount-1].dfGCPLine * adfGeoTransform[5];
+
+ nGCPCount ++;
+ pasGCPs[nGCPCount-1].dfGCPPixel = nRasterXSize;
+ pasGCPs[nGCPCount-1].dfGCPLine = 0;
+ pasGCPs[nGCPCount-1].dfGCPX = adfGeoTransform[0] +
+ pasGCPs[nGCPCount-1].dfGCPPixel * adfGeoTransform[1] +
+ pasGCPs[nGCPCount-1].dfGCPLine * adfGeoTransform[2];
+ pasGCPs[nGCPCount-1].dfGCPY = adfGeoTransform[3] +
+ pasGCPs[nGCPCount-1].dfGCPPixel * adfGeoTransform[4] +
+ pasGCPs[nGCPCount-1].dfGCPLine * adfGeoTransform[5];
+
+ nGCPCount ++;
+ pasGCPs[nGCPCount-1].dfGCPPixel = nRasterXSize;
+ pasGCPs[nGCPCount-1].dfGCPLine = nRasterYSize;
+ pasGCPs[nGCPCount-1].dfGCPX = adfGeoTransform[0] +
+ pasGCPs[nGCPCount-1].dfGCPPixel * adfGeoTransform[1] +
+ pasGCPs[nGCPCount-1].dfGCPLine * adfGeoTransform[2];
+ pasGCPs[nGCPCount-1].dfGCPY = adfGeoTransform[3] +
+ pasGCPs[nGCPCount-1].dfGCPPixel * adfGeoTransform[4] +
+ pasGCPs[nGCPCount-1].dfGCPLine * adfGeoTransform[5];
+
+ nGCPCount ++;
+ pasGCPs[nGCPCount-1].dfGCPPixel = 0;
+ pasGCPs[nGCPCount-1].dfGCPLine = nRasterYSize;
+ pasGCPs[nGCPCount-1].dfGCPX = adfGeoTransform[0] +
+ pasGCPs[nGCPCount-1].dfGCPPixel * adfGeoTransform[1] +
+ pasGCPs[nGCPCount-1].dfGCPLine * adfGeoTransform[2];
+ pasGCPs[nGCPCount-1].dfGCPY = adfGeoTransform[3] +
+ pasGCPs[nGCPCount-1].dfGCPPixel * adfGeoTransform[4] +
+ pasGCPs[nGCPCount-1].dfGCPLine * adfGeoTransform[5];
+ }
+
/* -------------------------------------------------------------------- */
/* Compute the source window from the projected source window */
/* if the projected coordinates were provided. Note that the */
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
