diff --git a/README.md b/README.md index ee9095381a2fed52fa6e2df5e5a3ed9c04ee869e..c9e68a671688932e544d4970a04716f7e18af0e4 100644 --- a/README.md +++ b/README.md @@ -11,7 +11,7 @@ This is the MPR tool of the > Permoserstr. 15
> 04318 Leipzig, Germany -- The current release is **[MPR v1.0.3][1]**. +- The current release is **[MPR v1.0.4][1]**. - The latest MPR release notes can be found in the file [RELEASES][3] or [online][4]. - General information can be found on the [MPR website](https://www.ufz.de/index.php?en=40126). - The MPR comes with a [LICENSE][6] agreement, this includes also the GNU Lesser General Public License. @@ -68,7 +68,7 @@ Please cite the MPR tool ([Schweppe et al. (2022)](https://doi.org/10.5194/gmd-1 You can reference the code by its [Zenodo ID](http://doi.org/10.5281/zenodo.4650513). -[1]: https://git.ufz.de/chs/mpr/tree/1.0.3 +[1]: https://git.ufz.de/chs/mpr/tree/1.0.4 [3]: doc/src/07_RELEASES.md [4]: https://git.ufz.de/chs/mpr/tags/ [5]: https://chs.pages.ufz.de/MPR diff --git a/doc/src/03_installation.md b/doc/src/03_installation.md index e41652175a16aaff8b52b410b6891c4f317d5918..4c45c885fc78d06c422b3c93081a0a3ef43550a2 100644 --- a/doc/src/03_installation.md +++ b/doc/src/03_installation.md @@ -126,7 +126,7 @@ cd mpr/ If you then want to compile a specific version (different from the latest development version), you can check that out with e.g.: ```bash -git checkout v1.0.3 +git checkout v1.0.4 ``` Afterwards you can continue with the compilation. diff --git a/doc/src/07_RELEASES.md b/doc/src/07_RELEASES.md index 93f3a2d7c9914016649a07b45970c487941245d8..aa5dfa3b7349a44d3401af85a93952f42ab7a577 100644 --- a/doc/src/07_RELEASES.md +++ b/doc/src/07_RELEASES.md @@ -1,5 +1,9 @@ # MPR Release Notes +## MPR v1.0.4 (September 2022) + +- fixed bug in upscaling to irregular grids, see (MR)[https://git.ufz.de/chs/MPR/-/merge_requests/93] + ## MPR v1.0.3 (September 2022) - updated documentation diff --git a/src/mo_mpr_data_array_upscale.F90 b/src/mo_mpr_data_array_upscale.F90 index 66a524004df046bd2dc1e40559fb8bfc0f4b11e2..41cfff762248b4010cc65b7c8dca66b3e9cc425e 100644 --- a/src/mo_mpr_data_array_upscale.F90 +++ b/src/mo_mpr_data_array_upscale.F90 @@ -1457,7 +1457,6 @@ contains subroutine compute_weights_poly(self, sourceCoord, targetCoord) !< compute the weights for remapping of two polygon coordinates (unstructured grids) !< the algorithm currently uses a simple ray casting algorithm - !< it does not work for polygons crossing the +180 degrees_east meridian class(CoordUpscaler), intent(inout) :: self type(Coordinate), intent(in) :: sourceCoord type(Coordinate), intent(in) :: targetCoord @@ -1516,7 +1515,7 @@ contains allocate(indices(sourceCoord%count)) indices = 0_i8 lastiPolygon = 1_i8 - orientation = 1_i4 + orientation = -1_i4 allocate(self%subcells(targetCoord%count)) self%subcells = 0_i8 @@ -1528,17 +1527,19 @@ contains .or. (sourceCoord%centersCoord1(iSubcell) > boundingBox(2)) & .or. (sourceCoord%centersCoord2(iSubcell) < boundingBox(3)) & .or. (sourceCoord%centersCoord2(iSubcell) > boundingBox(4))) then - log_trace(*) 'compute_weights_poly: point ', point, ' not in boundingBox: ', boundingBox + log_trace(*) 'compute_weights_poly: point ', & + [sourceCoord%centersCoord1(iSubcell), sourceCoord%centersCoord2(iSubcell)], & + ' not in boundingBox: ', boundingBox cycle end if end if - point = [sourceCoord%centersCoord1(iSubcell), sourceCoord%centersCoord2(iSubcell)] do iPolygon=1, targetCoord%count + point = [sourceCoord%centersCoord1(iSubcell), sourceCoord%centersCoord2(iSubcell)] ! check the index of last polygon first currentPolygon = mod(lastiPolygon + iPolygon - 2, targetCoord%count) + 1_i8 ! do not recalculate if same polygon is used - if (iSubcell == 1_i8 .or. currentPolygon /= lastiPolygon .or. orientation == 1_i4) then + if (iSubcell == 1_i8 .or. currentPolygon /= lastiPolygon) then nNodes = targetCoord%nodes(currentPolygon) ! select correct polygonNodes polygonNodes(1:nNodes,1) = targetCoord%cornersCoord1(1:nNodes, currentPolygon) @@ -1554,9 +1555,11 @@ contains else if (orientation == 1_i4) then ! the polygon is clockwise (this happens, if polygon crosses 180 meridian), shift coordinates polygonNodes(1:nNodes,1) = mod_shift(polygonNodes(1:nNodes,1)) - point(1) = mod_shift(point(1)) end if end if + if (orientation == 1_i4) then + point(1) = mod_shift(point(1)) + end if call inpoly(point, polygonNodes(1:nNodes,:), result) if (result > 0_i4) then log_trace(*) 'compute_weights_poly: subcell ', iSubcell, ' in target polygon ', currentPolygon diff --git a/src_python/pre_proc/convert_shp_to_SCRIP_nc.py b/src_python/pre_proc/convert_shp_to_SCRIP_nc.py index c6837ff5ebf5770040eb6148e2529251fb97ea5a..3eeb9c6aee398bfe6f5cf4e0902d3d49ac8cffc6 100755 --- a/src_python/pre_proc/convert_shp_to_SCRIP_nc.py +++ b/src_python/pre_proc/convert_shp_to_SCRIP_nc.py @@ -319,14 +319,14 @@ if __name__ == '__main__': for i_polygon in tqdm.trange(len(ds['grid_size'])): n = n_corners - (ds['grid_corner_lon'].isel(grid_size=i_polygon).values == -9999).sum() gdf.loc[i_polygon, 'geometry'] = Polygon(zip( - ds['grid_corner_lon'].isel(grid_size=i_polygon, grid_corners=slice(n)), - ds['grid_corner_lat'].isel(grid_size=i_polygon, grid_corners=slice(n)))) + ds['grid_corner_lon'].isel(grid_size=i_polygon, grid_corners=slice(n)).values, + ds['grid_corner_lat'].isel(grid_size=i_polygon, grid_corners=slice(n)).values)) else: # loop over polygons and set geometry item - for i_polygon in range(len(ds['grid_size'])): + for i_polygon in tqdm.range(len(ds['grid_size'])): gdf.loc[i_polygon, 'geometry'] = Polygon(zip( - ds['grid_corner_lon'].isel(grid_size=i_polygon), - ds['grid_corner_lat'].isel(grid_size=i_polygon))) + ds['grid_corner_lon'].isel(grid_size=i_polygon).values, + ds['grid_corner_lat'].isel(grid_size=i_polygon).values)) # set WGS84 by default gdf.crs = {'init': 'epsg:4326'} diff --git a/version.txt b/version.txt index e4c0d46e55ffb2237c9e900aa77172886f6c8aa5..a6a3a43c3a047ed2e5fc158b685fcbf65dbc5426 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -1.0.3 \ No newline at end of file +1.0.4 \ No newline at end of file