From eac094d0aa93bee2fec916a7c75788e7ef4bf790 Mon Sep 17 00:00:00 2001 From: Santeri Hiitola Date: Thu, 16 May 2024 08:46:25 +0300 Subject: [PATCH] Initial commit. Added command line program to calculate built area from landuse raster and zones --- Scripts/landuse/built_area.py | 103 ++++++++++++++++++++++++++++++ Scripts/landuse/cli_built_area.py | 41 ++++++++++++ 2 files changed, 144 insertions(+) create mode 100644 Scripts/landuse/built_area.py create mode 100644 Scripts/landuse/cli_built_area.py diff --git a/Scripts/landuse/built_area.py b/Scripts/landuse/built_area.py new file mode 100644 index 0000000..f328774 --- /dev/null +++ b/Scripts/landuse/built_area.py @@ -0,0 +1,103 @@ +import rasterstats +import geopandas as gpd +import pandas as pd + +from pathlib import Path + +class BuiltArea: + """ + A class to recalculate the built area of assignment zones based on landcover data and zone borders. + + ... + + Attributes + ---------- + landuse_filepath : str + Filepath of the original .lnd file as a string. + landcover_filepath : str + Filepath of the landcover file to be used in calculating the built area as a string. + In projects with new land use development, the landcover file can be modified. + CRS is expected to be EPSG 3067. + zones_filepath : str + Filepath of the shapefile representing the zones as a string, modified according to + the needs of the project. + CRS is expected to be EPSG 3879. + area_changes : dict + A dictionary describing zone area changes, where keys are original zone IDs, + and values are lists of new zone IDs. + + E.g., {292: [292, 295]} indicates that zone 292 has been split into two, + with one half retaining the original ID and the other assigned a new ID + that corresponds to a new centroid in EMME, following the numbering conventions + described in helmet-docs. + + Methods + ------- + calculate(year, output_path=""): + Calculates the share of built area for each assignment zone. + + Parameters: + year : int + Year of the examination in question. Used in naming the output .lnd file. + output_path : str, optional + Path of the output .lnd file as a string. If not provided, a default path is used. + + """ + + def __init__(self, landuse, landcover, zones, year, area_changes): + self.landuse = Path(landuse) + self.landcover = Path(landcover) + self.zones = Path(zones) + self.year = year + self.area_changes = area_changes + + def calculate(self, output_path=""): + # Check if file paths exist + if not self.landuse.exists() or not self.landcover.exists() or not self.zones.exists(): + print("One or more file paths do not exist.") + return + + print(f"Landuse filepath: {self.landuse}") + print(f"Landcover filepath: {self.landcover}") + print(f"Zones filepath: {self.zones}") + print(f"Target year: {self.year}") + print(f"Area changes received: {self.area_changes}") + + + output = output_path + str(self.year) + ".lnd" + + original_lnd = pd.read_csv(self.landuse, sep="\t", comment="#", index_col=0) + + sijoittelualueet = gpd.read_file(self.zones) + sijoittelualueet.crs = 'EPSG:3879' + sijoittelualueet = sijoittelualueet.to_crs('EPSG:3067') + stats = rasterstats.zonal_stats(sijoittelualueet.geometry, self.landcover, categorical=True) + + df = pd.DataFrame(data=stats) + + df = df.fillna(0) + df.index = sijoittelualueet['SIJ2019'].astype(int) + df.index.name = None + + # Corine ids representing built area, then convert squares to km^2 + df['builtar'] = (df[1] + df[2] + df[3] + df[4]+ df[5] + df[6] + df[7] + df[8] + df[9] + df[10] + df[11] + df[12] + df[13] + df[14] + df[15] + df[16])*0.0004 + df['detach'] = original_lnd['detach'] + df = df.sort_index() + + area_changes_mapped = {} + area_changes_mapped = {i: int(original) for original, new in self.area_changes.items() for i in new} + + df = self._calculate_detach_share_for_region(df, area_changes_mapped, original_lnd) + df = df[['builtar','detach']] + + f = open(output, 'a') + f.write('# Land use 2023\n#\n# builtar: area of built environment\n# detach: detached houses as share of total number of houses\n#\n') + + df.to_csv(f, float_format='%.4g', sep="\t", line_terminator="\n") + f.close() + return + + def _calculate_detach_share_for_region(self, df, area_changes_mapped, original_lnd): + for i in area_changes_mapped.keys(): + df.loc[i]['detach'] = original_lnd.loc[area_changes_mapped[i]]['detach'] + return df \ No newline at end of file diff --git a/Scripts/landuse/cli_built_area.py b/Scripts/landuse/cli_built_area.py new file mode 100644 index 0000000..7d33956 --- /dev/null +++ b/Scripts/landuse/cli_built_area.py @@ -0,0 +1,41 @@ +import argparse +import json +from built_area import BuiltArea +from pathlib import Path + +def get_data(): + parser = argparse.ArgumentParser(description='Recalculate built area of zones.') + parser.add_argument('-a', '--area_changes', type=json.loads, help='dictionary describing area changes') + parser.add_argument('-l', '--landuse', help='filepath of the original .lnd file') + parser.add_argument('-c', '--landcover', help='filepath of the landcover file to be used in calculating the built area') + parser.add_argument('-z', '--zones', help='filepath of the shapefile representing the zones') + parser.add_argument('-y', '--year', type=int, help='year of the examination in question', default=2024) + args = parser.parse_args() + + # Initialize a dictionary to hold area changes + area_changes = args.area_changes if args.area_changes else {} + + # Gather file paths from arguments or user input + landuse = Path(args.landuse) if args.landuse else Path(input("Enter filepath of the original .lnd file: ")) + landcover = Path(args.landcover) if args.landcover else Path(input("Enter filepath of the landcover file to be used in calculating the built area: ")) + zones = Path(args.zones) if args.zones else Path(input("Enter filepath of the shapefile representing the zones: ")) + + # If area_changes is not provided, prompt the user to input them + if not area_changes: + while True: + key = input("Enter original zone id: ") + values = input("Enter new zone ids separated by commas: ") + area_changes[int(key)] = [int(v) for v in values.split(',')] + if input("Do you want to add more area_changes? (y/n): ").lower() != 'y': + break + + # Return all gathered data + return area_changes, landuse, landcover, zones, args.year + +def main(): + area_changes, landuse, landcover, zones, year = get_data() + built_area = BuiltArea(landuse, landcover, zones, year, area_changes) + built_area.calculate() + +if __name__ == '__main__': + main() \ No newline at end of file