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

Jan Hartmann j.l.h.hartmann at uva.nl
Sun Feb 15 07:03:55 EST 2009


Yes, I think I understand the problem now: a GeoTiff header cannot have 
*both* SRS *and* GCP information. Even so, I think there is a point to 
be made for a gcp-transformation of an already georeferenced image. 
Would adding gcp's to the gdalwarp program be a solution?

Jan

Even Rouault wrote:
> Yes, once you have specified GCPs with the -gcp option, you remove the 
> original georeferencing made of the origin coordinates and pixel size.
>
> Several notes :
> - There is an exclusive test in gdal_translate.cpp that prevents you from 
> setting both both GCPs and (origin coordinates, pixel size) at the same time. 
> That can be easily disabled, by removing the '&& nGCPCount == 0 )' test at 
> line 692 of gdal_translate.cpp. (*)
> - ... but the GeoTIFF driver doesn't support setting both GCPs and (origin 
> coordinates, pixel size) at the same time. It wouldn't read both either. 
> Hacking the GTiff driver might be possible to make it support both. But my 
> understanding of the GeoTIFF specification 
> (http://www.remotesensing.org/geotiff/spec/geotiff2.6.html#2.6.1) is that it 
> is not a foreseen case. Although I don't see it is explicitely forbidden, I 
> would anticipate that it wouldn't be handled well by other software.
> - The VRT driver happens to support writing and reading both both GCPs and 
> (origin coordinates, pixel size) at the same time. But I doubt this is the 
> result of a design decision. It more looks like this works because that's 
> what required the minimum amount of coding !
> - Although I'm not 100% positive, I am not aware of any other GDAL driver that 
> would support both at the same time.
> - I don't think that you actually loose information by transforming the 
> (origin coordinates, pixel size) into 4 GCPs. With the raster dimension and 
> the coordinates, it's easy to retrieve the pixel size.
>
> Even
>
> (*) To tell the truth, by looking at the code, I've seen a workaround that 
> doesn't imply changing the code. Just use the -a_ullr option...
>
> Le Friday 13 February 2009 18:52:29 Brent Fraser, vous avez écrit :
>   
>> Jan,
>>
>> 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_082h0
>> 4_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",
>>     GEOGCS["NAD83",
>>         DATUM["North_American_Datum_1983",
>>             SPHEROID["GRS 1980",6378137,298.2572221010002,
>>                 AUTHORITY["EPSG","7019"]],
>>             AUTHORITY["EPSG","6269"]],
>>         PRIMEM["Greenwich",0],
>>         UNIT["degree",0.0174532925199433],
>>         AUTHORITY["EPSG","4269"]],
>>     PROJECTION["Transverse_Mercator"],
>>     PARAMETER["latitude_of_origin",0],
>>     PARAMETER["central_meridian",-111],
>>     PARAMETER["scale_factor",0.9996],
>>     PARAMETER["false_easting",500000],
>>     PARAMETER["false_northing",0],
>>     UNIT["metre",1,
>>         AUTHORITY["EPSG","9001"]],
>>     AUTHORITY["EPSG","26912"]]
>> Origin = (275804.277072605620000,5462603.410837069200000)
>> Pixel Size = (4.233604626074191,-4.233604626074191)
>> Metadata:
>>   AREA_OR_POINT=Area
>> Image Structure Metadata:
>>   INTERLEAVE=BAND
>> 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?
>>
>> Brent
>>
>> 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
>>>>>           
>> _______________________________________________
>> gdal-dev mailing list
>> gdal-dev at lists.osgeo.org
>> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>>     
>
>
>
>   
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20090215/fb282114/attachment.html


More information about the gdal-dev mailing list