[gdal-dev] Re: [Fwd: Re: Warping an already georeferenced image with control points
Jan Hartmann uva.nl> writes: > > No, this doesn't seem to work. Although the coordinate system *is* > retained when I add -a_srs to the gcps, the boundaries of the input file > are discarded: they are reset to the pixel dimensions of the file. Is > this a bug? > > Jan > > Frank Warmerdam wrote: > > Jan Hartmann wrote: > >> Hi, > >> > >> I am trying to georeference old maps in two steps: first > >> georeferencing them from the old projection to a modern one, and then > >> ameliorating the georeferenced image by rubber-sheeting it on the > >> base of control points, mostly triangulation points for which the > >> coordinates in the old and new maps are known. > >> > >> The second step causes trouble: if I add control points to the > >> georeferenced image with gdal_translate -gcp , the projection > >> information gets lost. Hi, I am not sure if I understood how you trying to do the job, but hopefully my task and work flow are close enough to be useful for you. 1. Let's say that the original image is in epsg:2393 and target is epsg:3067. 2. I insert some ground control points in epsg:3067 with gdal_translate and -gcp. Image format of the temporal output file can be GeoTIFF but also virtual raster(vrt). 3. I warp the temporary image, now containing gpc's in epsg:3067 projection, to final product using gdalwarp as gdalwarp -s_srs epsg:3067 -t_srs epsg:3067. -Jukka Rahkonen- ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: [Fwd: Re: Warping an already georeferenced image with control points
Jan Hartmann 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@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: [Fwd: Re: Warping an already georeferenced image with control points
Jukka Rahkonen wrote: Jan Hartmann uva.nl> writes: No, this doesn't seem to work. Although the coordinate system *is* retained when I add -a_srs to the gcps, the boundaries of the input file are discarded: they are reset to the pixel dimensions of the file. Is this a bug? Jan Frank Warmerdam wrote: Jan Hartmann wrote: Hi, I am trying to georeference old maps in two steps: first georeferencing them from the old projection to a modern one, and then ameliorating the georeferenced image by rubber-sheeting it on the base of control points, mostly triangulation points for which the coordinates in the old and new maps are known. The second step causes trouble: if I add control points to the georeferenced image with gdal_translate -gcp , the projection information gets lost. Hi, I am not sure if I understood how you trying to do the job, but hopefully my task and work flow are close enough to be useful for you. 1. Let's say that the original image is in epsg:2393 and target is epsg:3067. 2. I insert some ground control points in epsg:3067 with gdal_translate and -gcp. Image format of the temporal output file can be GeoTIFF but also virtual raster(vrt). 3. I warp the temporary image, now containing gpc's in epsg:3067 projection, to final product using gdalwarp as gdalwarp -s_srs epsg:3067 -t_srs epsg:3067. -Jukka Rahkonen- ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev Hi Jukka, You have described exactly what I want to do, but the problem is that in step 2, inserting gcps, the georeferenced boundaries of the file are lost. Eventually you get a somewhat displaced raw raster file. Jan ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: [Fwd: Re: Warping an already georeferenced image with control points
Jan Hartmann uva.nl> writes: Jukka Rahkonen wrote: >> Hi, >> >> I am not sure if I understood how you trying to do the job, but hopefully my >> task and work flow are close enough to be useful for you. >> >> 1. Let's say that the original image is in epsg:2393 and target is epsg:3067. >> 2. I insert some ground control points in epsg:3067 with gdal_translate and >> -gcp. Image format of the temporal output file can be GeoTIFF but also virtual >> raster(vrt). >> 3. I warp the temporary image, now containing gpc's in epsg:3067 projection, to >> final product using gdalwarp as >> gdalwarp -s_srs epsg:3067 -t_srs epsg:3067. > > >> Hi Jukka, >> You have described exactly what I want to do, but the problem is that >> in step 2, inserting gcps, the georeferenced boundaries of the file are >> lost. Eventually you get a somewhat displaced raw raster file. >> Jan Hi, That displacement could be avoided at least by inserting manually gcp to each corner of the image. Perhaps this article gives some ideas: http://www.scangis.org/scangis2007/papers/r3_rahkonen.pdf -Jukka- ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: [Fwd: Re: Warping an already georeferenced image with control points
Thanks Jukka, that is an interesting article in its own right, and it is relevant to historical maps too. Your solution is viable, but it takes more work (insert corner gcp in each image), which shouldn't be necessary, as this information is already available but is discarded by gdal_translate. If gdal_translate could keep this information in the target file, even when gcps have been added, everything would be fine (and my impression is it should). Jan Rahkonen Jukka wrote: Jan Hartmann uva.nl> writes: Jukka Rahkonen wrote: >> Hi, >> >> I am not sure if I understood how you trying to do the job, but hopefully my >> task and work flow are close enough to be useful for you. >> >> 1. Let's say that the original image is in epsg:2393 and target is epsg:3067. >> 2. I insert some ground control points in epsg:3067 with gdal_translate and >> -gcp. Image format of the temporal output file can be GeoTIFF but also virtual >> raster(vrt). >> 3. I warp the temporary image, now containing gpc's in epsg:3067 projection, to >> final product using gdalwarp as >> gdalwarp -s_srs epsg:3067 -t_srs epsg:3067. > > >> Hi Jukka, >> You have described exactly what I want to do, but the problem is that >> in step 2, inserting gcps, the georeferenced boundaries of the file are >> lost. Eventually you get a somewhat displaced raw raster file. >> Jan Hi, That displacement could be avoided at least by inserting manually gcp to each corner of the image. Perhaps this article gives some ideas: http://www.scangis.org/scangis2007/papers/r3_rahkonen.pdf -Jukka- ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: [Fwd: Re: Warping an already georeferenced image with control points
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 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@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
Index: apps/gdal_translate.cpp
===
--- apps/gdal_translate.cpp (revision 16283)
+++ apps/gdal_translate.cpp (working copy)
@@ -57,7 +57,7 @@
" [-scale [src_min src_max [dst_min dst_max]]]\n"
" [-srcwin xoff yoff xsize ysize] [-projwin ulx uly lrx lry]\n"
" [-a_srs srs_def] [-a_ullr ulx uly lrx lry] [-a_nodata value]\n"
-" [-gcp pixel line easting northing [elevation]]*\n"
+" [-gcp pixel line easting northing [elevation]]* [-add_corners_as_gcps]\n"
" [-mo \"META-TAG=VALUE\"]* [-quiet] [-sds]\n"
" [-co \"NAME=VALUE\"]*\n"
" src_dataset dst_dataset\n\n" );
@@ -115,6 +115,7 @@
int bSetNoData = FALSE;
double dfNoDataReal = 0.0;
int nRGBExpand = 0;
+int bAddCornersAsGCPs = FALSE;
anSrcWin[0] = 0;
@@ -206,6 +207,9 @@
else if( EQUAL(argv[i],"-sds") )
bCopySubDatasets = TRUE;
+
+else if( EQUAL(argv[i],"-add_corners_as_gcps") )
+bAddCornersAsGCPs = TRUE;
else if( EQUAL(argv[i],"-gcp") && i < argc - 4 )
{
@@ -489,6 +493,59 @@
bDefBands = FALSE;
}
+if (bAddCornersAsGCPs)
+{
+double adfGeoTransform[6];
+GDALGetGeoTransform( hDataset, adfGeoTransform );
+
+if (pszOutputSRS == NULL && GDALGetProjectionRef(hDataset) != NULL)
+pszOutputSRS = CPLStrdup(GDALGetProjectionRef( hDataset ));
+
+pasGCPs = (GDAL_GCP *)
+CPLRealloc( pasGCPs, sizeof(GDAL_GCP) * (nGCPCount + 4) );
+GDALInitGCPs( 4, pasGCPs + nGCPCount );
+
+nGCPCount ++;
+pasGCPs[nGCPCount-1].dfGCPPixel = 0;
+pasGCPs[nGCPCount-1].dfGCPLine = 0;
+pasGCPs[nGCPCount-1].dfGCPX = adfGeoTransform[0] +
+ pasGCPs[nGCPCount-1].dfGCPPixel * adfGeoTransform[1] +
+ pasGCPs[nGCPCount-1].dfGCPLine * adfGeoTransform[2];
+pasGCPs[nGCPCount-1].dfGCPY = adfGeoTransform[3] +
+ pasGCPs[nGCPCount-1].dfGCPPixel * adfGeoTransform[4] +
+ pasGCPs[nGCPCount-1].dfGCPLine * adfGeoTransform[5];
+
+nGCPCount ++;
+pasGCPs[nGCPCount-1].dfGCPPixel = nRasterXSize;
+pasGCPs[nGCPCount-1].dfGCPLine = 0;
+pasGCPs[nGCPCount-1].dfGCPX = adfGeoTransform[0] +
+ pasGCPs[nGCP
Re: [gdal-dev] Re: [Fwd: Re: Warping an already georeferenced image with control points
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 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@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
Re: [gdal-dev] Re: [Fwd: Re: Warping an already georeferenced image with control points
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 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@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
Re: [gdal-dev] Re: [Fwd: Re: Warping an already georeferenced image with control points
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 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
Re: [gdal-dev] Re: [Fwd: Re: Warping an already georeferenced image with control points
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",50], PARAMETER["false_northing",0], UNIT["metre",1, AUTHORITY["EPSG","9001"]], AUTHORITY["EPSG","26912"]] Origin = (275804.27707260562,5462603.41083706920) 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 mosai
Re: [gdal-dev] Re: [Fwd: Re: Warping an already georeferenced image with control points
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 '&& 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 (http://www.remotesensing.org/geotiff/spec/geotiff2.6.html#2.6.1) 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 : > 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_082h0 >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",50], > PARAMETER["false_northing",0], > UNIT["metre",1, > AUTHORITY["EPSG","9001"]], > AUTHORITY["EPSG","26912"]] > Origin = (275804.27707260562,5462603.41083706920) > 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,
Re: [gdal-dev] Re: [Fwd: Re: Warping an already georeferenced image with control points
Yes, I think I understand the problem now: a GeoTiff header cannot have *both* SRS *and* GCP information. Even so, I think there is a point to be made for a gcp-transformation of an already georeferenced image. Would adding gcp's to the gdalwarp program be a solution? Jan Even Rouault wrote: 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 '&& 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 (http://www.remotesensing.org/geotiff/spec/geotiff2.6.html#2.6.1) 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 : 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_082h0 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",50], PARAMETER["false_northing",0], UNIT["metre",1, AUTHORITY["EPSG","9001"]], AUTHORITY["EPSG","26912"]] Origin = (275804.27707260562,5462603.41083706920) 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 geo
