Convert PNG shape into KML or GeoJson

user_1012234 0
Viewed 2326

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

geojson kml opencv image-processing python
1 Answers

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
user_28366210
Related