Add z-value to shapefile made of multipatch data (roofsurfaces)

Jacob_JeffB · ‎01-27-2018

Hey guys, i am pretty new in 3D analysis in ArcMap/Scene (10.5.1.)

It might be pretty simple question ( and it hopefully is).

I am working with 3D buildings of my hometown to analyze roof solar potential.

My final goal is to produce slope and aspect data.

What kind of data i use:

- 109 multipatch (.gml) files (each contains GroundSurface, RoofSurface, Wallsurface)

>> i am interested in Roofsurfaces to analyze slope and aspect

- 1 shapefile (it contains the 109 roofsurfaces)

>> i produced it myself and used it for 2D analysis (for example: footprint to estimate and distinguish flat and sloping roof types)

What i am trying to do is:

- add z_value to the shapefile

- to create a raster elevation map

- and finally do slope and aspect analysis

My problem is:

- i dont know much about multipatch and gml. data

- i dont know how to create this z-value out of my existing data and add it the attribute table of the shapefile (if had the z-value i could use polygon to raster)

> Is it even possible? Can i use the shapefile? Or do i have to use the raw data (multipatch / gml.)

I hope you can understand what i am trying to do, and what kind of data i use.

Thanks in advance,

JacobB

XanderBakker · ‎01-31-2018

Find below the code that I wrote to process the data (DWG files):

import arcpy

def main():
    import os

    # settings
    ws = r'C:\GeoNet\CityGML\DWG_LoD2'
    fgdb = r'C:\GeoNet\CityGML\gdb\Hannover.gdb'
    dsm_name = 'DSM'
    slope_name = 'Slope'
    aspect_name = 'Aspect'
    cellsize = 1.0
    sr = arcpy.SpatialReference(25832)  # ETRS_1989_UTM_Zone_32N

    # necessary environment settngs
    arcpy.env.outputCoordinateSystem= sr
    arcpy.env.workspace = ws

    # create list of DWG file in folder (ws)
    lst = arcpy.ListFiles("*.DWG")

    # check out SA license
    arcpy.CheckOutExtension("Spatial")

    # loop through dwg file
    lst_ras = []
    lst_err = []
    for dwg in lst:
        dwg_mp = os.path.join(ws, dwg, 'MultiPatch')
        print("Processing: {}".format(dwg))

        # correct extent
        ext = arcpy.Describe(dwg_mp).extent
        ext_corr = CorrectExtent(ext, cellsize)
        arcpy.env.extent = ext_corr

        # define output raster name
        ras_name = 'ras{}'.format(os.path.splitext(dwg)[0])
        out_ras = os.path.join(fgdb, ras_name)

        # create raster from multipatch and add to lost
        try:
            # some files create a error. I will skip those
            arcpy.MultipatchToRaster_conversion(dwg_mp, out_ras, cellsize)
            lst_ras.append(out_ras)
        except Exception as e:
            print "ERROR: ", e
            lst_err.append(dwg)


    print "DWG that produced Errors:"
    print "\n - ".join(lst_err)

    # mosaic raster into new raster
    print("mosaic into new raster...")
    input_ras = ';'.join(lst_ras)
    arcpy.MosaicToNewRaster_management(input_ras, fgdb, dsm_name, sr, "32_BIT_FLOAT", cellsize, 1, "MEAN")

    # set extent to output DSM
    dsm_ras = os.path.join(fgdb, dsm_name)
    arcpy.env.extent = dsm_ras

    # calculate slope
    print("calculate slope...")
    slope_ras = os.path.join(fgdb, slope_name)
    out_slope = arcpy.sa.Slope(dsm_ras, "PERCENT_RISE")  # or use DEGREE
    out_slope.save(slope_ras)

    # calculate aspect
    print("calculate aspect...")
    aspect_ras = os.path.join(fgdb, aspect_name)
    out_aspect = arcpy.sa.Aspect(dsm_ras)
    out_aspect.save(aspect_ras)


