[gdal-dev] Re: Thin Plate Spline

[Jan Hartmann]

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 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
>>      
>
>