13

I have 6-band unregistered/unreferenced TIFF images that I have been trying to export as georeferenced GeoTIFFs. The GPS coordinates of the corners of each image is known and each image is 1024 height and 1280 width. I attempted the GDAL solution in https://stackoverflow.com/questions/55681995/how-to-georeference-an-unreferenced-aerial-image-using-ground-control-points-in but was receiving item has no get_attr_ errors at the setGCP command.

Attempt 1:

import rasterio as rio
from rasterio import Affine
from rasterio.warp import reproject, Resampling, calculate_default_transform
from rasterio.transform import from_origin
from rasterio.control import GroundControlPoint

im = rio.open('path to a tif')
arr = im.read()
dst_fn1 = outdir + 'out.tif'

tl = GroundControlPoint(0, 0, -83.70113410256013, 42.307951446432604)
bl = GroundControlPoint(1024, 0, -83.69940501521428, 42.307603183805234)
br = GroundControlPoint(1024, 1280, -83.698829074736, 42.3091785425499)
tr = GroundControlPoint(0, 1280, -83.70055820297041, 42.309526812647555)
gcps = [tl, bl, br, tr]

src_crs = im._get_crs()  ## This is a blank object. No crs 
dst_crs = {'init': 'epsg:4326'}
src_tran = im.transform
dst_transform, dst_width, dst_height = calculate_default_transform(a, crs, 1280, 1024, gcps=gcps)
dest_ds = np.empty(shape=(arr.shape[0], dst_width, dst_height))
reproject(sand, dest_ds, src_crs=dst_crs, dst_crs=dst_crs, 
      src_transform=src_tran, dst_transform=dst_transform, 
      resampling=Resampling.nearest)

This method does not throw any errors but the resulting array is only zeroes. Also exported array is not rotated to match where the corners should be. Perhaps this is the wrong expectation but I thought that the data would be inside the yellow box and the outside would be nodata values to square it off.

enter image description here

Attempt 2:

( This is essentially the method found here https://www.earthdatascience.org/courses/earth-analytics-python/lidar-raster-data/reproject-raster/)

with rio.open('path to a tif') as src:
    transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, gcps=gcps)
    kwargs = src.meta.copy()
    kwargs.update({'crs': dst_crs, 'transform': transform, 'width': width, 'height': height})
        with rio.open('path/output.tif', 'w', **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rio.band(src, i),
                    destination=rio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform, dst_crs=dst_crs,
                    resampling=Resampling.nearest)

Results in this error:

rasterio._err.CPLE_AppDefinedError: The transformation is already "north up" or a transformation between pixel/line and georeferenced coordinates cannot be computed for TIFFPATH. There is no affine transformation and no GCPs. Specify transformation option SRC_METHOD=NO_GEOTRANSFORM to bypass this check.
1
  • 5
    This is ill-posed: you should not reproject a file with no projection, but simply set a projection. Using only gdal you would need to specify the coordinates of the bottom left corner and the pixel size (in crs units) in the opened_band.geoTransform and the crs in opened_band.projection. See gis.stackexchange.com/a/315301/98816 Oct 6, 2020 at 15:07

1 Answer 1

6

Use rasterio.transform.from_gcps, then open the raster in append r+ mode and write the CRS and transform.

import rasterio as rio
from rasterio.transform import from_gcps
from rasterio.control import GroundControlPoint


tl = GroundControlPoint(0, 0, -83.70113410256013, 42.307951446432604)
bl = GroundControlPoint(1024, 0, -83.69940501521428, 42.307603183805234)
br = GroundControlPoint(1024, 1280, -83.698829074736, 42.3091785425499)
tr = GroundControlPoint(0, 1280, -83.70055820297041, 42.309526812647555)
gcps = [tl, bl, br, tr]

transform = from_gcps(gcps)
crs = 'epsg:4326'

with rio.open(filepath, 'r+') as ds:
    ds.crs = crs
    ds.transform = transform

enter image description here

1
  • 1
    Be aware that rasterio.transform.from_gcps tries to fit an affine linear transformation based on the GCP information. The precise transformation however, will usually not be linear. Instead of using the transform attribute of the dataset, you can store your list of GCPs in the gcps attribute to save the actual GCP information to the GeoTiff. Opening such a file with QGIS for example will then perform a more accurate reprojection based on thin plate splines.
    – Christoph
    Dec 1, 2023 at 9:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.