reconstruction.py

import laspy
import os
from owslib.wfs import WebFeatureService
import requests
import time
import zipfile
import io
from shapely.geometry import box
import shapely.wkt
import yaml
import subprocess
from osgeo import ogr
import json

f = open("settings.json", "r")
settings = json.load(f)

def download_bgt(b, bgt_out):
    postmsg = {
      "featuretypes": [
        "begroeidterreindeel",
        "kunstwerkdeel",
        "onbegroeidterreindeel",
        "ondersteunendwaterdeel",
        "ondersteunendwegdeel",
        "overbruggingsdeel",
        "overigbouwwerk",
        "pand",
        "scheiding",
        "tunneldeel",
        "waterdeel",
        "wegdeel"
      ],
      "format": "citygml",
      "geofilter": b.wkt
    }
    
    r = requests.post("https://api.pdok.nl/lv/bgt/download/v1_0/full/custom", json=postmsg)
    
    
    status = ""
    while(status != "COMPLETED"):
        r2 = requests.get("https://api.pdok.nl" + r.json()["_links"]["status"]["href"])
        status = r2.json()["status"]
        
        if status != "COMPLETED":
            print("not yet")
            time.sleep(5)
            
    print("done")
    
    bgtzip = requests.get("https://api.pdok.nl" + r2.json()["_links"]["download"]["href"]).content
    
    z = zipfile.ZipFile(io.BytesIO(bgtzip))
    
    z.extractall(bgt_out)

def find_ahn(bbox1):
    for root, dirs, files in os.walk("./ahn"):
        for file in files:
            if file[-3:].lower() == "laz":

                infile = laspy.file.File(root + "/" + file, mode="r-")
                bbox = infile.header.min[:2] + infile.header.max[:2]

                bbox2 = box(bbox[0], bbox[1], bbox[2], bbox[3])
                
                if bbox2.contains(bbox1):
                    print(root + "/" + file)
                    return(root + "/" + file)
                return None
                
    
def threedfy(pc_path):     
    with open("../config.yml") as f:
        config = yaml.load(f)
        
        pc_path = "data" + pc_path[1:]
        
        config["input_elevation"][0]["datasets"][0] = pc_path
        print(pc_path)
    
    with open("../config.yml", "w") as f:
        yaml.dump(config, f)
        
    fn = pc_path.split("/")[-1].split(".")[0]
    
    # go 1 directory up, after cmd go back again
    cwd = os.getcwd()
    os.chdir("..")
    cmd = "3dfier config.yml --CityJSON output/" + fn + ".json"
    
    # run 3dfier
    execute_cmd(cmd)
    os.chdir(cwd)
        

def execute_cmd(cmd):
    op = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    R = op.poll()
    print(R)
    if R:
        res = op.communicate()
        raise ValueError(res[1])
    re =  op.communicate()[0]
        

#wfs = WebFeatureService(url="https://geodata.nationaalgeoregister.nl/bag/wfs/v1_1?request=getCapabilities&service=WFS", version="1.1.0")

def create_terrain_polygon(bbox):
    bgt_geoms = []
    bgt_path = "./bgt/"
    for filename in os.listdir(bgt_path):
        if filename[-4:] == ".gml":
            bgt_geoms.append(ogr.Open(bgt_path + filename))
            
            gml = (open(bgt_path + filename, 'r')).read()
    
    
def create_cms(pc_in_file, pc_in_dir, downloadBGT, bgt_out):
    pc_path = "/".join([pc_in_dir, pc_in_file])
    
    pc = laspy.file.File(pc_path, mode="r-")
    bbox = pc.header.min[:2] + pc.header.max[:2]
    print(bbox)
    
    # create polygon from min/max
    bbox_poly = box(bbox[0], bbox[1], bbox[2], bbox[3])
    
    # find corresponding AHN file
    ahn_path = find_ahn(bbox_poly)
    
    """
    if downloadBAG == True:
        response = wfs.getfeature(typename="bag:pand", bbox=bbox_poly, outputFormat="geojson")
        out = open(file[:-4] + "_bag.geojson", "wb")
        out.write(response.read())
        out.close()
    """
        
    # only go on with the process if there as an AHN file to do the comparison with
    if ahn_path != None:
        if downloadBGT == True:
            # retrieve BGT from API and convert to GPKG
            download_bgt(bbox_poly, bgt_out)
        #p = Popen("./bgt/BGT_conversion.bat")
        #stdout, stderr = p.communicate()
        
        #bgt_path = "./bgt/"
        
        """
        for root, dirs, files in os.walk(bgt_path):
            for file in files:
                if file[-4:].lower() == "gpkg":
                    name = file.split(".")[0]
                    fout = name + "_subset.gpkg"
                    cmd = []
                    cmd.append("ogr2ogr -f gpkg -spat")
                    cmd.append(" ".join(map(str, bbox)))
                    cmd.append(bgt_path + fout)
                    cmd.append(bgt_path + file)
                    cmd = " ".join(cmd)
                    execute_cmd(cmd)
                    print(cmd)
        """

        threedfy(pc_path)
        print("First 3D reconstruction done")
        threedfy(ahn_path)
        print("Second 3D reconstruction done")

    # create polygon that can be used by 3dfier to create terrain
    #create_terrain_polygon(bbox)