Computing shadow of building footprint

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Antworten
100ri
User
Beiträge: 4
Registriert: Donnerstag 8. September 2016, 18:29

Dear All
I'm an international PhD student in the field of Climate change at Humboldt University. For my research I#m using ArcGIS 10.3. I need a tool which By luck I found the source as a Phyton code on the web:
http://gis.stackexchange.com/questions/ ... -arcgis-de
I am interested in importing this in ArcGIS as a toolbox. Its a code for computing the shadow of building footprint, considering the height of the building, at a certain date and time of the year. Unfortunately and sadly I have knowledge on scripting. Therefore I was wondering if someone could help me with that, or at least to check whether this code is correct? I appreciate your help in advance.

Cheers

CoDe:

Code: Alles auswählen

import arcpy
import os
import math

def message(msg, severity=0):
    # Adds a Message (in case this is run as a tool)
    # and also prints the message to the screen (standard output)
    #
    print msg

    # Split the message on \n first, so that if it's multiple lines,
    #  a GPMessage will be added for each line
    try:
        for string in msg.split('\n'):
            # Add appropriate geoprocessing message
            #
            if severity == 0:
                arcpy.AddMessage(string)
            elif severity == 1:
                arcpy.AddWarning(string)
            elif severity == 2:
                arcpy.AddError(string)
    except:
        pass

def main():
    arcpy.env.overwriteOutput=True

    # Read in parameters
    inputFC = arcpy.GetParameterAsText(0)
    outputFC = arcpy.GetParameterAsText(1)
    heightfield = arcpy.GetParameterAsText(2) #Must be in the same units as the coordinate system!
    azimuth = math.radians(float(arcpy.GetParameterAsText(3))) #Must be in degrees
    altitude = math.radians(float(arcpy.GetParameterAsText(4))) #Must be in degrees

    # Specify output field name for the original FID
    origfidfield = "ORIG_FID"

    # Get field names
    desc = arcpy.Describe(inputFC)
    shapefield = desc.shapeFieldName
    oidfield = desc.oidFieldName

    #Output
    message("Creating output feature class %s ..." % outputFC)
    arcpy.CreateFeatureclass_management(
        os.path.dirname(outputFC),
        os.path.basename(outputFC),
        'POLYGON', "", "", "",
        desc.spatialReference if not arcpy.env.outputCoordinateSystem else "")
    arcpy.AddField_management(outputFC, origfidfield, "LONG")
    inscur = arcpy.InsertCursor(outputFC)

    # Compute the shadow offsets.
    spread = 1/math.tan(altitude) #outside loop as it only needs calculating once

    count = int(arcpy.GetCount_management(inputFC).getOutput(0))
    message("Total features to process: %d" % count)

    searchFields = ",".join([heightfield, oidfield, shapefield])
    rows = arcpy.SearchCursor(inputFC, "", "", searchFields)

    interval = int(count/10.0) # Interval for reporting progress every 10% of rows

    # Create array for holding shadow polygon vertices
    arr = arcpy.Array()
    for r, row in enumerate(rows):
        pctComplete = int(round(float(r) / float(count) * 100.0))
        if r % interval == 0:
            message("%d%% complete" % pctComplete)
        oid = row.getValue(oidfield)
        shape = row.getValue(shapefield)
        height = float(row.getValue(heightfield))

        # Compute the shadow offsets.
        x = -height * spread * math.sin(azimuth)
        y = -height * spread * math.cos(azimuth)

        # Initialize a list of shadow polygons with the original shape as the first
        shadowpolys = [shape]

        # Compute the wall shadows and append them to the list
        for part in shape:
            for i,j in enumerate(range(1,part.count)):
                pnt0 = part[i]
                pnt1 = part[j]
                if pnt0 is None or pnt1 is None:
                    continue # skip null points so that inner wall shadows can also be computed

                # Compute the shadow offset points
                pnt0offset = arcpy.Point(pnt0.X+x,pnt0.Y+y)
                pnt1offset = arcpy.Point(pnt1.X+x,pnt1.Y+y)

                # Construct the shadow polygon and append it to the list
                [arr.add(pnt) for pnt in [pnt0,pnt1,pnt1offset,pnt0offset,pnt0]]
                shadowpolys.append(arcpy.Polygon(arr))
                arr.removeAll() # Clear the array so it can be reused

        # Dissolve the shadow polygons
        dissolved = arcpy.Dissolve_management(shadowpolys,arcpy.Geometry())[0]

        # Insert the dissolved feature into the output feature class
        newrow = inscur.newRow()
        newrow.setValue(origfidfield, oid) # Copy the original FID value to the new feature
        newrow.shape = dissolved
        inscur.insertRow(newrow)

        del row, newrow, shadowpolys, dissolved
    del inscur, rows

