Search code examples
pythonopencvimage-processinggeojsonkml

Convert PNG shape into KML or GeoJson


I have thousands of shapes stored as PNG files and boundaries' coordinates for each shape. Boundaries' coordinates are coordinates of 4 corners of the minimum enclosing rectangle of the shape (example below).

The goal is to use PNG images and their boundaries' coordinates to convert them into polygon (KML or GeoJSON).

I'm not sure even about the techs I can use to reach the result, so I'd appreciate any suggestions.

Input data (PNG):

Alt text

  • Coordinates of 4 corners of the minimum enclosing rectangle of the shape: 8.348236, 44.66804, 8.305321, 44.66829, 8.348579, 44.63507, 8.305492, 44.63507.

Desired output:

  • Polygon is a Gist that shows the result of interpreting the filled area of the PNG located in the right place on the map. Click on the Display the source blob to see the raw GeoJSON.

How do I imagine the process:

  • Step 1: we have a PNG image and 4 points. This let us place the PNG image on the map in the right place and scale it appropriately.
  • Step 2: we recognize the locations of shape's key points.
  • Step 3: we extract a set of recognized points into the polygon.

Alt text

I used simple PNG as an example but the shapes could be much more complex:

Alt text


Solution

  • Ok, I saved your image as "shape.png" and your GeoJSON enclosing rectangle as "boundaries.json". Then my method is as follows:

    • get the North, East, South and West limits in terms of latitude and longitude
    • load and trim the shape image to get rid of all black borders, threshold to pure black and white
    • work out the X and Y scaling from pixels to degrees by looking at the image width and height in pixels and degrees
    • use OpenCV findContours() to find the vertices in the shape image
    • translate all the vertices I find from image coordinates to latitude, longitude
    • write those points out to a JSON results file.

    #!/usr/bin/env python3
    
    import cv2
    import json
    import geojson
    import numpy as np
    from geojson import Feature, Point, FeatureCollection, Polygon, dump
    
    def getNESWextents(GeoJSONfile):
    
        # Load the enclosing rectangle JSON
        with open('boundaries.json','r') as datafile:
            data = json.load(datafile)
        feature_collection = FeatureCollection(data['features'])
    
        lats = []
        lons = []
        for feature in data['features']:
            coords = feature['geometry']['coordinates']
            lons.append(coords[0])
            lats.append(coords[1])
    
        # Work out N, E, S, W extents of boundaries
        Nextent = max(lats)
        Sextent = min(lats)
        Wextent = min(lons)
        Eextent = max(lons)
        return Nextent, Eextent, Sextent, Wextent
    
    def loadAndTrimImage(imagefilename):
        """Loads the named image and trims it to the extent of its content"""
        # Open shape image and extract alpha channel
        im = cv2.imread(imagefilename,cv2.IMREAD_UNCHANGED)
        alpha = im[...,3]
        # Find where non-zero, i.e. not black
        y_nonzero, x_nonzero = np.nonzero(alpha)
        # Crop to extent of non-black pixels and return
        res = alpha[np.min(y_nonzero):np.max(y_nonzero), np.min(x_nonzero):np.max(x_nonzero)]
    
        # Threshold to pure white on black
        _, res = cv2.threshold(res, 64, 255, cv2.THRESH_BINARY)
        return res
    
    def getVertices(im):
        """Gets the vertices of the shape in im"""
    
        _, contours, *_ = cv2.findContours(im, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    
        # Should probably sort by contour area here - and take contour with largest area
        perim = cv2.arcLength(contours[0], True)
        approx = cv2.approxPolyDP(contours[0], 0.01 * perim, True)
    
        print(f'DEBUG: Found shape with {approx.shape[0]} vertices')
        return approx
    
    if __name__ == "__main__":
    
        # Get N, E, S, W extents from JSON file
        Nextent, Eextent, Sextent, Wextent = getNESWextents('boundaries.json')
        print(f'DEBUG: Nextent={Nextent}, Eextent={Eextent}, Sextent={Sextent}, Wextent={Wextent}')
    
        # Load the image and crop to contents
        im = loadAndTrimImage('shape.png')
        print('DEBUG: Trimmed image is "trimmed.png"')
        cv2.imwrite('trimmed.png', im)
    
        # Get width and height in pixels
        Hpx, Wpx = im.shape
        # Get width and height in degrees
        Hdeg, Wdeg = Nextent-Sextent, Eextent-Wextent
        # Calculate degrees per pixel in East-West and North-South direction
        degppEW = Wdeg/Wpx
        degppNS = Hdeg/Hpx
        print(f'DEBUG: degppEW={degppEW}, degppNS={degppNS}')
    
        # Get vertices of shape and stuff into list of features
        features = []
        vertices = getVertices(im)
        for i in range(vertices.shape[0]):
           x, y = vertices[i,0]
           lon = Wextent + x*degppEW
           lat = Nextent - y*degppNS
           print(f'DEBUG: Vertex {i}: imageX={x}, imageY={y}, lon={lon}, lat={lat}')
           point = Point((lon,lat))
           features.append(Feature(geometry=point, properties={"key":"value"}))
    
        # Convert list of features into a FeatureCollection and write to disk
        featureCol = FeatureCollection(features)
        with open ('result.json', 'w') as f:
            dump(featureCol, f)
    

    Here is the trimmed image:

    enter image description here

    Here is the debug output:

    DEBUG: Nextent=44.66828662253787, Eextent=8.348579406738281, Sextent=44.63507036301143, Wextent=8.305320739746094
    DEBUG: Trimmed image is "trimmed.png"
    DEBUG: degppEW=8.634464469498503e-05, degppNS=6.0503204966194347e-05
    DEBUG: Found shape with 6 vertices
    DEBUG: Vertex 0: imageX=211, imageY=2, lon=8.323539459776736, lat=44.668165616127936
    DEBUG: Vertex 1: imageX=2, imageY=224, lon=8.305493429035483, lat=44.654733904625445
    DEBUG: Vertex 2: imageX=81, imageY=472, lon=8.312314655966388, lat=44.63972910979383
    DEBUG: Vertex 3: imageX=374, imageY=548, lon=8.337613636862018, lat=44.63513086621639
    DEBUG: Vertex 4: imageX=500, imageY=392, lon=8.348493062093587, lat=44.64456936619112
    DEBUG: Vertex 5: imageX=484, imageY=155, lon=8.347111547778466, lat=44.65890862576811