Commit ae00f0d4 authored by Sarah Lohr's avatar Sarah Lohr
Browse files

Initial commit

parents
# Default ignored files
/shelf/
/workspace.xml
<?xml version="1.0" encoding="UTF-8"?>
<module type="PYTHON_MODULE" version="4">
<component name="NewModuleRootManager">
<content url="file://$MODULE_DIR$" />
<orderEntry type="inheritedJdk" />
<orderEntry type="sourceFolder" forTests="false" />
</component>
<component name="TestRunnerService">
<option name="PROJECT_TEST_RUNNER" value="pytest" />
</component>
</module>
\ No newline at end of file
<component name="InspectionProjectProfileManager">
<settings>
<option name="USE_PROJECT_PROFILE" value="false" />
<version value="1.0" />
</settings>
</component>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectRootManager" version="2" project-jdk-name="Python 3.7" project-jdk-type="Python SDK" />
</project>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectModuleManager">
<modules>
<module fileurl="file://$PROJECT_DIR$/.idea/OSM_ATKIS.iml" filepath="$PROJECT_DIR$/.idea/OSM_ATKIS.iml" />
</modules>
</component>
</project>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="VcsDirectoryMappings">
<mapping directory="$PROJECT_DIR$/.." vcs="Git" />
</component>
</project>
\ No newline at end of file
import json
import ogr
from functions.functions import split_by_fclass, diss_rel_fclasses, get_files_by_ext, zonal_ndvi_stats
import geopandas as gpd
from modules.calculate_acc_NDVI import ndvi_per_cell
import os
from modules.blue_share import calc_water_per_cell
from modules.landuse_share import calc_lu_per_cell
def main_osmatkis():
# Input parameters
input_path = "/home/sarah/PycharmProjects/blue_green_landuse/"
config_file = "./config/config.config"
roi = "study_area"
driver = ogr.GetDriverByName('ESRI Shapefile')
with open(config_file, 'r') as src:
config = json.load(src)
ds_grid = driver.Open(input_path + config[roi]["grid"], 0) # open the grid file
grid = ds_grid.GetLayer()
ds = driver.Open(input_path + config[roi]["admin_study"], 0) # open the administrative area
admin_area = ds.GetLayer()
index_file = gpd.read_file(input_path + config[roi]["grid"]) # open grid file
# MODULES TO BE EXECUTED
do_split_fclasses = False
do_dissolve = False
do_calc_acc_ndvi = False
do_spatial_ndvi_stats = True
do_calc_lu_percentage = False
do_calc_water_share = False
# Split OSM classes
if do_split_fclasses:
split_by_fclass(config[roi]["features"]["landuse"],config[roi]["outpath"], "fclass")
split_by_fclass(config[roi]["features"]["water"],config[roi]["outpath"], "fclass")
# dissolve OSM land use classes
if do_dissolve:
if not os.path.exists(config[roi]["outpath"]+ 'dissolved/'):
os.makedirs(config[roi]["outpath"] + 'dissolved/')
file_paths, names, col_names = get_files_by_ext("shp", config[roi]["outpath"])
diss_rel_fclasses(file_paths, names, config[roi]["outpath"] + 'dissolved/', "id", "fclass")
# Accumulate NDVI values per grid cell
if do_calc_acc_ndvi:
index_file = gpd.read_file(input_path + config[roi]["grid"])
ndvi_per_cell(input_path, index_file, config, roi)
# Calculate spatial statistics of NDVI values per land use class for OSM and ATKIS
if do_spatial_ndvi_stats:
# OSM
dict_osm = zonal_ndvi_stats(input_path, input_path + config[roi]["outpath"] + "dissolved/",
input_path + config[roi]['ndvi'], config, roi)
for fclass in dict_osm.keys():
print(fclass,dict_osm[fclass][0]['percentile_10'], dict_osm[fclass][0]['percentile_90'],dict_osm[fclass][0]['median'])
# ATKIS
dict_atkis =zonal_ndvi_stats(input_path, input_path + config[roi]["atkis_fixed"],
input_path + config[roi]['ndvi'], config, roi)
for fclass in dict_atkis.keys():
print(fclass,dict_atkis[fclass][0]['percentile_10'], dict_atkis[fclass][0]['percentile_90'],dict_atkis[fclass][0]['median'])
# Calculate the percentage of each land use class for each grid cell
if do_calc_lu_percentage:
out_path_diss = input_path + config[roi]["outpath"] + "dissolved/"
file_path_list_diss, names, col_names = get_files_by_ext("shp", out_path_diss)
nameslist = []
for i in names:
nameslist.append(i[:-4] + "i")
calc_lu_per_cell(out_path_diss, grid, admin_area, file_path_list_diss, nameslist, driver, index_file)
ds = None
print("-------------------------------------------------------------------------------------------------------")
# Calculate the percentage of water for each grid cell (water + river + wetland + reservoir)
if do_calc_water_share:
file_path = config[roi]["features"]["water"]
output_path = config[roi]["outpath"] + "blue_idx.shp"
calc_water_per_cell(output_path, file_path, grid, index_file)
print("-------------------------------------------------------------------------------------------------------")
if __name__ == "__main__":
main_osmatkis()
\ No newline at end of file
{"study_area": {
"name": "Whole study area",
"features":{
"landuse": "input/study_area/landuse_utm_region.shp",
"buildings":"input/study_area/buildings_utm_region.shp",
"water": "input/study_area/water_utm_region.shp",
"new_ids": "input/study_area/dissolved_new_groups.shp"},
"grid": "input/study_area/grid_region.shp",
"admin_study": "input/study_area/study_area_utm.shp",
"ndvi": "input/mannheim-heidelberg_2019-01-01_2019-12-31.NDVI.tif",
"outpath_new": "output/new/",
"outpath": "output/study_area/",
"ndvi_outpath": "output/study_area/acc_ndvi_studyarea.shp",
"ndvi_to_points": "output/study_area/ndvi_points_studyarea.asc",
"plots": "output/study_area/plots/",
"atkis_fixed":"input/study_area/atkis_fixed.shp"}
}
\ No newline at end of file
import geopandas as gpd
#import functions.functions as utils
import fiona
from shapely.geometry import Polygon
from shapely.geometry import mapping
from shapely.geometry import MultiLineString
from shapely.geometry import LineString
from fiona.crs import from_epsg
import os
import gdal
import os
import rasterio
import geopandas as gpd
from rasterstats import zonal_stats
import os
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import itertools
import geopandas as gpd
def zonal_ndvi_stats(input_path, input_files, infile_raster, config, roi):
ndvi = rasterio.open(infile_raster)
array = ndvi.read(1)
affine = ndvi.transform
#osm_output = "Projects/NDVI_Quality/zonalstats_osm.png"
if os.path.isdir(input_files): # if it is a directory --> OSM
print("Input is a directory.")
vector_infiles, names, colnames = get_files_by_ext("shp", input_files)
nameslist = []
for i in names:
nameslist.append(i[:-5])
dict = {}
for infile_vector, name in zip(vector_infiles, nameslist):
vector = gpd.read_file(infile_vector)
file = vector.to_crs(crs=ndvi.crs.data)
dict[name] = zonal_stats(file, array, affine=affine,
stats=['min', 'max', 'mean', 'median', 'range', 'std', "percentile_25",
"percentile_75","percentile_10","percentile_90"])
data = "OSM"
order = print_lists_def_order(dict)
plot_stats(order, input_path + config[roi]["plots"], data, dict)
if os.path.isfile(input_files): # if it is a file --> ATKIS
print("Input is a file.")
vector = gpd.read_file(input_files)
nameslist = []
for name in vector["type"]:
nameslist.append(name)
dict = {}
file = vector.to_crs(crs=ndvi.crs.data)
for row, name in zip(file.iterrows(), nameslist):
dict[name] = zonal_stats(row[1][11], array, affine=affine, stats=['min', 'max', 'mean', 'median', 'range',
"percentile_25", "percentile_75",
'std', "percentile_10", "percentile_90"])
data = "ATKIS"
order = print_lists_def_order(dict)
plot_stats(order, input_path + config[roi]["plots"], data, dict)
return dict
def print_lists_def_order(dict):
mean_list = []
for fclass in dict.keys():
mean_list.append((fclass, dict[fclass][0]['mean']))
mean_list = sorted(mean_list, key=lambda x: x[1], reverse=True)
print("NDVI means per land use class:", mean_list)
range_list = []
for fclass in dict.keys():
range_list.append((fclass, dict[fclass][0]['range']))
range_list = sorted(range_list, key=lambda x: x[1], reverse=True)
print("NDVI ranges per land use class:", range_list)
order = []
for pair in mean_list:
order.append(pair[0])
return order
def plot_stats(nameslist, outpath, name, dict):
for fclass in dict.keys():
dict[fclass][0]['stats'] = [{
"label": fclass, # not required
# "mean": dict[fclass][0]['mean'], # not required
"med": dict[fclass][0]['median'],
"q1": dict[fclass][0]['percentile_25'],
"q3": dict[fclass][0]['percentile_75'],
# "cilo": 5.3 # not required
# "cihi": 5.7 # not required
"whislo": dict[fclass][0]['min'], # required
"whishi": dict[fclass][0]['max'], # required
"fliers": [] # required if showfliers=True
}]
inputdata = [range(0, 4), range(0, 6)]
result = list(itertools.product(*inputdata))
fs = 10 # fontsize
fig, axes = plt.subplots(4, 6, figsize=(10, 10), constrained_layout=True, sharey=True)
# for position,fclass in zip(result, dict.keys()):
for position, fclass in zip(result, nameslist):
axes[position[0], position[1]].bxp(dict[fclass][0]['stats'], widths=0.5, patch_artist=True)
axes[position[0], position[1]].set_title(fclass, fontsize=fs)
fig.suptitle('max. NDVI of '+ name + ' landuse classes', fontsize=14)
plt.savefig(outpath + name + "_" + "zonalstats.png")
# plt.show()
return dict
def get_files_by_ext(fileext, outpath):
filepathlist = []
colnames = []
names = []
for item in os.listdir(outpath):
if os.path.isfile(outpath+item):
item_ext = item.split(".")[-1]
if item_ext == fileext:
filepathlist.append(outpath+item)
item = item[:-4]
print(item)
colnames.append(item[:8] + '_i')
names.append(item)
return filepathlist, names, colnames
def diss_rel_fclasses(file_paths,names,outpath_diss, id_col, attr_col):
print("Start dissolving feature classes...")
for file_path,name in zip(file_paths, names):
f = gpd.read_file(file_path)
# dissolve the state boundary by region
f_diss = f[[id_col, attr_col, 'geometry']].dissolve(by='fclass')
filename = outpath_diss + name + '_diss.shp'
f_diss.to_file(driver='ESRI Shapefile', filename=filename)
print("Files have been created and saved in: ", outpath_diss)
return
def write_shp_selection(outfile, file, feature, class_col, id_col):
id_list = []
geom_list = []
fclass_list = []
for index, row in file.iterrows():
if row[class_col] == feature:
# Update the value in 'area' column with area information at index
geom_list.append(row['geometry'])
id_list.append(row[id_col])
fclass_list.append(row[class_col])
newlist = []
for geom in geom_list:
newlist.append(type(geom))
polygon_schema = {'geometry': 'Polygon', 'properties': {'id': 'int', 'fclass': 'str'}}
with fiona.open(outfile, 'w', crs=from_epsg(32632), driver='ESRI Shapefile', schema=polygon_schema) as featureclass:
for geom, id, fclass in zip(geom_list, id_list, fclass_list):
featureclass.write({'geometry': mapping(geom), 'properties': {'id': id, 'fclass': fclass}})
return
def split_by_fclass(filepath, outpath, split_col):
print("Start splitting shapefiles by feature classes...")
file = gpd.read_file(filepath)
class_list = file[split_col].unique()
for feature in class_list:
outfile = outpath + str(feature) + str('.shp')
print(outfile)
write_shp_selection(outfile, file, feature, split_col, "osm_id")
print('Files have been created and saved in: ', outpath)
\ No newline at end of file
UTF-8
\ No newline at end of file
PROJCS["ETRS89_UTM_zone_32N",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]
\ No newline at end of file
PROJCS["ETRS89 / UTM zone 32N",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","25832"]]
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment