diff --git a/NATURF/Binary.py b/NATURF/Binary.py index 6fd7c982fb..24413c4ce3 100644 --- a/NATURF/Binary.py +++ b/NATURF/Binary.py @@ -1,7 +1,7 @@ from sys import exit from math import ceil from numpy import array, zeros, uint8, uint16, indices,empty, float64, \ - isnan, ones, save, float32, double + isnan, ones, save, float32, double, amax from skimage.draw import polygon import osr import pickle @@ -344,7 +344,7 @@ def make_binary(name, path, tifDir, cenlat, cenlon, outputDir, id_field, height_ for i in range(len(fad_out2)): for j in range(len(bldgsh)): dbh[j][bldgsh[j] == double(names2[i])] += float(1) / float(len(fad_out2)) - + print(amax(dbh)) master = zeros((132, IMAGE_SIZE_X, IMAGE_SIZE_Y), dtype=float64) master[0:15] = fadn[0:15] @@ -453,11 +453,11 @@ def make_binary(name, path, tifDir, cenlat, cenlon, outputDir, id_field, height_ index.write(' dx=100.0\n') index.write(' known_x=1\n') index.write(' known_y=1\n') - index.write(' known_lat=38.794034\n') - index.write(' known_lon=-77.129474\n') + index.write(' known_lat=35.6936372\n') + index.write(' known_lon=-115.4739874\n') index.write(' truelat1=45.5\n') index.write(' truelat2=29.5\n') - index.write(' stdlon=-77.0\n') + index.write(' stdlon=-114.75\n') index.write(' wordsize=4\n') index.write(' endian=big\n') index.write(' signed=no\n') diff --git a/NATURF/Parameter_Calculations.py b/NATURF/Parameter_Calculations.py index 0e9545b6ff..f3a5d8b5da 100644 --- a/NATURF/Parameter_Calculations.py +++ b/NATURF/Parameter_Calculations.py @@ -2,7 +2,7 @@ from Parameter_Definitions import * from sys import exit import math -from numpy import array, zeros, uint8, mean +from numpy import array, zeros, uint8, mean, set_printoptions, inf from PIL import Image import csv import warnings @@ -84,8 +84,10 @@ def calculate_parameters(path, tifDir, outputDir, name, height_field, id_field): buils_ids = tifDir +'/Building_IDs.tif' ids = Image.open(buils_ids) - + #print(ids) ids = array(ids) + #set_printoptions(threshold=inf) + #print(ids) filename = "%s.csv" % name @@ -143,13 +145,15 @@ def calculate_parameters(path, tifDir, outputDir, name, height_field, id_field): bs21 = mean(bs2par_out[i]) # 1 for qw in range(15): - + #print("Before", fad_out[i]['n'][qw]) fad_out[i]['n'][qw] = mean(fad_out[i]['n'][qw]) + #print("After", fad_out[i]['n'][qw]) fad_out[i]['w'][qw] = mean(fad_out[i]['w'][qw]) fad_out[i]['s'][qw] = mean(fad_out[i]['s'][qw]) fad_out[i]['e'][qw] = mean(fad_out[i]['e'][qw]) fad1 = fad_out[i] # 15 * 4 + #print(fad1) fai_out[i]['n'] = mean(fai_out[i]['n']) fai_out[i]['w'] = mean(fai_out[i]['w']) @@ -175,8 +179,9 @@ def calculate_parameters(path, tifDir, outputDir, name, height_field, id_field): rrl1 = rrl_out[i] # 4 rdh1 = rdh_out[i] # 4 mrl1 = mrl_out[i] # 4 - - paf = mean(builfrac_out[i]) + #print(builfrac_out[i]) + paf = mean(builfrac_out[i]) # Averages the building fraction over all vertical heights + #print(paf) mdh1 = mean(mdh_out[i]) # 1? diff --git a/NATURF/Parameter_Definitions.py b/NATURF/Parameter_Definitions.py index 8642bb1eb0..9eec813c7c 100644 --- a/NATURF/Parameter_Definitions.py +++ b/NATURF/Parameter_Definitions.py @@ -1,67 +1,67 @@ import math -from numpy import array, zeros, uint8, uint16, where, empty, sum, unique, rint +from numpy import array, zeros, uint8, uint16, where, empty, sum, unique, rint, array_equiv from skimage.draw import polygon from scipy import ndimage def get_pixels2_c(IMAGE_SIZE_X, IMAGE_SIZE_Y, x1, y1, x2, y2): - mx = (float(y1) - float(y2)) / (float(x1) - float(x2) + 0.0000000000000000000000000000001) - yint = float(y1) - (mx * float(x1)) + mx = (float(y1) - float(y2)) / (float(x1) - float(x2) + 0.0000000000000000000000000000001) # Finds the slope between two points + yint = float(y1) - (mx * float(x1)) # Finds the slope of the line px_func = [] py_func = [] - if abs(x1 - x2) > abs(y1 - y2) or abs(x1 - x2) == abs(y1 - y2): + if abs(x1 - x2) > abs(y1 - y2) or abs(x1 - x2) == abs(y1 - y2): # If the absolute value of the slope is >= 1 if x1 > x2: - for i in range(x2, x1 + 1): - yp = int(round(mx * i + yint)) + for i in range(x2, x1 + 1): # All x integers between x1 and x2 + yp = int(round(mx * i + yint)) # The y integers for each x integer - bx = i - by = yp + bx = i # The current x coordinate + by = yp # The current y coordinate - if 0 <= bx < IMAGE_SIZE_Y and 0 <= by < IMAGE_SIZE_X: - px_func.append(i) - py_func.append(yp) + if 0 <= bx < IMAGE_SIZE_Y and 0 <= by < IMAGE_SIZE_X: # Checks if the current coordinate is in quad 1 and within the tile + px_func.append(i) # Appends the x coordinate + py_func.append(yp) # Appends the y coordinate if x1 < x2: - for i in range(x1, x2 + 1): - yp = int(round(mx * i + yint)) + for i in range(x1, x2 + 1): # All x integers between x1 and x2 + yp = int(round(mx * i + yint)) # the y integers for the x integers - bx = i - by = yp + bx = i # The current x coordinate + by = yp # The current y coordinate - if 0 <= bx < IMAGE_SIZE_Y and 0 <= by < IMAGE_SIZE_X: - px_func.append(i) - py_func.append(yp) + if 0 <= bx < IMAGE_SIZE_Y and 0 <= by < IMAGE_SIZE_X: # Checks if the current coordinate is in quad 1 and within the tile + px_func.append(i) # Appends the x coordinate + py_func.append(yp) # Appends the y coordinate - if abs(x1 - x2) < abs(y1 - y2): + if abs(x1 - x2) < abs(y1 - y2): # If the absolute value of the slope is > 1 if y1 > y2: - for i in range(y2, y1 + 1): - xp = int(round((i - yint) / mx)) + for i in range(y2, y1 + 1): # All y integers between y1 and y2 + xp = int(round((i - yint) / mx)) # The x integers for the y integers - bx = xp - by = i + bx = xp # The current x coordinate + by = i # The current y coordinate - if 0 <= bx < IMAGE_SIZE_Y and 0 <= by < IMAGE_SIZE_X: + if 0 <= bx < IMAGE_SIZE_Y and 0 <= by < IMAGE_SIZE_X: # Checks if the current coordinate is in quad 1 and within the tile px_func.append(xp) py_func.append(i) if y1 < y2: - for i in range(y1, y2 + 1): - xp = int(round((i - yint) / mx)) + for i in range(y1, y2 + 1): # All y integers between y1 and y2 + xp = int(round((i - yint) / mx)) # The x integers for the y integers - bx = xp - by = i + bx = xp # The current x coordinate + by = i # The current y coordinate - if 0 <= bx < IMAGE_SIZE_Y and 0 <= by < IMAGE_SIZE_X: + if 0 <= bx < IMAGE_SIZE_Y and 0 <= by < IMAGE_SIZE_X: # Checks if the current coordinate is in quad 1 and within the tile px_func.append(xp) py_func.append(i) - return px_func, py_func + return px_func, py_func # Returns the coordinates for all integers within each wall def inter4p(x1, y1, x2, y2, rad1): - dx = x1-x2 - dy = y1-y2 - dist = math.sqrt(dx*dx + dy*dy) - dx /= (dist + 0.00000000000000000000001) - dy /= (dist + 0.00000000000000000000001) + dx = x1-x2 # The distance between the x coordinates + dy = y1-y2 # The distance between the y coordinates + dist = math.sqrt(dx*dx + dy*dy) # The distance between the two points + dx /= (dist + 0.00000000000000000000001) # The ratio of the horizontal distance to the hypotenuse + dy /= (dist + 0.00000000000000000000001) # The ratio of the vertical distance to the hypotenuse if y1 > y2: x1r = x1 - rad1*dy @@ -107,11 +107,11 @@ def ang2lines(linea, lineb): def ang2points(x1, y1, x2, y2): - dx = x2 - x1 - dy = y2 - y1 + dx = x2 - x1 # The distance between the x coordinates + dy = y2 - y1 # The distance between the y coordinates - r = math.atan2(dy, dx) - degr = math.degrees(r) + r = math.atan2(dy, dx) # The angle facing the y coordinate (theta-y) in radians + degr = math.degrees(r) # Theta-y in degrees if degr < 0: degr += 360 @@ -215,39 +215,43 @@ def avg_building_dist(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, ids, PIXEL_SIZE, heig dils.fill(0) narr.fill(0) if len(ids[ids == cntr]) != 0: - i2arr[ids == cntr] = 1 - builarea = sum(i2arr) * PIXEL_SIZE**2 - i2arr = ndimage.binary_dilation(i2arr, structure=struct, iterations=rad).astype(i2arr.dtype) + i2arr[ids == cntr] = 1 # Puts a 1 wherever the current counter corresponds with the IDs file + builarea = sum(i2arr) * PIXEL_SIZE**2 # Sums up the 1s and multiplies by the pixel area to get the current building area + i2arr = ndimage.binary_dilation(i2arr, structure=struct, iterations=rad).astype(i2arr.dtype) # Dilates the i2arr to get a square with radius 100 pixels # dilarea = sum(i2arr) * PIXEL_SIZE**2 - dibuils = unique(ids[where(i2arr == 1)]) - + dibuils = unique(ids[where(i2arr == 1)]) # A list of the unique building IDs found in the current dilated area + #print(dibuils) + #print(heights) sumareas += builarea # print dibuils if len(dibuils) != 0: - currcx = cents[cntr][0] # current building centroid - currcy = cents[cntr][1] + currcx = cents[cntr][0] # current building centroid x coordinate + currcy = cents[cntr][1] # current building centroid y coordinate # print currcx, currcy te, tn, ts, tw = 0., 0., 0., 0. ce, cn, cs, cw = 0., 0., 0., 0. - - for asdf in dibuils: - if asdf != 0 and asdf in heights: - hts1.append(heights[asdf]) - dils[(ids == asdf) & (i2arr == 1)] = 1 - - sumarhts += ars[asdf] * heights[asdf] - sumareas2 += ars[asdf] - - if asdf != cntr: - cx = cents[asdf][0] - cy = cents[asdf][1] - - d = math.hypot((currcx - cx), (currcy - cy)) - angle = ang2points(currcx, currcy, cx, cy) - + + for asdf in dibuils: # loops through every building and calculates the distance and angle from that building to all other buildings + #print("asdf", asdf) + if asdf != 0 and asdf in heights: # Checks that the unique ID is not 0 (because height indexing begins at 1) and that it is in the heights list + hts1.append(heights[asdf]) # Appends the current building height to the hts1 list + dils[(ids == asdf) & (i2arr == 1)] = 1 # Creates an array with 1s where the current ID is within the dilated area + + sumarhts += ars[asdf] * heights[asdf] # Multiplies the current building area and height together, and then adds it to a rolling sum for the tile + sumareas2 += ars[asdf] # Adds the current building height to a rolling sum + #print("cntr", cntr) + if asdf != cntr: # Loops through every building except for the current one + cx = cents[asdf][0] # The centroid x coordinate of the other building + cy = cents[asdf][1] # The centroid y coordinate of the other building + #print("currcx",currcx) + #print("cx", cx) + d = math.hypot((currcx - cx), (currcy - cy)) # The distance between the two centroids + angle = ang2points(currcx, currcy, cx, cy) # The angle for the y coordinate + + # The angles correspond to a circle where 0/360 degrees is directly east, and the degrees increase counter-clockwise if (315 <= angle <= 360 or 0 <= angle < 45) and d != 0: # east te += d ce += 1 @@ -260,27 +264,28 @@ def avg_building_dist(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, ids, PIXEL_SIZE, heig tn += d cn += 1 - if 225 <= angle < 315 and d != 0: # south + if 225 <= angle < 315 and d != 0: # south ts += d cs += 1 if cn != 0 or cs != 0: - avgns = (tn + ts) / (cn + cs) + avgns = (tn + ts) / (cn + cs) # The mean north and south distance from the current building to all other buildings else: avgns = 0 if ce != 0 or cw != 0: - avgew = (te + tw) / (ce + cw) + avgew = (te + tw) / (ce + cw) # The mean west and east distance from the current building to all other buildings else: avgew = 0 # print avgns, avgew - narr[(dils == 1) & (i2arr == 1)] = 1 - parea = sum(narr) * PIXEL_SIZE**2 - - pareas[cntr] = [ht, parea] # ####################################################### + narr[(dils == 1) & (i2arr == 1)] = 1 # Creates an array with 1s only where the buildings are + #print(array_equiv(dils, narr)) + parea = sum(narr) * PIXEL_SIZE**2 # Calculates the area of all buildings within the plan area + pareas[cntr] = [ht, parea] # Creates a dictionary where the counter corresponds to the height and plan area of each building + #print(pareas) # if parea > dilarea: # print parea, dilarea @@ -400,6 +405,7 @@ def parameters1(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, ids, PIXEL_SIZE, height_fiel i2arr = ndimage.binary_dilation(i2arr, structure=struct, iterations=rad).astype(i2arr.dtype) # i2arr[(ids != bid) & (iarr == 127)] = 0 dilarea = sum(i2arr) * PIXEL_SIZE**2 + #print(dilarea) sumarhts += builarea * nht sumareas += builarea @@ -467,6 +473,7 @@ def parameters1(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, ids, PIXEL_SIZE, height_fiel if newbarea[counter][0] > ht: builfrac = newbarea[counter][1] / dilarea + else: builfrac = 0 @@ -512,8 +519,10 @@ def parameters1(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, ids, PIXEL_SIZE, height_fiel if newbarea[counter][0] > ht: builfrac = newbarea[counter][1] / dilarea + else: builfrac = 0 + print(builfrac) if (cents_ns[counter] * cents_ew[counter]) != 0: fai = (breadth * zh) / (cents_ns[counter] * cents_ew[counter]) @@ -602,6 +611,7 @@ def parameters1(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, ids, PIXEL_SIZE, height_fiel if newbarea[counter][0] > ht: builfrac = newbarea[counter][1] / dilarea + #print(builfrac, newbarea[counter][1], dilarea) else: builfrac = 0 @@ -749,13 +759,16 @@ def parameters1(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, ids, PIXEL_SIZE, height_fiel i2arr.fill(0) builarea = 0 dilarea = 0 - # print counter + #print(counter) if len(ids[ids == counter]) != 0: - i2arr[ids == counter] = 1 - builarea = sum(i2arr) * PIXEL_SIZE**2 + i2arr[ids == counter] = 1 # Puts a 1 wherever the current counter corresponds with the IDs file + #print(i2arr) + builarea = sum(i2arr) * PIXEL_SIZE**2 # Sums up the 1s and multiplies by the pixel area to get the current building area # i2arr[(ids != counter) & (iarr == 127)] = 0 - dilarea = sum(i2arr) * PIXEL_SIZE**2 + i2arr = ndimage.binary_dilation(i2arr, structure=struct, iterations=rad).astype(i2arr.dtype) # Dilates the i2arr to get a square with radius 100 pixels + dilarea = sum(i2arr) * PIXEL_SIZE**2 # Sums up all the 1s in i2arr and multiplies by the pixel area to get the dilated area for each building + #print(sum(i2arr)) sumarhts += builarea * nht sumareas += builarea @@ -764,42 +777,45 @@ def parameters1(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, ids, PIXEL_SIZE, height_fiel dilarea = builarea - nx, ny, nz = ring.GetPoint(0) - - parrxc = [[nx]] - parryc = [[ny]] - + nx, ny, nz = ring.GetPoint(0) # The minimum x and y for each building + #print(ring.GetPoint(0)) + #print(ring.GetPointCount()) + parrxc = [[nx]] # Puts the minimum x into a list of lists + parryc = [[ny]] # Puts the minimum y into a list of lists + for bi in range(1, ring.GetPointCount()): - bpt = ring.GetPoint(bi) - - parrxc.append([bpt[0]]) - parryc.append([bpt[1]]) + #print(bi) + bpt = ring.GetPoint(bi) # Finds the coordinates for each point in the building polygon + parrxc.append([bpt[0]]) # Appends the x coordinates to the parrxc list + parryc.append([bpt[1]]) # Appends the y coordinates to the parrxc list + parrxc = array(parrxc) parryc = array(parryc) - + for na in range(0, len(parrxc)-1): - p1x = parrxc[na] - p1y = parryc[na] + p1x = parrxc[na] # Extracts the x coordinate corresponding to na + p1y = parryc[na] # Extracts the y coordinate corresponding to na - p2x = parrxc[na+1] - p2y = parryc[na+1] + p2x = parrxc[na+1] # Extracts the x coordinate corresponding to na + 1 + p2y = parryc[na+1] # Extracts the y coordinate corresponding to na + 1 if (p1x, p1y) != (p2x, p2y): - perpx0 = parrxc[na] - + perpx0 = parrxc[na] # Sets the origin x coordinate to p1x + #print("parryc", parryc[na]) + #print("IMAGE_SIZE_X", IMAGE_SIZE_X) if parryc[na] + 1 < IMAGE_SIZE_X: - perpy0 = parryc[na] + 1 + perpy0 = parryc[na] + 1 # Sets the origin y coordinate to p1y + 1 if it is less than the width of the study area else: - perpy0 = IMAGE_SIZE_X - - line1 = [[p1x, p1y], [p2x, p2y]] - line0deg = [[p1x, p1y], [perpx0, perpy0]] + perpy0 = IMAGE_SIZE_X # Otherwise, sets the origin y coordinate to the width of the study area - deg = ang2lines(line1, line0deg) - - dist = math.sqrt((p2x - p1x)**2 + (p2y - p1y)**2) + line1 = [[p1x, p1y], [p2x, p2y]] # Creates a line from the na point to the na + 1 point + line0deg = [[p1x, p1y], [perpx0, perpy0]] # Creates a line from the na point to + #print(perpy0) + deg = ang2lines(line1, line0deg) # Returns the direction of the wind normal to the wall + #print(deg) + dist = math.sqrt((p2x - p1x)**2 + (p2y - p1y)**2) # The length of the wall breadth = dist @@ -819,8 +835,8 @@ def parameters1(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, ids, PIXEL_SIZE, height_fiel if 315 <= deg <= 360 or 0 <= deg < 45: if ((dist / 10) - 1) > 0: - edist = dist - wallarea = dist * nht + edist = dist # Used for calculation of complete aspect ratio + wallarea = dist * nht # The area of the wall at the current 5m slice if dilarea == 0: break @@ -861,7 +877,7 @@ def parameters1(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, ids, PIXEL_SIZE, height_fiel if 45 <= deg < 135: if ((dist / 10) - 1) > 0: - wallarea = dist * nht + wallarea = dist * nht # The area of the wall at the current 5m slice if asdf == 0: dilarea = 10000 @@ -904,8 +920,8 @@ def parameters1(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, ids, PIXEL_SIZE, height_fiel if 135 <= deg < 225: if ((dist / 10) - 1) > 0: - wdist = dist - wallarea = dist * nht + wdist = dist # Used for calculation of complete aspect ratio + wallarea = dist * nht # The area of the wall at the current 5m slice if dilarea == 0: break @@ -946,8 +962,9 @@ def parameters1(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, ids, PIXEL_SIZE, height_fiel if (direc.lower() == 'south') and (225 <= deg < 315): if ((dist / 10) - 1) > 0: - sdist = dist - wallarea = dist * nht + sdist = dist # Used for calculation of complete aspect ratio + wallarea = dist * nht # The area of the wall at the current 5m slice + if dilarea == 0: break @@ -957,8 +974,9 @@ def parameters1(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, ids, PIXEL_SIZE, height_fiel # print newbarea[counter][0], float(asdf), float(newbarea[counter][0]) > float(asdf) - if float(newbarea[counter][0]) > float(asdf): - builfrac = newbarea[counter][1] / dilarea + if float(newbarea[counter][0]) > float(asdf): # Checks if the height of the current building is greater than the current 5m slice + builfrac = newbarea[counter][1] / dilarea # Divides the plan area by the dilated area to calculate the ratio of building area to dilated area + #print(builfrac, newbarea[counter][1], dilarea) else: builfrac = 0 @@ -999,7 +1017,7 @@ def parameters1(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, ids, PIXEL_SIZE, height_fiel mrls.append(mrl) bs2pars.append(bs2par) - + #print(newbarea[counter][1], dilarea) fad_out['n'].append(fadn) fad_out['w'].append(fadw) fad_out['s'].append(fads) @@ -1044,9 +1062,9 @@ def parameters1(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, ids, PIXEL_SIZE, height_fiel # print sdist, edist, wdist if edist != 0: - car_out_inc.append(float(sdist) / float(edist)) + car_out_inc.append(float(sdist) / float(edist)) # The ratio of the length of the south wall to the east wall elif edist == 0 and wdist != 0: - car_out_inc.append(float(sdist) / float(wdist)) + car_out_inc.append(float(sdist) / float(wdist)) # The ratio of the length of the south wall to the west wall else: car_out_inc.append(1) @@ -1056,6 +1074,7 @@ def parameters1(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, ids, PIXEL_SIZE, height_fiel feature_buil = layer2.GetNextFeature() awmh = sumarhts / sumareas + #print(builfrac_out_inc) return fad_out_inc, builfrac_out_inc, fai_out_inc, rdh_out_inc, rrl_out_inc, mdh_out_inc, mrl_out_inc, \ bs2par_out_inc, zo_out_inc, zd_out_inc, car_out_inc @@ -1081,33 +1100,33 @@ def h2w(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, xOrigin, yOrigin, pixelWidth, pixelH numpoints = ring.GetPointCount() if numpoints != 0: - nx, ny, nz = ring.GetPoint(0) + nx, ny, nz = ring.GetPoint(0) # The minimum x and y coordinates for the current building - xOffset = int((nx - xOrigin) / pixelWidth) - yOffset = int((ny - yOrigin) / pixelHeight) + xOffset = int((nx - xOrigin) / pixelWidth) # The distance from the minimum x coordinate to the left extent of the tile + yOffset = int((ny - yOrigin) / pixelHeight) # The distance from the minimum y coordinate to the top extent of the tile - parr = [[xOffset, yOffset]] + parr = [[xOffset, yOffset]] parrx = [[xOffset]] parry = [[yOffset]] for bi in range(1, ring.GetPointCount()): - bpt = ring.GetPoint(bi) + bpt = ring.GetPoint(bi) # Finds the coordinates for each point in the building polygon - xoff = int((bpt[0] - xOrigin) / pixelWidth) - yoff = int((bpt[1] - yOrigin) / pixelHeight) + xoff = int((bpt[0] - xOrigin) / pixelWidth) # Calculates the xoffset for each point in the building polygon + yoff = int((bpt[1] - yOrigin) / pixelHeight) # Calculates the yoffset for each point in the building polygon parr.append([xoff, yoff]) parrx.append([xoff]) parry.append([yoff]) - parrx[:] = [x[0] - start_x_p for x in parrx] + parrx[:] = [x[0] - start_x_p for x in parrx] parry[:] = [y[0] - start_y_p for y in parry] parrx = array(parrx) parry = array(parry) - rr, cc = polygon(parrx, parry, [IMAGE_SIZE_Y, IMAGE_SIZE_X]) + rr, cc = polygon(parrx, parry, [IMAGE_SIZE_Y, IMAGE_SIZE_X]) # Uses the offsets to creates arrays of polygon coordinates - buils[cc, rr] = 127 + buils[cc, rr] = 127 # Creates an array with 127 where the polygon coordinates dictate feature_buil.Destroy() feature_buil = layer2.GetNextFeature() @@ -1135,10 +1154,10 @@ def h2w(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, xOrigin, yOrigin, pixelWidth, pixelH cntr1 = 0 summ1 = 0 - nx, ny, nz = ring.GetPoint(0) + nx, ny, nz = ring.GetPoint(0) # The minimum x and y coordinates for the current building - xOffset = int((nx - xOrigin) / pixelWidth) - yOffset = int((ny - yOrigin) / pixelHeight) + xOffset = int((nx - xOrigin) / pixelWidth) # The distance from the minimum x coordinate to the left extent of the tile + yOffset = int((ny - yOrigin) / pixelHeight) # The distance from the minimum y coordinate to the top extent of the tile parr = [[xOffset, yOffset]] parrx = [[xOffset]] @@ -1148,10 +1167,10 @@ def h2w(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, xOrigin, yOrigin, pixelWidth, pixelH parryc = [[ny]] for bi in range(1, ring.GetPointCount()): - bpt = ring.GetPoint(bi) + bpt = ring.GetPoint(bi) # Finds the coordinates for each point in the building polygon - xoff = int((bpt[0] - xOrigin) / pixelWidth) - yoff = int((bpt[1] - yOrigin) / pixelHeight) + xoff = int((bpt[0] - xOrigin) / pixelWidth) # Calculates the xoffset for each point in the building polygon + yoff = int((bpt[1] - yOrigin) / pixelHeight) # Calculates the yoffset for each point in the building polygon parr.append([xoff, yoff]) parrx.append([xoff]) @@ -1175,16 +1194,16 @@ def h2w(IMAGE_SIZE_X, IMAGE_SIZE_Y, layer2, xOrigin, yOrigin, pixelWidth, pixelH xarr, yarr = get_pixels2_c(IMAGE_SIZE_X, IMAGE_SIZE_Y, int(p1x), int(p1y), int(p2x), int(p2y)) - if ((len(xarr)/20)-1) > 0: + if ((len(xarr)/20)-1) > 0: # If there are more than 20 points in the wall: - for nb in range(0, int(rint((len(xarr)/20.)-1))): + for nb in range(0, int(rint((len(xarr)/20.)-1))): tempa = 0 - x1 = xarr[(nb*20)+1] - y1 = yarr[(nb*20)+1] + x1 = xarr[(nb*20)+1] # Find a starting x coordinate + y1 = yarr[(nb*20)+1] # Find a starting y coordinate - x2 = xarr[((nb+1)*20)+1] - y2 = yarr[((nb+1)*20)+1] + x2 = xarr[((nb+1)*20)+1] # Find the x coordinate 20 away from the starting coordinate + y2 = yarr[((nb+1)*20)+1] # Find the y coordinate 20 away from the starting coordinate for tg in range(1, 200): perp1x, perp1y = inter4p(x1, y1, x2, y2, tg) diff --git a/NATURF/im3Workflow.py b/NATURF/im3Workflow.py index 17526176bb..e5f5fd9cc0 100644 --- a/NATURF/im3Workflow.py +++ b/NATURF/im3Workflow.py @@ -11,7 +11,7 @@ def processFile(shapeFile): scenario = shapeFile.parent.parent.stem # first part of output directory name buildingTile = shapeFile.stem # filename without extension outputDir = folder/'{}_{}'.format(scenario, buildingTile) # create output directory named according to tile - tifDir = outputDir/'Tifs' + tifDir = outputDir/'{}_{}_Tifs'.format(scenario, buildingTile) tifDir.mkdir(parents=True, exist_ok=True) # automatically also creates outputDir indexfile = filenames[buildingTile][0] # binary file name @@ -27,9 +27,9 @@ def processFile(shapeFile): print(buildingTile + " has finished") -basePath = pathlib.Path('/mnt/data/UrbanMorphology/data/32x32/dilarea_testing/Shapefiles') # Path to directory containing shapefiles +basePath = pathlib.Path('/mnt/data/UrbanMorphology/data/WDC/100m/test') # Path to directory containing shapefiles shapeFiles = list(basePath.glob('*.shp')) # Creates a list of all the shapefile paths -grid = pd.read_csv(pathlib.Path('/mnt/data/UrbanMorphology/data/32x32/LA_Indices.csv')) # Path to geogrid csv +grid = pd.read_csv(pathlib.Path('/mnt/data/UrbanMorphology/data/WDC/100m/DC_Indices.csv')) # Path to geogrid csv filenames = grid.set_index('GRID_ID').T.to_dict("list") # Creates a dictionary from the csv where grid IDs and the names are values try: