Skip to content

Allow calculate_k to use an optional buffer #67

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
patricoferris merged 1 commit into from
Oct 26, 2023
Merged

Allow calculate_k to use an optional buffer #67

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
151 changes: 86 additions & 65 deletions methods/matching/calculate_k.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,17 +1,18 @@
import argparse
import glob
import os
import tempfile
from collections import namedtuple
from itertools import product
from typing import List
from typing import List, Optional

import pandas as pd
from geopandas import gpd # type: ignore
from yirgacheffe.layers import TiledGroupLayer, RasterLayer, VectorLayer # type: ignore
from yirgacheffe.window import PixelScale # type: ignore

from methods.common import LandUseClass
from methods.common.geometry import area_for_geometry
from methods.common.geometry import area_for_geometry, expand_boundaries
from methods.common.luc import luc_range

HECTARE_WIDTH_IN_METERS = 100
Expand Down Expand Up @@ -117,83 +118,95 @@ def calculate_k(
slope_directory_path: str,
access_directory_path: str,
countries_raster_filename: str,
buffer: Optional[int],
result_dataframe_filename: str,
) -> None:

project = gpd.read_file(project_boundary_filename)

# Assumption: even if a buffer is specified, we use the unbuffered project
# boundary as the project area to work out pixel_skip
project_area_in_metres_squared = area_for_geometry(project)
project_area_in_hectares = project_area_in_metres_squared / 10_000
pixel_skip = PIXEL_SKIP_LARGE_PROJECT if (project_area_in_hectares > 250_000) else PIXEL_SKIP_SMALL_PROJECT

# everything is done at JRC resolution, so load a sample file from there first to get the ideal pixel scale
example_jrc_filename = glob.glob("*.tif", root_dir=jrc_directory_path)[0]
example_jrc_layer = RasterLayer.layer_from_file(os.path.join(jrc_directory_path, example_jrc_filename))
with tempfile.TemporaryDirectory() as tmpdir:
if buffer is not None:
expanded_project = expand_boundaries(project, buffer)
expanded_project_filename = os.path.join(tmpdir, "expanded-boundaries.geojson")
expanded_project.to_file(expanded_project_filename, driver="GeoJSON")
project_boundary_to_use = expanded_project_filename
else:
project_boundary_to_use = project_boundary_filename

project_collection = build_layer_collection(
example_jrc_layer.pixel_scale,
example_jrc_layer.projection,
list(luc_range(start_year, evaluation_year)),
[start_year, start_year - 5, start_year - 10],
project_boundary_filename,
jrc_directory_path,
cpc_directory_path,
ecoregions_directory_path,
elevation_directory_path,
slope_directory_path,
access_directory_path,
countries_raster_filename,
)
# everything is done at JRC resolution, so load a sample file from there first to get the ideal pixel scale
example_jrc_filename = glob.glob("*.tif", root_dir=jrc_directory_path)[0]
example_jrc_layer = RasterLayer.layer_from_file(os.path.join(jrc_directory_path, example_jrc_filename))

results = []
project_collection = build_layer_collection(
example_jrc_layer.pixel_scale,
example_jrc_layer.projection,
list(luc_range(start_year, evaluation_year)),
[start_year, start_year - 5, start_year - 10],
project_boundary_to_use,
jrc_directory_path,
cpc_directory_path,
ecoregions_directory_path,
elevation_directory_path,
slope_directory_path,
access_directory_path,
countries_raster_filename,
)

project_width = project_collection.boundary.window.xsize
for yoffset in range(0, project_collection.boundary.window.ysize, pixel_skip):
row_boundary = project_collection.boundary.read_array(0, yoffset, project_width, 1)
row_elevation = project_collection.elevation.read_array(0, yoffset, project_width, 1)
row_ecoregion = project_collection.ecoregions.read_array(0, yoffset, project_width, 1)
row_slope = project_collection.slope.read_array(0, yoffset, project_width, 1)
row_access = project_collection.access.read_array(0, yoffset, project_width, 1)
row_countries = project_collection.countries.read_array(0, yoffset, project_width, 1)
row_luc = [
luc.read_array(0, yoffset, project_width, 1) for luc in project_collection.lucs
]
# For CPC, which is at a different pixel_scale, we need to do a little math
coord = project_collection.boundary.latlng_for_pixel(0, yoffset)
_, cpc_yoffset = project_collection.cpcs[0].pixel_for_latlng(*coord)
row_cpc = [
cpc.read_array(0, cpc_yoffset, project_collection.cpcs[0].window.xsize, 1)
for cpc in project_collection.cpcs
]
results = []