if __name__ == "__main__":
    main()
Benutzeravatar
noisefloor
User
Beiträge: 3843
Registriert: Mittwoch 17. Oktober 2007, 21:40
Wohnort: WW
Kontaktdaten:

Hi,

well, in your initial post the actual problem respectively the question is missing.

Did you run the code you posted? If so, what's the result? Any error messages? Did your computer explode (probably and hopefully not ;-) )? ...?

Besides that, the language of this board is German. Basically it's not forbidden to post requests in other languages (AFAIK), but the number of answers you get me be limited.

Regards, noisefloor
BlackJack

@100ri: I think more limiting than the language is using ArcGIS. I'm pretty confident more users here have no problem reading and answering questions in English, than ones who have experience with ArgGIS.

Any general programming (with Python) or Python specific questions are likely to get you answers, but if it involves ArgGIS specific stuff chances of answers are higher in a ArcGIS related forum or community. Like GIS StackExchange where you found that solution.

Regarding the correctness: That's a really hard question to answer. Proving a programm correct is difficult, often even if all the external functionality, in this case `arcpy`, is known.

It compiles. The answer on GIS StackExchange got a 100 points bounty award. There's no comment to that answer saying it doesn't work. So give it a try.

From a Python programmer's point of view there are a couple of ”bad” things in that program:

”Naked” ``except:``\s without one or more specific expected Exception is a bad idea. The handling part is just a ``pass`` so this is a sink of errors. They just get ignored, even ones you might not want to ignore.

The ``del`` statements are all unnecessary. Even worse: they might mislead a reader not fluent in Python into thinking there are objects deleted and memory freed by those statements. ``del`` on bare names deletes those *names*, not the *objects* that are assigned to the names. Such an object may or may not be freed after a ``del`` of the (only) name referencing it. There's no guarantees as to when or even if that happens. All deleted names in the loop are reassigned anyway and deleting names at the end of a function when they fall out of scope anyway just doesn't make any sense at all.

The best way of ”managing memory” in Python is to write functions and limit the scope of names this way. Everything not bound to a name or referenced from a data structure is fair game to the garbage collector.

List comprehensions (LCs) are not meant to write an ordinary loop, used for side effects, in one line. Additionally the abuse of that language feature would be even unnecessary if one wants to write the loop in a single line because this would be possible without the LC syntax.

The `main()` function is much too long. Over 30 local names are used there. That's more than the average person can keep in mind at once. It's much easier to read, reason and talk about a couple of simple, well named functions that are doing a small subtask.

Well named is something that's also important for other objects besides functions. Try to follow the Style Guide for Python Code. Avoid cryptic abbrevations and smushed together words without underscores. `inscur` is an example of both points. Better would be `insert_cursor`. Don't think about saving keystrokes when writing but about readability. Source code is much more often read than written. :-)

The `severity` argument of `message()` isn't used, so the function could be written a bit simpler.

Comments explaining what a function is doing are better set as DocStrings instead of just comments.

Getting at `pnt0` and `pnt1` is a bit weird. Instead of a simple counter and accessing via the counter and the counter + 1, or using some elegant `itertools` solution, there is this hybrid of iterators and old school indexing.

`pctComplete` is computed for each row but used only after a batch worth 10% of rows has been processed.

I personally would avoid ``continue`` as much as possible. It's an unconditional jump to the loop start not reflected in the code structure. And it may make it harder to factor out parts or the whole loop body into function(s). In the given program it would be very easy to change the code to get rid of the ``continue``.
100ri
User
Beiträge: 4
Registriert: Donnerstag 8. September 2016, 18:29

noisefloor hat geschrieben:Hi,

well, in your initial post the actual problem respectively the question is missing.

Did you run the code you posted? If so, what's the result? Any error messages? Did your computer explode (probably and hopefully not ;-) )? ...?

Besides that, the language of this board is German. Basically it's not forbidden to post requests in other languages (AFAIK), but the number of answers you get me be limited.

Regards, noisefloor
Thank you very much for your response. I understand that my options will be limited, but studying in Germany I though posting my question here as well :) The initial question as mentioned, would have been, to see if the script is good and well written or not...
100ri
User
Beiträge: 4
Registriert: Donnerstag 8. September 2016, 18:29

BlackJack hat geschrieben:@100ri: I think more limiting than the language is using ArcGIS. I'm pretty confident more users here have no problem reading and answering questions in English, than ones who have experience with ArgGIS.

Any general programming (with Python) or Python specific questions are likely to get you answers, but if it involves ArgGIS specific stuff chances of answers are higher in a ArcGIS related forum or community. Like GIS StackExchange where you found that solution.

Regarding the correctness: That's a really hard question to answer. Proving a programm correct is difficult, often even if all the external functionality, in this case `arcpy`, is known.

It compiles. The answer on GIS StackExchange got a 100 points bounty award. There's no comment to that answer saying it doesn't work. So give it a try.

From a Python programmer's point of view there are a couple of ”bad” things in that program:

”Naked” ``except:``\s without one or more specific expected Exception is a bad idea. The handling part is just a ``pass`` so this is a sink of errors. They just get ignored, even ones you might not want to ignore.

The ``del`` statements are all unnecessary. Even worse: they might mislead a reader not fluent in Python into thinking there are objects deleted and memory freed by those statements. ``del`` on bare names deletes those *names*, not the *objects* that are assigned to the names. Such an object may or may not be freed after a ``del`` of the (only) name referencing it. There's no guarantees as to when or even if that happens. All deleted names in the loop are reassigned anyway and deleting names at the end of a function when they fall out of scope anyway just doesn't make any sense at all.

The best way of ”managing memory” in Python is to write functions and limit the scope of names this way. Everything not bound to a name or referenced from a data structure is fair game to the garbage collector.

List comprehensions (LCs) are not meant to write an ordinary loop, used for side effects, in one line. Additionally the abuse of that language feature would be even unnecessary if one wants to write the loop in a single line because this would be possible without the LC syntax.

The `main()` function is much too long. Over 30 local names are used there. That's more than the average person can keep in mind at once. It's much easier to read, reason and talk about a couple of simple, well named functions that are doing a small subtask.

Well named is something that's also important for other objects besides functions. Try to follow the Style Guide for Python Code. Avoid cryptic abbrevations and smushed together words without underscores. `inscur` is an example of both points. Better would be `insert_cursor`. Don't think about saving keystrokes when writing but about readability. Source code is much more often read than written. :-)

The `severity` argument of `message()` isn't used, so the function could be written a bit simpler.

Comments explaining what a function is doing are better set as DocStrings instead of just comments.

Getting at `pnt0` and `pnt1` is a bit weird. Instead of a simple counter and accessing via the counter and the counter + 1, or using some elegant `itertools` solution, there is this hybrid of iterators and old school indexing.

`pctComplete` is computed for each row but used only after a batch worth 10% of rows has been processed.

I personally would avoid ``continue`` as much as possible. It's an unconditional jump to the loop start not reflected in the code structure. And it may make it harder to factor out parts or the whole loop body into function(s). In the given program it would be very easy to change the code to get rid of the ``continue``.
Thanks A lot indeed for your complete response. I don't know if you have ever been eager to proceed with something, but missing the tool to go further...? Well currently I'm at that place :D at the moment there is not source for me other than the code which I found by luck on the web to proceed with my research. And I understand your perceptive point of view regarding ArcGIS.
As I have, unfortunately, not background on scripting at all, I was wondering if I could ask you kindly to apply your comments with the code I have posted? Then at least I have the chance to compare the output of the revised code as well and integrate the best result. I appreciate if that will be an option... However I understand one's limit as well. Cheers
BlackJack

@100ri: I don't think its a good idea to change the code unless you've confirmed that it actually solves your problem. And even then the points I've made don't change the functionality but ”just” the quality of the code regarding readability and maintainability. Unless you want to work on the program, which requires programming skills that you don't have as you said, it doesn't make much sense to refactor the code.
100ri
User
Beiträge: 4
Registriert: Donnerstag 8. September 2016, 18:29

BlackJack hat geschrieben:@100ri: I don't think its a good idea to change the code unless you've confirmed that it actually solves your problem. And even then the points I've made don't change the functionality but ”just” the quality of the code regarding readability and maintainability. Unless you want to work on the program, which requires programming skills that you don't have as you said, it doesn't make much sense to refactor the code.


Vielen vielen dank für Ihren Antwort anyways :)
BlackJack

FYI: Packt offers the e-book „Programming ArcGIS 10.1 with Python Cookbook“ (PDF/EPub/Mobi) for the next 11 hours for free: https://www.packtpub.com/packt/offers/free-learning
Antworten