brainteaser viewshed wind turbines - the more you see the worse it gets...

[Anonymous User]

Original User: firefry82

working with arcgis 10.1 Sp1

Hi everybody,
I have a little brainteaser for you: I am trying to make a visibility analysis of wind turbines.  But instead of the simple viewshed thing,  I am trying to include stuff like partial visibility and multiple visibilities. The latter one is driving me insane and I hope some of you guys (or girls of course) maybe have an idea how to approach this.
I have zones around each single turbine, which indicate the level of visual distress in relation to the distance.
[ATTACH=CONFIG]26444[/ATTACH]
I can now merge them to a big viewshed which works fine.

To include the multiple visibilities I need to consider the folowing idea from "Paul et al. 2004"

[ATTACH=CONFIG]26442[/ATTACH]
Paul et al 2004 & ECOGIS
Even though this is for tension towers I think it is about the same for wind turbines. It says that the visibility  increases when more when multiple turbines are visible. So he set up a factor which is 1 for the closest turbine, 1/2 for the second closest,  1/3 for the third closest and so on (1/i). Now how do I find out which turbine is closest to each cell from my raster? and then how do I get a hold of the zone-id?


This is the idea sketch:

[ATTACH=CONFIG]26447[/ATTACH]


What I tried so far:
1) converted viewsheds to polygones, made union analysis to try to intersect each one with the other to find overlapping areas. Then could go through each row and find the biggest zone-id which would be taken in account by the factor 1, the second biggest zone-id with 1/2 and so on. --> failed because union analysis can'T handle the data load

2) tried intersect the the polygon-viewsheds without  the union,  but then only a small part in which they ALL intersect is produced. I the searched for a intersect method that maybe finds and visualizes all the intersections from the overlaps...but I didn't find anything

3) tried to work it out with the raster calculator, but then I wouldn't know how to tell him that he should only take the zone-value when it is the closest to the cell...

I just can't get my head around it...

I am open for any ideas, phyton-scripts, ...whatever...

[Anonymous User]

Original User: firefry82

@csny490
Wow! It's really scary to read this, because I have done this almost the EXACT same way! Exept I used an addition where you suggested the combine-tool. I've pretty much done it in python. here are some screenshots:

the single viewsheds
[ATTACH=CONFIG]26549[/ATTACH]