for xoffset in range(0, project_width, pixel_skip):
if not row_boundary[0][xoffset]:
continue
lucs = [x[0][xoffset] for x in row_luc]
project_width = project_collection.boundary.window.xsize
for yoffset in range(0, project_collection.boundary.window.ysize, pixel_skip):
row_boundary = project_collection.boundary.read_array(0, yoffset, project_width, 1)
row_elevation = project_collection.elevation.read_array(0, yoffset, project_width, 1)
row_ecoregion = project_collection.ecoregions.read_array(0, yoffset, project_width, 1)
row_slope = project_collection.slope.read_array(0, yoffset, project_width, 1)
row_access = project_collection.access.read_array(0, yoffset, project_width, 1)
row_countries = project_collection.countries.read_array(0, yoffset, project_width, 1)
row_luc = [
luc.read_array(0, yoffset, project_width, 1) for luc in project_collection.lucs
]
# For CPC, which is at a different pixel_scale, we need to do a little math
coord = project_collection.boundary.latlng_for_pixel(0, yoffset)
_, cpc_yoffset = project_collection.cpcs[0].pixel_for_latlng(*coord)
row_cpc = [
cpc.read_array(0, cpc_yoffset, project_collection.cpcs[0].window.xsize, 1)
for cpc in project_collection.cpcs
]

coord = project_collection.boundary.latlng_for_pixel(xoffset, yoffset)
cpc_xoffset, _ = project_collection.cpcs[0].pixel_for_latlng(*coord)
cpcs = [x[0][cpc_xoffset] for x in row_cpc]
for xoffset in range(0, project_width, pixel_skip):
if not row_boundary[0][xoffset]:
continue
lucs = [x[0][xoffset] for x in row_luc]

results.append([
xoffset,
yoffset,
coord[0],
coord[1],
row_elevation[0][xoffset],
row_slope[0][xoffset],
row_ecoregion[0][xoffset],
row_access[0][xoffset],
row_countries[0][xoffset],
] + lucs + cpcs)
coord = project_collection.boundary.latlng_for_pixel(xoffset, yoffset)
cpc_xoffset, _ = project_collection.cpcs[0].pixel_for_latlng(*coord)
cpcs = [x[0][cpc_xoffset] for x in row_cpc]

luc_columns = [f'luc_{year}' for year in luc_range(start_year, evaluation_year)]
cpc_columns = ['cpc0_u', 'cpc0_d', 'cpc5_u', 'cpc5_d', 'cpc10_u', 'cpc10_d']
output = pd.DataFrame(
results,
columns=['x', 'y', 'lat', 'lng', 'elevation', 'slope', 'ecoregion', 'access', 'country'] \
+ luc_columns + cpc_columns
)
output.to_parquet(result_dataframe_filename)
results.append([
xoffset,
yoffset,
coord[0],
coord[1],
row_elevation[0][xoffset],
row_slope[0][xoffset],
row_ecoregion[0][xoffset],
row_access[0][xoffset],
row_countries[0][xoffset],
] + lucs + cpcs)

luc_columns = [f'luc_{year}' for year in luc_range(start_year, evaluation_year)]
cpc_columns = ['cpc0_u', 'cpc0_d', 'cpc5_u', 'cpc5_d', 'cpc10_u', 'cpc10_d']
output = pd.DataFrame(
results,
columns=['x', 'y', 'lat', 'lng', 'elevation', 'slope', 'ecoregion', 'access', 'country'] \
+ luc_columns + cpc_columns
)
output.to_parquet(result_dataframe_filename)

def main():
parser = argparse.ArgumentParser(description="Calculates sample pixels in project, aka set K")
Expand Down Expand Up @@ -274,6 +287,13 @@ def main():
dest="output_filename",
help="Destination parquet file for results."
)
parser.add_argument(
"--buffer",
type=int,
required=False,
dest="buffer",
help="The size, in metres, of the buffer to apply to the edge of the project boundary."
)
args = parser.parse_args()

calculate_k(
Expand All @@ -287,6 +307,7 @@ def main():
args.slope_directory_path,
args.access_directory_path,
args.countries_raster_filename,
args.buffer,
args.output_filename
)

Expand Down