[gdal-dev] Re: [Fwd: Re: Warping an already georeferenced image with control points

Brent Fraser bfraser at geoanalytic.com
Fri Feb 13 12:52:29 EST 2009


I stepped through your procedure with a Canadian topo map to slightly adjust its positioning, using GDAL 1.6.0:

1. Warp raster image of topo map to a projected coordinate system
No Need.  Canadian topos already are georeferenced so I simply downloaded http://ftp2.cits.rncan.gc.ca/pub/canmatrix/50k_300dpi/082/h/canmatrix_082h04_tif.zip
2. Attach my new GCPs to the original topo map:
gdal_translate --optfile gcps.opt -a_srs "EPSG:26912" 082h04_1_1.tif 082h04_gcps.tif
	(the GCPs are in UTM zone 12, as indicated by the -a_srs parameter)

3. Warp/resample image using GCPs
gdalwarp -tps 082h04_gcps.tif 082h04_warped.tif

4. Examine SRS of output image:
gdalinfo  082h04_warped.tif

Driver: GTiff/GeoTIFF
Files: 082h04_warped.tif
Size is 11786, 8881
Coordinate System is:
PROJCS["NAD83 / UTM zone 12N",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
Origin = (275804.277072605620000,5462603.410837069200000)
Pixel Size = (4.233604626074191,-4.233604626074191)
Image Structure Metadata:
Corner Coordinates:
Upper Left  (  275804.277, 5462603.411) (114d 4'56.70"W, 49d16'30.13"N)
Lower Left  (  275804.277, 5425004.768) (114d 3'41.57"W, 48d56'14.30"N)
Upper Right (  325701.541, 5462603.411) (113d23'49.60"W, 49d17'28.67"N)
Lower Right (  325701.541, 5425004.768) (113d22'51.13"W, 48d57'12.15"N)
Center      (  300752.909, 5443804.089) (113d43'49.87"W, 49d 6'53.14"N)
Band 1 Block=11786x1 Type=Byte, ColorInterp=Palette
  Color Table (RGB with 256 entries)

The result of step 2 (gdal_translate) DOES remove the image's original georeferencing.  Actually, it replaces it with a set of GCPs (and their SRS).  To leave the original would resulting in two sets of conflicting georeferencing information.  

And to me it doesn't matter anyway, as that "GCP-ed" image is just temporary.  After the gdalwarp step I get an image with more traditional/well-supported georeferencing, I can delete it.

Should gdalinfo (and possibly other programs) use the GCPs to calculate the image extents?  Maybe.  I'd find it useful as a quality-control check to make sure I didn't mis-code a GCP resulting in a huge raster after warping.  Do you have another reason?


Jan Hartmann wrote:
> Yes exactly! Brent highlights a point I never mentioned, and that may be 
> the source of the confusion: when I add control points to an already 
> georeferenced image, I gdalwarp it with the -tps option: a thin spline 
> transformer, that transforms the control points to the exact 
> georeferenced position, but stretches the areas in between. AFAIK, this 
> is done by triangulation. This is different from a traditional warp, in 
> which the eventual position of the control points on the georeferenced 
> map has errors, although minimized ones. So in matching roads on two 
> adjacent images, the road edges on both images are assigned the same 
> georeferenced coordinates and then, after gdalwarp -tps, they exactly 
> align.
> So my workflow is:
> gdalwarp raw image --> georeferenced image with correct projection, 
> based on the corner points of the original map
> gdal_translate : add control points to projected image, essentially 
> small corrections of road edges and triangulation points for
>                         which we know the exact coordinates from later 
> sources
> gdalwarp -tps georeferenced image --> georeferenced image with small 
> corrections, where roads, churcht towers are on their exact locations
> The whole procedure is something I do regularly with ArcGIS, not only 
> with historical maps but also with things like aerial photographs, just 
> like Brent describes. It can be done with gdal too, as described above, 
> if only the second step would work. I still don't see why adding control 
> points to an image that already has a projection defined, should destroy 
> that projection information, including its georeferenced boundaries and 
> dimensions. Even when this is not the case: when  a -a_srs flag is 
> allowed for gdal_translate, it should be possible to add the complete 
> geolocation information for that image too: bounds, pixel dimensions and 
> geotranformation matrix, else the projection information doesn't make 
> sense.
> Have I made my intention a bit clearer now?
> Jan
> Brent Fraser wrote:
>> Jan,
>>  I think what you want gdalwarp to do is [Delauney] image 
>> triangulation.  Not an uncommon (or unreasonable) request.  I 
>> occasionally have a need for triangulation when mosaicking several 
>> satellite images.  Its visually pleasing to have the roads match up 
>> precisely where the images overlap...
>> Brent Fraser
>> Even Rouault wrote:
>>> Jan,
>>> (I'm CC'ing the list)
>>> I'm not sure what you mean with adapting the pixelsize, now that the 
>>> output has only GCPs and no more geotransform matrix.
>>> As far as including this option in baseline gdal_translate.cpp, I'm 
>>> currently not really enthousiastic, as there are lots of way to get 
>>> the result achieved (what is done with that patch could be done 
>>> similarly with GDAL python API for instance) and the choice of 
>>> corners is something rather arbitrary. Why should you trust the 
>>> georeferenced coordinates of the corner points as reliable GCPs ? Why 
>>> not some regularly discretized points along the edges ? Or a regular 
>>> grid over the whole raster ? etc. etc. So, I'd prefer hearing from 
>>> other GDAL developers or users that this patch is a sensical addition 
>>> before commiting it.
>>> Even
>>> Le Thursday 12 February 2009 11:47:00, vous avez écrit :
>>>> Hi Even,
>>>> Thanks a lo. Just one small addition: the pixelsize has also to be
>>>> adapted, else you get an image with false dimension. If it's not too
>>>> much to add that, I'm going to try it out and will let you know the
>>>> results. If it works, I have indeed a practical solution for my 
>>>> problem.
>>>> Even so, I am not convinced that this is an uncommon problem. 
>>>> Consider a
>>>> very basic GIS functionality: edge matching.  You have digitized a map
>>>> from many sheets, and at the end the georeferenced images don't match
>>>> precisely at the edges. Edge-matching functions (ArcGIS has lots of
>>>> them) ensure that the edge of a feature on one map perfectly matches 
>>>> the
>>>> feature edge on the adjacent map. This is exactly the kind of thing you
>>>> can do with gdalwarp using control points on a georeferenced image. It
>>>> is also more or less the thing I want to do with my historical maps.
>>>> Would this be an argument to preserve the boundaries of a georeferenced
>>>> image, when adding control points?
>>>> Jan
>>>> Even Rouault wrote:
>>>>> 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 at lists.osgeo.org
>>>>>> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>>> _______________________________________________
>>> gdal-dev mailing list
>>> gdal-dev at lists.osgeo.org
>>> http://lists.osgeo.org/mailman/listinfo/gdal-dev

More information about the gdal-dev mailing list