the zoned-viewsheds with zeroes and big extent (btw. I don't understand what would you use the snap extent for?)
[ATTACH=CONFIG]26550[/ATTACH]

here is the python-code I used to sum (addition) the rasters together

rasters2 = arcpy.ListRasters("Raster*")
for raster in rasters:
    print "processing raster:"+raster
    #big extent
    arcpy.env.extent = "C:/studium/00_Master/clean/clipper.shp"
    #convert nodata to zero
    out1 = Con(IsNull(raster), 0, raster)

    #sum rasters together
    if i == 0:
        out2 = out1
        i += 1
    else:
        out2 = out2 + out1
        i += 1
head, tail = ntpath.split(raster)
#save final output
out2.save("C:/studium/project/script6/output.tif")

*the head and tail are for naming purposes when this script is looped.

after this I used another con to only put the values from the just generated (big) raster, which belong to the viewshed of turbine1 and so on. I have to do this because, I used the near tool to dtermine the distances from all 155 turbines to turbine 1, then 2 and so on. I sorted them asending by distance and numbered them (for the division). The only problem is: this order of appearance (which turbine comes first and which second...) is only correct for turbine1! So I do have to calculate this for all 155 turbines (155*155). In addition I need to erase the parts from the second Raster which I already calculated for Raster1. With my resolution (cell size=5 meters) a single run takes 18 minutes (on a i7 3.5GHZ with 16GB RAM and SSD). But since you would use the same method I think it won't work faster. This is why I was looking for a different solution with that euclidean distance and all.

@msayler : I don't know these tools, let me have a look at them and get back to you asap

[ChrisSnyder]

How about something like:

arcpy.env.extent = "C:/studium/00_Master/clean/clipper.shp"
comboGrdList = []
rastersList = arcpy.ListRasters("Raster*")
for raster in rasterList:
    out1 = Con(IsNull(raster), 0, raster)
    out1.save("twr_" + str(raster.split("_")[-1])) #assuming that raster is something like "tower_12"
    comboGrdList.append(out1)
comboGrd = arcpy.sa.Combine(comboGrdList)

I don't really understand the tower distance order weighting you are trying to implement, but the GenerateNearTable tool will give you a matrix for all 155 ^ 2 combinations and the distances between them. A dictionary that stores the distances (generated from the near table) plus an update cursor on the combine table that accesses the dictionary should be what you need I think.

[Anonymous User]

Original User: msayler

I'm wondering, do you really need to know the actual next towers out? I'm wondering if you could do something like:

1. make a raster where each cell is a count of the towers visible




• create viewsheds for each tower, make 0=not visible and 1=visible (can't remember if this is what the output looks like by default)


• add up the rasters per csny490's instructions, which I believe keeps the cells lined up?

create a heat map or the like (I think that's basically what you have)

normalize the values for the heat map so that they range from 0 to 1

multiply the rasters together (again, some processing may need to be done to make sure cells are coincident)

[nb]

hello again, well I tried everything for the last few days but still no good results...

the combogrid

comboGrd = arcpy.sa.Combine(comboGrdList)

sounded good but it seems that combine tool doesn't work for more than 20 files, because the algorithm writes the extra info in the raster bands and 20 seems to be the maximum. 

create a heat map or the like (I think that's basically what you have)
normalize the values for the heat map so that they range from 0 to 1

I tried that (not sure if I did it right), but with the multiplication small values (e.g. 0,5) keep shrinking instead of slowly growing.

I tried to resample my data to 25mx25 m then cut out the cells where the value is 0 (so no turbine visible) and then converted raster to points (see pciture1).
[ATTACH=CONFIG]26655[/ATTACH]
picture 1

I thought I could then run the "generate near table tool" on the point shapefile and the turbines shapefile (155 points or "turbines") . It took a few minutes, but actually gnereated the table. But I figure running any calculation on
this huge table (1,9 million cells * 155 locations / turbines) will take like foreever

so I am looking for a different solution on how I could accomplish this here:
[ATTACH=CONFIG]26656[/ATTACH]
picture 2: red dots are locations, grid represents (49) cells, green lines are to show the distance. Of course this is just a small scratch. Real values should be 1,9 million cells and 155 locations...


does anybody think this is even doable or do you think its time to give up on it?

[Anonymous User · Accepted solution]

Original User: csny490

"If something is hard you should just give up, since it probably wasn't worth doing it in the first place."

NOT!

Question: What if two towers were visible: One at 200m distance, and the other at 201m distance (shouldn't that distance rank variable be a bit more scaler?)

Okay - How about this (and BTW, this would go way faster if you can coarsen the cell size). Also, make sure the cells are all alligned by setting a snap raster. Also, be aware you might fill up your RAM depending on how many pixels are involved:

1. For each tower make a viewshed/euclidean distance raster so that the output raster name has the tower id in it "tower_12", and the raster values are the distances away from the tower.
2. Convert that raster to a point featureclass.
3. Using a search cursor, load all the point FCs into a Python dictionary. The dictionary key will be a tuple of the X,Y coordinate pair.

So maybe some thing like this (warning UNTESTED!):

pointDict = {} pointFcList = ["tower_1.shp","tower_2.shp","tower_3.shp"] for pointFc in pointfcList:    searchRows = arcpy.da.SearchCursor(pointFc, ["SHAPE@XY","DISTANCE"])    towerId = pointFc.split("_")[-1]    for searchRow in searchRows:       xyKey, distance = searchRow       if xyKey not in pointDict:          pointDict[xyKey] = [(distance, towerId)]       else:          pointDict[xyKey].append((distance, towerId)) for keyKey in pointDict:    pointDict[xyKey].sort()

So now (I think) you have everything you need to compute the score:

How many towers are visible at a given x,y coordinate? len(pointDict[xyKey])
What is the closest tower's distance to a given x,y coordinate? pointDict[xyKey][0][0]
What is the closest tower id to a given x,y coordinate? pointDict[xyKey][0][1]
How about the 2nd closest tower distance? Or the 3rd? pointDict[xyKey][1][0], pointDict[xyKey][2][0]
How about the 2nd closest tower id? Or the 3rd? pointDict[xyKey][1][1], pointDict[xyKey][2][1]

[nb]

Hey csny490, fist of all a big thank you! For the magnificent piece of code which you (if it was nothing) pulled out of your sleeve. When I try to code something its always more like a try and error thing. And then of course for the motivation. I was about to let it all go.
So here is what I did (like you advised is)
1. Set the 0-cells in my vewsheds to nodata, then resample (for testing purposes to 90x90m, 25x25m, 5x5m)
2. create euclidean distance in the same resolutions (90, 25 ,5)
3. Use con to insert the eucl distances in my modified viewshed raster
4. converted these to point shapes (snap to raster, very important!! otherwise the lenght of the dict will always be 1 and the dict will have multiple times more points
5. ran your script with these minor modifications (typos):

for pointFc in pointfcList:

must be pointFcList

for keyKey in pointDict:
   pointDict[xyKey].sort()

must be pointDict[keyKey].sort()

towerId = pointFc.split("_")[-1]

I added [:-4] to cut off the ".shp"

and I have to say it works like a charm! resolution 90x90m takes only 1:30 minutes, 25x25 only 18 minutes. Sadly the 5x5 resolution takes up a sweet time and diskspace. But I guess that would have happend with any solution.

Well I think I still owe you an answer:

Question: What if two towers were visible: One at 200m distance, and the other at 201m distance (shouldn't that distance rank variable be a bit more scaler?)

No, I think. lets assume tower 1 (turbine or whatever) is in zone 6, and tower 2 in zone 5. Then the calculation would be 6 + 5/2= 8,5. And since if you have one huge wind turbine in front of you, one more, a little further away won't make too much difference.
I thought about making my zones more as you say "scale" by using the distance but I think its not worth it.

[Anonymous User]

Original User: csny490

Glad it worked! I replied to your message too - but yeah, go right ahead!