4

I have some scan files with spatial information.
My Original Image Original Image
In the example image, the rectangular area drawn with the blue line is my area containing spatial information. In the four corners of the rectangle, there is information that gives the coordinate value of the corner it is in.
Top Left Of Raster
Top Left
Top Right Of Raster
Top Right
Bottom Left Of Raster
Bottom Left
Bottom Right Of Raster
Bottom Right
As you can see, the coordinates of the EPSG:2320 Coordinate System are written on all 4 corners. I tried to make georeference at Python code below. But it gave me a very different result. In the program it is going to the real-world domain of the raster. However, it rotates and inverts the image even though I enter the position of the corners in the image and the coordinate information in the real world correctly.

My code:

    """ PINDEX CORNER X,Y(IMAGE) AND EPSG:2320 COORDİNATES 

                    IN IMAGE  | EPSG:2320
                    ----------|----------------------
    TOP LEFT:       796, 569  | 443893.7, 4420713.72
    TOP RIGHT       5005, 598 | 444428.05,4420709.81
    BOTTOM RIGHT    4965,6065 | 444423,4420015.85
    BOTTOM LEFT     754, 6035 | 443888.6,4420019.76
    """

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

    file_path="ciktilar/"+i

    tl = GroundControlPoint(796, 569, 443893.7, 4420713.72)   #top left
    bl = GroundControlPoint(5005, 598, 444428.05,4420709.81)    #top right
    br = GroundControlPoint(4965,6065, 444423,4420015.85)     #bottom right
    tr = GroundControlPoint(754, 6035, 443888.6,4420019.76)    #bottom left
    gcps = [tl, bl, br, tr]

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

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


My Result GeoTIF
My Result GeoTIF
The result image was added to the GIS program and taken from the screen.
How do I do the correct Georeference?
Any ideas?

3
  • Please use the native StackExchange image system for images; folks aren't likely to follow links to external sites, and automated repair can be performed when something changes in the SE stack
    – Vince
    Commented Sep 6, 2022 at 19:58
  • Thank you. I rearranged the images Commented Sep 6, 2022 at 20:09
  • Doesn't it seem that rasterio takes the pixel coordinates in another order than GDAL? Compare rasterio.readthedocs.io/en/latest/topics/… row - col and gdal.org/programs/… pixel - line.
    – user30184
    Commented Sep 6, 2022 at 20:58

2 Answers 2

4

You have two alternatives:

  1. Change the order of pixel coordinates into row-column order for rasterio as documented in https://rasterio.readthedocs.io/en/latest/topics/transforms.html#using-ground-control-points
    tl = GroundControlPoint(569, 796, 443893.7, 4420713.72)   #top left
        bl = GroundControlPoint(598, 5005, 444428.05,4420709.81)    #top right
        br = GroundControlPoint(6065, 4965, 444423,4420015.85)     #bottom right
        tr = GroundControlPoint(6035, 754, 443888.6,4420019.76)    #bottom left
  1. Keep the order of pixel coordinates as they are but add ground control points with GDAL, for example https://gdal.org/programs/gdal_translate.html#cmdoption-gdal_translate-gcp

gdal_translate -gcp 796 569 443893.7 4420713.72 -gcp 5005 598 444428.05 420709.81 -gcp 4965 6065 444423 4420015.85 -gpc 754 6035 443888.6 420019.76 input.tif with_gcp.tif

Then you can use gdalwarp

gdalwarp -t_srs epsg:2320 with_gcp.tif final.tif

1
  • Thank you for the information you provided. As a result of the answer I improved my code and added the solution as the answer Commented Sep 6, 2022 at 23:47
3

My Solution
user30184 ' Answer
thank you for your ideas. Your answer helped me a lot to find the solution. As you said, I added the ground control points to a list with the gdal.GCP function and used that list in the transformation.
Here is the final version of my code

""" PINDEX CORNER X,Y(IMAGE) AND EPSG:2320 COORDİNATES 

                IN IMAGE  | EPSG:2320
                ----------|----------------------
TOP LEFT:       796, 569  | 443893.7, 4420713.72
TOP RIGHT       5005, 598 | 444428.05,4420709.81
BOTTOM RIGHT    4965,6065 | 444423,4420015.85
BOTTOM LEFT     754, 6035 | 443888.6,4420019.76
"""
from osgeo import gdal 
import osr
from rasterio.control import GroundControlPoint
dataset = gdal.Open(file_path) 
I = dataset.ReadAsArray(0,0,dataset.RasterXSize,dataset.RasterYSize)

outdataset = gdal.GetDriverByName('GTiff') 
output_SRS = osr.SpatialReference() 
output_SRS.SetProjection(2320)
outdataset = outdataset.Create("ciktilar/donusmus/"+i,dataset.RasterXSize,dataset.RasterYSize,I.shape[0]) 
for nb_band in range(I.shape[0]):
    outdataset.GetRasterBand(nb_band+1).WriteArray(I[nb_band,:,:])
    
    #My Ground Control Points
topleft = GroundControlPoint(796, 569, 443893.7, 4420713.72)   #top left
topright = GroundControlPoint(5005, 598, 444428.05,4420709.81)    #top right
bottomright = GroundControlPoint(4965,6065, 444423,4420015.85)     #bottom right
bottomleft = GroundControlPoint(754, 6035, 443888.6,4420019.76)    #bottom left

gcp_list = [] 
    #Repeat this step for as many GCP as you want (there is a limit but I don't know the number of GCPs)

gcp_list.append(gdal.GCP(443893.7, 4420713.72,0,569,796))   #sorted x , y , z , col , row
gcp_list.append(gdal.GCP(444423,4420015.85, 0,4965,6065))#sorted x , y , z , col , row
gcp_list.append(gdal.GCP(444428.05,4420709.81, 0,5005,598 ))#sorted x , y , z , col , row
gcp_list.append(gdal.GCP(443888.6,4420019.76 ,0,754,6035 ))#sorted x , y , z , col , row

outdataset.SetProjection(output_SRS.ExportToWkt()) 
wkt = outdataset.GetProjection() 
outdataset.SetGCPs(gcp_list,wkt)
outdataset = None

My Result GeoTIF with Google Satellite on GIS Program enter image description here Don't mind the bad image. I reduced the size of the photo as it has a maximum 2mb limit. The result is excellent.

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.