<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
  <meta content="text/html;charset=UTF-8" http-equiv="Content-Type">
</head>
<body bgcolor="#ffffff" text="#000000">
Thanks to Brent Fraser, I think I know now were the problem lies. He
even sent me some suggestions for a workaround, without even exactly
completely comprehending the problem (isn't that brilliant?). I want to
rubbersheet maps (gdalwarp -tps), using already georeferenced maps and
using control points in georeferenced coordinates. I do this to align
exactly historical maps, but it is also something that is useful in
edge-matching adjacent sheets of recent topographical maps. This
*cannot* be done by the gdal_translate / gdalwarp suite. The problem is
that control points have to be added to a tiff-file with
gdal_translate, and when that is done, the file has lost its
georeferencing information: SRS and dimensions. This is a restriction
of the geotiff format, so the problem is unsolvable.<br>
<br>
However, computationally everything is available in GDAL to implement
this functionality. With gdaltransform, I can do exactly what I want:
here, I translate a coordinate in the Dutch SRS (epsg:28992) one meter
east, based on four control points, also in epsg:28992:<br>
<br>
gdaltransform -s_srs epsg:28992 -t_srs epsg:28992                      
<br>
                       -gcp 179624 318062 179625 318062
<br>
                      -gcp 181371 318062 181372 318062
<br>
                      -gcp 181371 320751 181372 320751
<br>
                      -gcp 179624 320751 179625 320751 &lt;&lt;eof
<br>
 180000 320000
<br>
 180001 320000 0
<br>
eof<br>
<br>
So the problem is not that GDAL can't do it (GDAL can do practically
anything:-)), it is just that the Geotiff file format isn't up to the
problem.<br>
<br>
Now, this geotiff file with gcps is just an intermediate file. When
warping with control points, you need first to gdal_translate the
original image into a second file with gcps added, and than you can
gdalwarp that second file into the final map, a third file. I have
always wondered why this intermediate gcp-file was needed: you could
just as well transform directly from image to map, *if* gcps could be
indicated with gdalwarp, or am I overlooking something? Adding a -gcp
option to gdalwarp would AFAICS not be very difficult, it wouldn't
break existing applications, and it would make warping by control
points more efficient, because no intermediate files have to be
generated any more. I think I understand how the gcp option came to be
added to gdal_translate: some image providers, especially satellite
images, provide their images with control points within the file, so
the -gcp option has to be retained for gdal_translate.<br>
<br>
So, even if I have been unable to make clear my exotic problems with
historical maps, could this feature, adding a -gcp parameter to
gdalwarp, be implemented with a reasonable amount of programming effort?<br>
<br>
Thanks,<br>
<br>
Jan<br>
<br>
<br>
Even Rouault wrote:
<blockquote cite="mid:200902132002.44935.even.rouault@mines-paris.org"
 type="cite">
  <pre wrap="">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 '&amp;&amp; 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 
(<a class="moz-txt-link-freetext" href="http://www.remotesensing.org/geotiff/spec/geotiff2.6.html#2.6.1">http://www.remotesensing.org/geotiff/spec/geotiff2.6.html#2.6.1</a>) 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 :
  </pre>
  <blockquote type="cite">
    <pre wrap="">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
<a class="moz-txt-link-freetext" href="http://ftp2.cits.rncan.gc.ca/pub/canmatrix/50k_300dpi/082/h/canmatrix_082h0">http://ftp2.cits.rncan.gc.ca/pub/canmatrix/50k_300dpi/082/h/canmatrix_082h0</a>
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:
    </pre>
    <blockquote type="cite">
      <pre wrap="">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 --&gt; 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 --&gt; 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:
      </pre>
      <blockquote type="cite">
        <pre wrap="">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:
        </pre>
        <blockquote type="cite">
          <pre wrap="">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 :
          </pre>
          <blockquote type="cite">
            <pre wrap="">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:
            </pre>
            <blockquote type="cite">
              <pre wrap="">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 :
              </pre>
              <blockquote type="cite">
                <pre wrap="">Jan Hartmann &lt;j.l.h.hartmann &lt;at&gt; uva.nl&gt; writes:
                </pre>
                <blockquote type="cite">
                  <pre wrap="">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
                  </pre>
                </blockquote>
                <pre wrap="">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
<a class="moz-txt-link-abbreviated" href="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</a>
<a class="moz-txt-link-freetext" href="http://lists.osgeo.org/mailman/listinfo/gdal-dev">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a>
                </pre>
              </blockquote>
            </blockquote>
          </blockquote>
          <pre wrap="">_______________________________________________
gdal-dev mailing list
<a class="moz-txt-link-abbreviated" href="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</a>
<a class="moz-txt-link-freetext" href="http://lists.osgeo.org/mailman/listinfo/gdal-dev">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a>
          </pre>
        </blockquote>
      </blockquote>
    </blockquote>
    <pre wrap="">_______________________________________________
gdal-dev mailing list
<a class="moz-txt-link-abbreviated" href="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</a>
<a class="moz-txt-link-freetext" href="http://lists.osgeo.org/mailman/listinfo/gdal-dev">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a>
    </pre>
  </blockquote>
  <pre wrap=""><!---->


  </pre>
</blockquote>
</body>
</html>