|
| 1 | +import json |
| 2 | +import os |
| 3 | +import numpy as np |
| 4 | +import logging |
| 5 | +import rasterio |
| 6 | +import argparse |
| 7 | +from rasterio.transform import rowcol |
| 8 | +from rasterio.windows import Window |
| 9 | +from pyproj import Proj, Transformer |
| 10 | + |
| 11 | +from config import LOG_LEVEL, LOG_FORMAT, LOG_DIR |
| 12 | + |
| 13 | +logging.basicConfig(level=LOG_LEVEL, format=LOG_FORMAT) |
| 14 | + |
| 15 | +logfile = "assign_field_from_raster-log.txt" |
| 16 | + |
| 17 | +def main(): |
| 18 | + |
| 19 | + description_str = """ |
| 20 | + Assigns spatial field to properties, spatial_fields from raster to geojson. If spatial field allready exists, then append values. |
| 21 | + """ |
| 22 | + parser = argparse.ArgumentParser(prog="assign_field_from_raster.py", |
| 23 | + description=description_str) |
| 24 | + parser.add_argument('input_geojson', type=str, |
| 25 | + help='geojson containing features to assign values to.') |
| 26 | + parser.add_argument('raster', type=str, |
| 27 | + help='raster file') |
| 28 | + parser.add_argument('assigned_geojson', type=str, |
| 29 | + help='out file') |
| 30 | + parser.add_argument('field_name', type=str, |
| 31 | + help='Name of the assigned field.') |
| 32 | + parser.add_argument('-c','--categorical', action='store_true', |
| 33 | + help="Raster contains cathegorical value to be assigned to each feature") |
| 34 | + args = parser.parse_args() |
| 35 | + |
| 36 | + # Add new file handler to logger. |
| 37 | + file_handler = logging.FileHandler(filename=os.path.join(LOG_DIR, logfile)) |
| 38 | + log_formatter = logging.Formatter(fmt=LOG_FORMAT) |
| 39 | + file_handler.setFormatter(log_formatter) |
| 40 | + file_handler.setLevel(LOG_LEVEL) |
| 41 | + logging.getLogger().addHandler(file_handler) |
| 42 | + |
| 43 | + with open(args.input_geojson, 'r') as file: |
| 44 | + logging.info("Reads elements from file: {}".format(args.input_geojson)) |
| 45 | + elements = json.load(file) |
| 46 | + |
| 47 | + assigned = { |
| 48 | + "type": "FeatureCollection", |
| 49 | + "features": [], |
| 50 | + } |
| 51 | + |
| 52 | + with rasterio.open(args.raster) as dataset: |
| 53 | + logging.info("Reads data from raster: {}".format(dataset.name)) |
| 54 | + rastercoords_from_lonlat = Transformer.from_proj( |
| 55 | + Proj('epsg:4326'), # source coordinates (lonlat) |
| 56 | + Proj(dataset.crs), # target coordinates |
| 57 | + always_xy=True # Use easting-northing, longitude-latitude order of coordinates. |
| 58 | + ) |
| 59 | + # transforms lat-lon to raster row-col of raster |
| 60 | + # array = dataset.read(masked=True, out_dtype=np.float32) |
| 61 | + |
| 62 | + # rowcol(region_dataset.transform, -9.73363, 36.94755) |
| 63 | + nr_of_assigned_features = 0 |
| 64 | + for feature in elements['features']: |
| 65 | + if nr_of_assigned_features%100 == 0: |
| 66 | + logging.info("Assigned features: {}".format(nr_of_assigned_features)) |
| 67 | + |
| 68 | + coords = feature["geometry"]["coordinates"] |
| 69 | + xs, ys = rastercoords_from_lonlat.transform(*zip(*coords)) |
| 70 | + rows, cols = rowcol(dataset.transform, xs, ys) |
| 71 | + # rows, cols = rowcol(dataset.transform, *zip(*coords)) |
| 72 | + window, window_rows, window_cols = get_window(rows, cols) |
| 73 | + array = dataset.read(out_dtype=np.float64, window=window) |
| 74 | + |
| 75 | + try: |
| 76 | + spatial_field_list = np.round(array[:, window_rows, window_cols], 3).tolist() |
| 77 | + |
| 78 | + except IndexError as error: |
| 79 | + logging.warning("{} - OSM Segment is outside of raster bounds.") |
| 80 | + contained_in_raster = [0 <= row < dataset.shape[0] and 0 <= col < dataset.shape[1] for (row, col) in |
| 81 | + zip(rows, cols)] |
| 82 | + rows = [row for (contained, row) in zip(contained_in_raster, rows) if contained] |
| 83 | + cols = [col for (contained, col) in zip(contained_in_raster, cols) if contained] |
| 84 | + window, window_rows, window_cols = get_window(rows, cols) |
| 85 | + array = dataset.read(out_dtype=np.float64, window=window) |
| 86 | + |
| 87 | + # append zero values outside of raster bounds. |
| 88 | + padded_array = np.zeros([dataset.count, len(contained_in_raster)]) |
| 89 | + padded_array[:, contained_in_raster] = array[:, window_rows, window_cols] |
| 90 | + spatial_field_list = np.round(padded_array, 3).tolist() |
| 91 | + |
| 92 | + if not args.categorical: |
| 93 | + if args.field_name in feature["properties"]["spatial_fields"]: |
| 94 | + # Append to existing values |
| 95 | + feature["properties"]["spatial_fields"][args.field_name].extend(spatial_field_list) |
| 96 | + else: |
| 97 | + # create new property |
| 98 | + feature["properties"]["spatial_fields"][args.field_name] = spatial_field_list |
| 99 | + else: |
| 100 | + # categorical value. Assign most frequent value as property. |
| 101 | + feature["properties"][args.field_name] = int(np.bincount(spatial_field_list[0]).argmax()) |
| 102 | + assigned["features"].append(feature) |
| 103 | + nr_of_assigned_features += 1 |
| 104 | + logging.info("Done processing features. Updated {} features".format(nr_of_assigned_features)) |
| 105 | + with open(args.assigned_geojson, 'w') as outfile: |
| 106 | + json.dump(assigned, outfile) |
| 107 | + logging.info("Wrote to file: {}".format(args.assigned_geojson)) |
| 108 | + |
| 109 | +def get_window(rows, cols): |
| 110 | + # find window |
| 111 | + col_off = min(cols) |
| 112 | + row_off = min(rows) |
| 113 | + width = max(cols) - min(cols) + 1 |
| 114 | + height = max(rows) - min(rows) + 1 |
| 115 | + |
| 116 | + window_rows = [row - row_off for row in rows] |
| 117 | + window_cols = [col - col_off for col in cols] |
| 118 | + |
| 119 | + return Window(col_off, row_off, width, height), window_rows, window_cols |
| 120 | + |
| 121 | + |
| 122 | +if __name__ == "__main__": |
| 123 | + main() |
0 commit comments