Got the @#$% problem. Nothing to do with GDAL, it's the way QGIS handles
WMS requests. I wrote a small PHP mapscript to retrieve raw raster files
as a WMS service, to be able to view them in QGIS and to sample control
points. For some reason I put a 100 pixel margin in the extents of the
GetCapabilities request. QGIS caches those extents and uses them with
each following GetMap request, disregarding the actual extents of the
raster file. Not sure if it should do this, but even if not, it was a
bad idea to create that margin.
Anyway, I'm glad I got that devil's error. It's treacheous: the
displacements are very small and don't look like regular rotations or
translations (when you use thin plate warping). I've been showing those
georeferenced map overlaid over Google satellite to people over here in
Amsterdam, and they immediately noticed the displacements in their back
gardens. Being Dutch they then started talking enthusiastically about
suing the government over their land tax assessment. Can you imagine the
horror when "tout Amsterdam" would be disputing their taxes on the basis
on this magnificent piece of Open Source software?
Mind, I'm reasonably sure there are lots of errors in the present
digitized cadastral maps compared to their historical sources, as soon
you are going to look on meter-level accuracy. See for comparable
problems in Great Brittain
http://www.ordnancesurvey.co.uk/oswebsite/pai/ (thanks, Peter Hall, for
the link and other comments). However, before we open that can of worms
over here, I want to be absolutely sure that at least my own technology
is right, preferably at sub-pixel level, as it now is.
Thanks too Jukka and Even for confirming so quickly that the problem
couldn't be with gdalwarp. If you hadn't done so, I probably would be
still debugging gdal source code. Fun, of course, but not what I'm
supposed to do.
Cheers,
Jan
On 05/04/10 21:09, Even Rouault wrote:
Jan,
you'll need to provide the exact source data and command line you use because
I can't reproduce any problem.
I've generated a 1000x1000 artificial dataset with the following features :
* background is black
* a white grid made of horizontal lines and vertical lines spaced with 100
pixels
* 5 GCPs :
- red square at (150, 150). Target coords are (500000 + 100, 4500000 -
100)
- green square at (250,450). Target coords are (500000 + 100, 4500000 -
500)
- yellow square at (50, 850). Target coords are (500000 + 100, 4500000
- 900)
- cyan square at (950, 50). Target coords are (500000 + 900, 4500000 -
100)
- magenta square at (850, 950). Target coords are (500000 + 900,
4500000 -
900)
So the expected result is a warped grid, but such as the red,yellow,cyan and
magenta points form a perfect square (and the green square is in the middle
of the red-yellow segment). And that's exactly what I get. With a precision
of sub-pixel.
Attached a ZIP with the source image and C source code used to generate it.
Command line used to warp : gdalwarp -tps src_grid.tif dst_grid.tif
Best regards,
Even
Le Tuesday 04 May 2010 15:14:34 Jan Hartmann, vous avez écrit :
Thanks, Jukka, that's exactly what I did. I used QGIS to test the
coordinates in both images.
Jan
On 05/04/10 15:08, Jukka Rahkonen wrote:
Jan Hartmann<j.l.h.hartmann<at> uva.nl> writes:
you can see a screenshot of the original image (right) and the
georeferenced one (left, rotated almost 90 degrees). The red markers
with the number 3 inside them show one of the control points: to the
right the scan pixel (5023/3421), to the left the targeted
georeferenced coordinate in EPSG:28992 (121527/487174). As you see, the
georeferenced point is not exactly on the border of parcel 7, as I
expected. Above the map you see the four coordinates I used as control
points; I added four extra points half way the original ones, so the
complete warp was done with 8 control points.
Hi,
I can try to repeat your test some day. Is this your workflow:
1. gdal_translate -of VRT -gcp p1 l1 e1 n1 -gcp p2 l2 e2 n2 -gcp p3 l3 e3
n3 -gcp p4 l4 e4 n4 -gcp p5 l5 e5 n5 input.tif temp_with_gcp.vrt
2. gdalwarp -s_srs EPSG:28992 -t_srs EPSG:28992 -tps temp_with_gcp.vrt
warped.tif
3. Measure the coordinates of the ground control points from warped.tif
with some GIS software like QGis and find out disappointed that
they have shifted.
-Jukka Rahkonen-
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