def CorrectExtent(ext, cellsize):
    ext.XMin = divmod(ext.XMin, cellsize)[0]
    ext.YMin = divmod(ext.YMin, cellsize)[0]
    ext.XMax = divmod(ext.XMax, cellsize)[0] + cellsize
    ext.YMax = divmod(ext.YMax, cellsize)[0] + cellsize
    return ext

if __name__ == '__main__':
    main()
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

The file geodatabase with the resulting rasters (DSM, slope, aspect), is about 1.5GB (compressed a little over 210MB). You can download it here: HannoverCityGML.zip - Google Drive

View solution in original post

DanPatterson_Retired · ‎01-27-2018

You don't happen to have the Data Interoperability extension or FME do you?

'gml' keeps hitting those, for example http://desktop.arcgis.com/en/arcmap/latest/extensions/data-interoperability/exercise-2a-importing-da...

I will move this to the https://community.esri.com/community/gis/3d?sr=search&searchId=80bf69cf-9493-4690-8280-71d3ae3693b5&...‌ space since Geonet Help is about getting help on how to use GeoNet

XanderBakker · ‎01-27-2018

If you don't have DataInterop you could start here: How To: Access GML data without a Data Interoperability license

Jacob_JeffB · ‎01-27-2018

Hey Dan,

thanks for your response. I do have the Data Interoperability extension. Ive alreadyimported and converted the gml. to shapefile to do my 2d analysis. I am now interested in the 3D data (z.values) to create a raster terrain map and do the slope and aspect analysis.

XanderBakker · ‎01-27-2018

Can you share the data to have al look and determine the best way to generate the DSM?

Jacob_JeffB · ‎01-30-2018

Hey Xander. I am pretty new to GeoNet (whats the best way to share the data in a question?).

As far as i figured out, i have to use the "Multipatch to Raster" - Tool to convert the multipatch (gml.) data to raster.

It didnt work yet. (ERROR 010092: Invalid output extent.)

XanderBakker · ‎01-30-2018

Hi Jacob Bernhardt ,

When you reply to a question, you can find a link in the upper right corner of where you are editing that states "Use advanced editor". When you click on the link the advanced editor will open and in the lower left section you will find the link "Attach" to attach the file to the thread.

Just ZIP the files that you want to share into a single ZIP and attach it to the thread.

Jacob_JeffB · ‎01-30-2018

As you can see in the picture the data format is .gml (catalog). Every dataset contains different building types. I am interested in the "Roof Surface Multipatch" features. You can also the the Error in the result window, that pop up after i used the "Multipatch to Raster"-Tool. I already tried to change the output folder / gdb / server and harddrive. Any ideas?

This is the pythoncode: (the code isself seems to work)

# Name: MultipatchToRaster_Ex_02.py
# Description: Converts multipatch features to a raster dataset.

# Import system modules
import arcpy
from arcpy import env

# Set environment settings
env.workspace = "I:\NeueEnergie2050\Stadtmodell_Hannover_CityGML_LoD2\GML_CityGML_LoD2_12_12_17\5410_5806.gml"

# Set local variables
inFeatures = "RoofSurface Multipatch"
outRaster = "I:\NeueEnergie2050\PV_Analyse\GML_Stadtgebiet Hannover\RoofSurface_raster/5410_5806.tif"
cellSize = 1.0

# Execute MultipatchToRaster
arcpy.MultipatchToRaster_conversion(inFeatures, outRaster, cellSize)

XanderBakker · ‎01-30-2018

The error message that you shared refers to the extent not being valid. That should be the first part where you should look and also the coordinate system in which the data is stored, which I believe is ETRS89 (UTM), EPSG-Code 25832 if the data on this page is correct: Digitales 3D-Stadtmodell | 3D-Stadtmodell und Geländemodell | Open GeoData | Geoinformation | Fachbe...

Since the data is available for download, I will download it do some tests to see if it allows me to create a DSM.

XanderBakker · ‎01-30-2018

I did a quick test and the most important aspect is to define the processing extent and coordinate system in the geoprocessing environment and then you can run the tool:

Extract of result:

A script that processes all the gml files and creates a single DSM will require some more programming.