[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
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_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",
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
>>>
>>
>
More information about the gdal-dev
mailing list