From d87863682c0ef9725185d7b126dd56c64ea75ab0 Mon Sep 17 00:00:00 2001 From: Lukas Date: Mon, 17 Jul 2023 13:52:16 +0200 Subject: [PATCH] Update diag_shapeselect.py to work with shapely v2 (#3283) --- .../shapeselect/diag_shapeselect.py | 61 ++++++++++++------- esmvaltool/recipes/recipe_shapeselect.yml | 2 +- 2 files changed, 39 insertions(+), 24 deletions(-) diff --git a/esmvaltool/diag_scripts/shapeselect/diag_shapeselect.py b/esmvaltool/diag_scripts/shapeselect/diag_shapeselect.py index 403bc3e630..94ab3b14ec 100644 --- a/esmvaltool/diag_scripts/shapeselect/diag_shapeselect.py +++ b/esmvaltool/diag_scripts/shapeselect/diag_shapeselect.py @@ -11,8 +11,11 @@ from shapely.geometry import MultiPoint, shape from shapely.ops import nearest_points -from esmvaltool.diag_scripts.shared import (run_diagnostic, ProvenanceLogger, - get_diagnostic_filename) +from esmvaltool.diag_scripts.shared import ( + ProvenanceLogger, + get_diagnostic_filename, + run_diagnostic, +) logger = logging.getLogger(os.path.basename(__file__)) @@ -48,16 +51,22 @@ def main(cfg): xname = name + '_table' writexls(cfg, filename, ncts, nclon, nclat) caption = 'Selected gridpoints within shapefile.' - get_provenance_record( - cfg, xname, caption, 'xlsx', ancestor_files=[filename]) + get_provenance_record(cfg, + xname, + caption, + 'xlsx', + ancestor_files=[filename]) path = os.path.join( cfg['work_dir'], name + '.nc', ) write_netcdf(path, ncts, nclon, nclat, cube, cfg) caption = 'Selected gridpoints within shapefile.' - get_provenance_record( - cfg, name, caption, 'nc', ancestor_files=[filename]) + get_provenance_record(cfg, + name, + caption, + 'nc', + ancestor_files=[filename]) def write_keyvalue_toxlsx(worksheet, row, key, value): @@ -93,17 +102,16 @@ def writexls(cfg, filename, ncts, nclon1, nclat1): workbook = xlsxwriter.Workbook( os.path.join( cfg['work_dir'], - os.path.splitext(os.path.basename(filename))[0] + '_polygon_table' - + '.xlsx')) + os.path.splitext(os.path.basename(filename))[0] + + '_polygon_table' + '.xlsx')) worksheet = workbook.add_worksheet('Data') worksheet.write(0, 0, 'Date') worksheet.write(0, 1, 'Lon/Lat') worksheet.write_column(2, 0, wtime) for row in range(ncts.shape[1]): worksheet.write( - 1, row + 1, - str("%#.3f" % round(float(nclon1[row]), 3)) + '_' + str( - "%#.3f" % round(float(nclat1[row]), 3))) + 1, row + 1, f"{round(float(nclon1[row]), 3):.3f}_\ + {round(float(nclat1[row]), 3)}") worksheet.write_column(2, row + 1, np.around(np.squeeze(ncts[:, row]), decimals=8)) worksheet.set_column(0, row + 1, 20) @@ -135,7 +143,8 @@ def shapeselect(cfg, cube): coordpoints[i] = (coordpoints[i][0] - 360., coordpoints[i][1]) else: raise ValueError("Support for 2-d coords not implemented!") - points = MultiPoint(coordpoints) + multipoint = MultiPoint(coordpoints) + points = list(multipoint.geoms) with fiona.open(shppath) as shp: gpx = [] gpy = [] @@ -149,15 +158,16 @@ def shapeselect(cfg, cube): if wgtmet == 'mean_inside': gpx, gpy = mean_inside(gpx, gpy, points, multi, cube) if not gpx: - gpx, gpy = representative(gpx, gpy, points, multi, cube) + gpx, gpy = representative(gpx, gpy, multipoint, multi, + cube) elif wgtmet == 'representative': - gpx, gpy = representative(gpx, gpy, points, multi, cube) + gpx, gpy = representative(gpx, gpy, multipoint, multi, cube) if len(gpx) == 1: ncts[:, ishp] = np.reshape(cube.data[:, gpy, gpx], (cube.data.shape[0], )) else: ncts[:, ishp] = np.mean(cube.data[:, gpy, gpx], axis=1) - gxx, gyy = representative([], [], points, multi, cube) + gxx, gyy = representative([], [], multipoint, multi, cube) nclon[ishp] = cube.coord('longitude').points[gxx] nclat[ishp] = cube.coord('latitude').points[gyy] return ncts, nclon, nclat @@ -179,10 +189,10 @@ def mean_inside(gpx, gpy, points, multi, cube): return gpx, gpy -def representative(gpx, gpy, points, multi, cube): +def representative(gpx, gpy, multipoint, multi, cube): """Find representative point in shape.""" reprpoint = multi.representative_point() - nearest = nearest_points(reprpoint, points) + nearest = nearest_points(reprpoint, multipoint) npx = nearest[1].coords[0][0] npy = nearest[1].coords[0][1] if npx < 0: @@ -235,19 +245,24 @@ def write_netcdf(path, var, plon, plat, cube, cfg): polys.setncattr_string('standard_name', 'polygon') polys.setncattr_string('long_name', 'polygon') polys.setncattr_string('shapefile', cfg['shapefile']) - lon = ncout.createVariable( - cube.coord('longitude').var_name, 'f8', 'polygon', zlib=True) + lon = ncout.createVariable(cube.coord('longitude').var_name, + 'f8', + 'polygon', + zlib=True) lon.setncattr_string('standard_name', cube.coord('longitude').standard_name) lon.setncattr_string('long_name', cube.coord('longitude').long_name) lon.setncattr_string('units', cube.coord('longitude').units.origin) - lat = ncout.createVariable( - cube.coord('latitude').var_name, 'f8', 'polygon', zlib=True) + lat = ncout.createVariable(cube.coord('latitude').var_name, + 'f8', + 'polygon', + zlib=True) lat.setncattr_string('standard_name', cube.coord('latitude').standard_name) lat.setncattr_string('long_name', cube.coord('latitude').long_name) lat.setncattr_string('units', cube.coord('latitude').units.origin) - data = ncout.createVariable( - cube.var_name, 'f4', ('time', 'polygon'), zlib=True) + data = ncout.createVariable(cube.var_name, + 'f4', ('time', 'polygon'), + zlib=True) data.setncattr_string('standard_name', cube.standard_name) data.setncattr_string('long_name', cube.long_name) data.setncattr_string('units', cube.units.origin) diff --git a/esmvaltool/recipes/recipe_shapeselect.yml b/esmvaltool/recipes/recipe_shapeselect.yml index f37a302dff..0fb22c0d5d 100644 --- a/esmvaltool/recipes/recipe_shapeselect.yml +++ b/esmvaltool/recipes/recipe_shapeselect.yml @@ -11,7 +11,7 @@ documentation: - berg_peter maintainer: - - unmaintained + - ruhe_lukas projects: - c3s-magic