Re: [gdal-dev] Warping with GCPs at high latitudes

2013-01-25 Thread Knut-Frode Dagestad

On 24. jan. 2013 20:43, Even Rouault wrote:

Did you try to do your suggestion at hand to actually check that your idea
leads to the expected result ? This is basically a matter of running :

1) gdaltransform -s_srs EPSG:4326 -t_srs +proj=stere +lat_0=90 +lon_0=0 on
the coordinates of the GCPs of the source image
2) gdal_translate -a_srs  +proj=stere +lat_0=90 +lon_0=0  [-gcp x y X Y]*
src.tif src_reprojected_gcp.tif where x y come from the source image and X Y
are the reprojected coordinates
3) gdalwarp src_reprojected_gcp.tif out.tif


I tried it now with Python API, and it works perfectly, all the way to 
the pole.


Thus the conclusion should be that GCPs can have any spatial reference 
system, but preferably one with no singularity within or near the 
destination area.


I used code like this:
# Make tranformer from GCP SRS to destination SRS
dstSRS = osr.SpatialReference()
dstSRS.ImportFromProj4(srsString)
srcSRS = osr.SpatialReference()
srcGCPProjection = srcDataset.GetGCPProjection()
srcSRS.ImportFromWkt(srcGCPProjection)
transformer = osr.CoordinateTransformation(srcSRS, dstSRS)

# Reproject all GCPs
srcGCPs = self.vrt.dataset.GetGCPs()
dstGCPs = []
for srcGCP in srcGCPs:
	(x, y, z) = transformer.TransformPoint(srcGCP.GCPX, srcGCP.GCPY, 
srcGCP.GCPZ)
	dstGCP = gdal.GCP(x, y, z, srcGCP.GCPPixel, srcGCP.GCPLine, 
srcGCP.Info, srcGCP.Id)

dstGCPs.append(dstGCP)

# Update dataset
self.vrt.dataset.SetGCPs(dstGCPs, dstSRS.ExportToWkt())

___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev


[gdal-dev] Warping with GCPs at high latitudes

2013-01-24 Thread Knut-Frode Dagestad

Hi,

I am having trouble to perform accurate warping of images with GCPs 
(lon,lat) in the high Arctic (to +proj=stere +lat_0=90 +lon_0=0); 
problems starting already around 80 degrees of latitude.
Thin Plate Spline (-tps) performs better than polynomial fits, but 
output is often severely dragged at one or two edges, and sometimes it 
looks like the destination image is mirrored across a line.
Trimming images at the edges and adjusting warp options SAMPLE_STEPS and 
SAMPLE_GRID helps a little, but not sufficiently.


For the polynomial fits I assume the problem is related to precision due 
to the fact that lon and lat become singular at the pole. However, it 
should then probably be better if the algorithm could:
- convert first GCPs to the output projection (in my case stereographic, 
which has no singularity for the relevant images)

- make the polynomial fit GCP_stereo - pixel, line
- use this polynomial for the warping

I guess most GDAL users have a long wishlist of switches, but how about 
a switch to control whether GCP polynomial should rather be made in 
destination projection? Or a switch to gdal_translate to change the 
GCPProjection and GCP coordinates. A simpler and more pragmatic solution 
could be to make a separate Python script gdal_translate_gcps.py?


Are there other suggestion to deal with these problems?


Best regards from Knut-Frode

___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev


Re: [gdal-dev] Warping with GCPs at high latitudes

2013-01-24 Thread Even Rouault
Le jeudi 24 janvier 2013 13:18:50, Knut-Frode Dagestad a écrit :
 Hi,
 
 I am having trouble to perform accurate warping of images with GCPs
 (lon,lat) in the high Arctic (to +proj=stere +lat_0=90 +lon_0=0);
 problems starting already around 80 degrees of latitude.
 Thin Plate Spline (-tps) performs better than polynomial fits, but
 output is often severely dragged at one or two edges, and sometimes it
 looks like the destination image is mirrored across a line.
 Trimming images at the edges and adjusting warp options SAMPLE_STEPS and
 SAMPLE_GRID helps a little, but not sufficiently.
 
 For the polynomial fits I assume the problem is related to precision due
 to the fact that lon and lat become singular at the pole. However, it
 should then probably be better if the algorithm could:
 - convert first GCPs to the output projection (in my case stereographic,
 which has no singularity for the relevant images)
 - make the polynomial fit GCP_stereo - pixel, line
 - use this polynomial for the warping
 
 I guess most GDAL users have a long wishlist of switches, but how about
 a switch to control whether GCP polynomial should rather be made in
 destination projection? Or a switch to gdal_translate to change the
 GCPProjection and GCP coordinates. A simpler and more pragmatic solution
 could be to make a separate Python script gdal_translate_gcps.py?
 
 Are there other suggestion to deal with these problems?

In your use case, the source image is more or less very close to its final 
shape in Stereographics projection (the warping doesn't warp much the image), 
right ?

Did you try to do your suggestion at hand to actually check that your idea 
leads to the expected result ? This is basically a matter of running :

1) gdaltransform -s_srs EPSG:4326 -t_srs +proj=stere +lat_0=90 +lon_0=0 on 
the coordinates of the GCPs of the source image
2) gdal_translate -a_srs  +proj=stere +lat_0=90 +lon_0=0  [-gcp x y X Y]* 
src.tif src_reprojected_gcp.tif where x y come from the source image and X Y 
are the reprojected coordinates
3) gdalwarp src_reprojected_gcp.tif out.tif

 
 
 Best regards from Knut-Frode
 
 ___
 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