Get all GWR addresses on a given parcel using api3.geo.admin.ch services

get_gwr_from_parcel.py

#!/usr/bin/env python

# Get all GWR (ch.bfs.gebaeude_wohnungs_register) on a given parcel (ch.kantone.cadastralwebmap-farbe)


import sys
import re
import requests
import json
import logging
try:
    import esrijson
except ImportError:
    pass


logging.basicConfig(format='%(levelname)s:%(message)s', level=logging.INFO)


url = 'https://api3.geo.admin.ch'

parcel_name = 'CH573575634627'
# parcel_name = 'CH408463626576'


egid_name= '1230414,'

regex = r"BOX\((.*)\)"

# Given an EGID number, find the approximate location (bbox and center) using the `SearchService`
def get_parcel_location(parcel_name):
    logging.info("######## Search for Parcel <{}> location ##############".format(parcel_name))
    params = {
            'sr': 2056,
            'searchText': 'parcel ' + parcel_name,
            'type': 'locations'
            }
    
    
    resp = requests.get( url + '/rest/services/ech/SearchServer', params= params )
    
    req = requests.Request('GET',url + '/rest/services/ech/SearchServer', params= params)
    prepared = req.prepare()
    logging.debug("url={}".format(prepared.url))
    
    res = resp.json()['results'][0]
    
    x,y = res['attrs']['x'], res['attrs']['y']
    box = res['attrs']['geom_st_box2d']
    
    matches = re.findall(regex, box, re.MULTILINE)
    
    extent = matches[0].replace(" ", ",")
    
    logging.info("parcel: {},{}".format(y,x))
    logging.info("bbox: {}".format(extent))

    return ((y,x), extent)


# Get the parcel geometry using the `Identify` service
def get_geometry(point, extent, layer='ch.kantone.cadastralwebmap-farbe'):
    logging.debug("=========== Get geometry of Parcel ================")
    y,x = point
    params = {
            'geometry': '{},{}'.format(y,x),
            'geometryFormat': 'geojson',
            'geometryType': 'esriGeometryPoint',
            'imageDisplay': '833,949,96',
            'layers': 'all:' + layer,
            'mapExtent': extent,
            'returnGeometry': True,
            'sr': 2056, 
            'tolerance':10
    
            }
    
    req = requests.Request('GET',url + '/rest/services/all/MapServer/identify', params=params)
    prepared = req.prepare()
    logging.debug("url={}".format(prepared.url))
    resp = requests.get(url + '/rest/services/all/MapServer/identify', params=params)
    
    res = resp.json()['results'][0]
    
    logging.debug(res)

    return res['geometry']


# Get all GWR inside a given geometry

def get_gwr(geom, extent):
    logging.debug("############# Get GWR on geometry ####################")
    
    
    rings = geom['coordinates']
    
    # Request are mad with EsriGeometriey
    esri_geom_ = { 'rings': rings }
    esri_geom_["spatialReference"] = {"wkid" : 2056}
    esri_geom = json.dumps(esri_geom_)
    logging.debug(esri_geom)
    
    # Just to check esri geom is valid
    if 'esrijson' in sys.modules:
        es = esrijson.loads(esri_geom)
        logging.info("Translating to EsriJson: {}".format(es))
        logging.info("Translating again to Shape: {}".format(esrijson.to_shape(es)))
    
    params = {
            'geometry': esri_geom,
            'geometryFormat': 'geojson',
            'geometryType': 'esriGeometryPolygon',
            'imageDisplay': '833,949,96',
            'layers': 'all:ch.bfs.gebaeude_wohnungs_register',
            'mapExtent': extent,
            'returnGeometry': True,
            'sr': 2056, 
            'tolerance':10
            }
    
    resp = requests.get(url + '/rest/services/api/MapServer/identify', params=params)
    
    req = requests.Request('GET',url + '/rest/services/api/MapServer/identify', params=params)
    prepared = req.prepare()
    logging.debug("url={}".format(prepared.url))



    return resp.json()['results']

def usage():
    print("Get EGID on a given parcel")
    print("{} <parcel nr>".format(__file__))
    print("   for instance {}".format("CH573575634627"))


def main():
    print(sys.argv)
    if len(sys.argv) != 2:
        usage()
        sys.exit(2)
    parcel_name = sys.argv[1]
    point, extent = get_parcel_location(parcel_name)

    geom = get_geometry(point, extent,)
    results = get_gwr(geom, extent)
    
    logging.info("===EGID for Parcel {}====".format(parcel_name))
    for res in results:
        #print(res)
        attrs = res['properties']
        logging.info("egid={}, strname={} {}".format(attrs['egid'], attrs['strname1'], attrs['deinr']))
    

if __name__ == '__main__':
    main()