4

I am trying to georeference an image using many (>1,000,000) GCPs (derived by matching two images using triangulation). However, after many attempts doing this in different ways I think I found that there is a limit to the number of GCPs that I can use.

ds = gdal.Open("tmp_GCP.TIF", gdal.GA_Update)
wkt = ds.GetProjection()
ds.SetGCPs(gcp_list_gdal[:10923], wkt)
print(ds.GetGCPCount()) # output 10923
ds.FlushCache()
ds = None
gdal.Open("tmp_GCP.TIF").GetGCPCount() # output 0

When I set 10923 GCPs I obtain a file that seemingly has zero GCPs and when I set this to 10922 I do get a file with 10922 GCPs. This seems to be independent of the width/height of the raster.

Also, I tried to set the GDAL cache to a higher value, but this also didn't have an effect.

gdal.SetConfigOption('GDAL_CACHEMAX', '1000000')

Edit: Just found out that I am getting the following error in the terminal when running this code in a Jupyter Notebook (which I don't get when I use gcp_list_gdal[:10922]).

Warning 1: TIFFFetchNormalTag:Incorrect count for "GeoTiePoints"; tag ignored
Warning 1: TIFFFetchNormalTag:Incorrect count for "GeoTiePoints"; tag ignored

I have tried this on Ubuntu 16 and Ubuntu 18, using Gdal 2.2, 2.4 and 3.0, but I always get the same error.

Edit 2: Because of the error I thought it might be a problem with the amount of GCP information that can be stored in a GeoTIFF file, so I tried another approach in python (trying not to store the DS before applying a warping), but it seems the same problem occurs (since it does produce an output using 10922 GCPs but not when using 10923 GCPs).

ds = gdal.Open('tmp_GCP.TIF', gdal.GA_Update) 
gcp_srs = osr.SpatialReference()
gcp_srs.ImportFromEPSG(4326)
wkt = gcp_srs.ExportToWkt()
ds.SetGCPs(gcp_list_gdal[:10922], wkt)
dstDS = gdal.Warp('TEST.TIF', ds, format='GTiff', tps=True, dstNodata=0,
                  srcNodata=0, resampleAlg=gdal.GRA_Bilinear)

Edit 3: I found out about the option to use geolocation arrays https://gdal.org/development/rfc/rfc4_geolocate.html which seems very interesting. As I understand it you can basically create two bands of the same size as the image and indicate for each pixel the desired longitude and latitude of the pixel. This should work better with a large number of GCPs so I am going to try this.

3
  • 1
    Out of curiosity, why do you need >1,000,000 GCPs? Nov 29 '19 at 11:18
  • We want our image to match as much as possible with our calculated height map (calculated using stereo processing).
    – David
    Nov 29 '19 at 11:46
  • This is actually related to a previous question I posted about not being able to orthorectify my WorldView-3 images using the RPC model (gis.stackexchange.com/questions/328366/…). So that's why I am trying to use GCP points directly.
    – David
    Nov 29 '19 at 12:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.