I am trying to rectify a tif of of sonar images collected from a moving ship. The tifs are in the form of long ribbons. I have ground control points (GCPs) that match coordinates on the tif (row and column of a particular pixel) to the coordinates of the ship track. Every fifth row of pixels in the image has one GCP

Examples of the GCP points:

(1904.00P ,23086.00L) -> (-67.1203008E,42.0984950N,0.00)
(0.00P ,23081.00L) -> (-67.1203325E,42.0993950N,0.00)
(1904.00P ,23076.00L) -> (-67.1202966E,42.0984951N,0.00)
(0.00P ,23066.00L) -> (-67.1202984E,42.0993966N,0.00)

The tifs themselves have no GCS and are completely rectangular. However, the track to which I am trying to georeference them is curved, occasionally turns 90 degrees, and in some cases will double back on itself.

What i've been trying to do is the following (adapted from this stack overflow thread) :

# open tif
ds = gdal.Open(sonar_tif, gdal.GA_Update)

# Get Source Projections (PCS)
ds_prj = ds.GetProjection() 
if ds_prj:
    print("Source Projection:",ds_prj,"\n")
    print("Source Projection: None\n")

# new instance of SpatialReference class from osr module
drs = osr.SpatialReference() 
print("drs is",drs)

# Set Geographic Coordinate Reference System (GCS)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326) # WGS84 GCS

if srs.IsProjected:
    print ("Target Projection:",srs.GetAttrValue('projcs'))

print ("Target GCS:",srs.GetAttrValue('geogcs'))

# Attach GCP points to tif
ds.SetGCPs(gcp_list, ds_prj) 

ds = None # closes dataset, unlocks files
d_out = None

If I understand correctly, this should perform the same action as gdal_translate, and I should have embedded GCP points and a geotiff. However, the output I get is a big black square.

Am I doing something fundamentally incorrect?

Is there a coordinate system mismatch that I might be unaware of?

Do I have to run gdal_warp before I can see a fully georeferenced image?

| improve this question | | | | |
  • This seems like a hard task or atleast a complicated transformation... I do think you need to run gdalwarp to produce a new raster. Otherwise you still have the same structure of the ribbon, which I am guessing is not what you want. You might need to use the -tps option of gdalwarp to achieve this complex transformation, but maybe a polynomial will work? A few questions: is d_out the output dataset? Where is it defined & when is data added? Why is drs created? If the rasters don't have geolocation information then srs should be provided to SetGCPs not ds_prj I suspect. – Logan Byers Sep 1 '17 at 21: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.