1

I have a binary mask with buildings' footprints. What I have is a latitude and longitude of the center of my image and I want to georeference each polygon first. And here is what I have done before:

  1. I get the contours (pixel coordinates) for each polygon and store them in a data frame using OpenCV contours = measure.find_contours(self.bin_mask, 0.5).

  2. Now I need to get the georeferenced coordinate. I have already done that before with rasterio: (px, py) = rio.transform.xy(self.tif_img.transform, x, y, offset='center')

Where I get the georeferenced coordinates for each pixel in a polygon Everything was working great and I was able to extract the building, the issue now is that I don't have the tiff image.

Is it possible to get the geo affine (transform) without the referenced TIFF?

In short, I need a way to georeferenced the image pixels and all I have is:

1- a lat/long coordinates of the center of the image.

2- The images size 650X650 and 1 pixel = 0.297 m

3- coordinates is (-95.407189656922 29.893852182227928)

Update

1- I know that I need to get the bounding box for my image. From what I understood so far, this will help me get the geo transform

2- I also tried this solution but it didn't work, I think I need more than one GCP to implement this correctly.

5
  • I found out a solution, however, I have to edit my answer because I had to modify the code for considering all detected contours. Your updated data (center of the image and resolution) was enough for getting a desired result.
    – xunilk
    Feb 11 at 21:55
  • Thanks! I tried georeferencing my image using GDAL in qgis, but I need more than one GCP point to complete the process. So I'm looking for other methods at the moment
    – salRad
    Feb 12 at 4:43
  • It is not necessary because your provided data is enough. On the other hand, I edited completely my answer with a modified code and it works. I hope it helps.
    – xunilk
    Feb 13 at 3:05
  • Thanks a lot! I have to say I learned a lot from this :)
    – salRad
    Feb 13 at 7:25
  • 1
    Glad to help you.
    – xunilk
    Feb 13 at 13:24
3

OpenCV is a powerful and useful tool for 'contour detection' and convert them in vector layers. However, they are produced in image coordinates (i, j indices); not in map coordinates. For overcoming this issue you need to have your binary mask with buildings' footprints georeferenced by using your specific reference building.

With your data (center of the image and resolution) it is possible to find out which is your interest area. Assuming a projection in meters with EPSG:32615, the center of the image is (267550.7822539393 3309457.9713751674). So, xmin = 267550.7822539393 - 0.297*325 = 267454.257253939 and ymax = 3309457.9713751674 + 0.297*325 = 3309554.496375167. Now, we have all necessary geotransform parameters: xmin, ymax, xsize, and ysize

For my approach and additional corroboration, I determined all points of bounding box of complete area and digitized over them its respective polygon. By using rasterize (vector to raster) tool of QGIS with this polygon, I got a raster with 650x650 size as expected.

xmin, ymax point(267454.257253939 3309554.496375167)
xmin, ymin point(267454.257253939 3309361.446375167)
xmax, ymax point(267647.307253939 3309554.496375167)
xmax, ymin point(267647.307253939 3309361.446375167)

By using an image captured from Google Satellite, I clipped it with referred polygon as follows.

enter image description here

Afterward, I digitized arbitrarily three building in that area; as it can be observed in following image.

enter image description here

Deselecting raster layer; as in following image:

enter image description here

I used gdal_translate command for producing my georeferenced mask (prov.tif) with three reference buildings.

import os

layer = iface.activeLayer()
epsg = layer.crs().postgisSrid()
        
mapcanvas = iface.mapCanvas()
output_raster = 'prov.tif'
input_raster = output_raster[:-4] + ".png"
mapcanvas.saveAsImage(input_raster)
        
extent = mapcanvas.extent()
xmin, ymin, xmax, ymax = extent.toRectF().getCoords()

cmd = "gdal_translate -a_srs EPSG:" + str(epsg) + " -a_ullr " + str(xmin) + " " \
                                                              + str(ymax) + " " \
                                                              + str(xmax) + " " \
                                                              + str(ymin) + " " \
                                                              + input_raster + " " \
                                                              + output_raster

os.system(cmd)
os.remove(input_raster)

I directly used my reference building mask (named prov.tif) for detecting its contours with following code (adapted from here) and to produce its respective polygon vector layer. Observe that this image is used for obtaining geotransform parameters with gdal python module and converting image coordinates in map coordinates.

from osgeo import gdal
import cv2

raster_file = "prov.tif"

data = gdal.Open(raster_file)

geo_transform = data.GetGeoTransform()

image=cv2.imread(raster_file)
#cv2.imshow('input image',image)

orig_image = image.copy()

gray=cv2.cvtColor(image,cv2.COLOR_BGR2GRAY)
ret, thresh=cv2.threshold(gray,127,255,cv2.THRESH_BINARY_INV)

_, contours, hierarchy=cv2.findContours(thresh.copy(),cv2.RETR_LIST,cv2.CHAIN_APPROX_NONE)

areas = []

for c in contours:
    area = cv2.contourArea(c)
    areas.append(area)

max_area = max(areas)

xmin = geo_transform[0]
ymax = geo_transform[3]
xsize = geo_transform[1]
ysize = geo_transform[5]

polygons = []

for i, area in enumerate(areas):
    for c in contours:
        area = cv2.contourArea(c)
        x,y,w,h=cv2.boundingRect(c)
        cv2.rectangle(orig_image,(x,y),(x+w,y+h),(0,0,255),2)
        #calculate accuracy as a percent of contour perimeter
        accuracy=0.01*cv2.arcLength(c,True)
        approx=cv2.approxPolyDP(c,accuracy,True)
        if area == areas[i]:
            pol = c
            cv2.drawContours(image,[approx],0,(0,255,0),2)

#cv2.imshow('Approx polyDP', image)

        points = []
        points2 = []

        for point in pol:
            points2.append(QgsPointXY(xmin + xsize*point[0][0], ymax + ysize*point[0][1]))

        polygon = [points2]
        
        polygons.append(polygon)

epsg = 32615

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'polygon_32615',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(polygons)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPolygonXY(polygons[i]))

prov.addFeatures(feats)

QgsProject.instance().addMapLayer(mem_layer)

registry = QgsProject.instance()

layer = registry.mapLayersByName('polygon_32615')

print("Done!")

After running above code, it was produced a vector layer with adequate contour in expected positions. It has three times detected features because prov.tif is a RGB image.

enter image description here

8
  • Thanks a lot for the great answer! I'll try this out
    – salRad
    Feb 11 at 4:07
  • One question, you started by saying that you used a georeferenced image. But in my case, I don't have any so I'm not sure how to utilize this.
    – salRad
    Feb 11 at 4:10
  • 1
    You need to produce them by using gdal_translate command. I will to edit my answer to put an example.
    – xunilk
    Feb 11 at 4:17
  • 2
    It will be tomorrow. It's late in my country.
    – xunilk
    Feb 11 at 4:23
  • 1
    Can I have latitude and longitude of the center of your image and the spatial resolution?
    – xunilk
    Feb 11 at 12:47

Your Answer

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

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