From 6099b62a84c0d4245cff47f12ad32e512b5f43be Mon Sep 17 00:00:00 2001 From: Robert Schweppe Date: Mon, 26 Sep 2022 12:42:42 +0200 Subject: [PATCH 1/9] - always set point in inner loop as it may be overwritten later by mod_shift in case orientation == 1 --- src/mo_mpr_data_array_upscale.F90 | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/mo_mpr_data_array_upscale.F90 b/src/mo_mpr_data_array_upscale.F90 index 66a5240..5b115f6 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 @@ -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) -- GitLab From afb7011a9077389ed32dbddab8c5d033b7c868cc Mon Sep 17 00:00:00 2001 From: Robert Schweppe Date: Mon, 26 Sep 2022 14:07:56 +0200 Subject: [PATCH 2/9] - always set point in inner loop as it may be overwritten later by mod_shift in case orientation == 1 --- README.md | 4 ++-- doc/src/03_installation.md | 2 +- doc/src/07_RELEASES.md | 4 ++++ version.txt | 2 +- 4 files changed, 8 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index ee90953..c9e68a6 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 e416521..4c45c88 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 93f3a2d..aa5dfa3 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/version.txt b/version.txt index e4c0d46..a6a3a43 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 -- GitLab From 2c7ba6a90c61facbfc436420edf0f3934252d422 Mon Sep 17 00:00:00 2001 From: Robert Schweppe Date: Wed, 28 Sep 2022 15:17:18 +0200 Subject: [PATCH 3/9] - set debug settings, as there are still problems --- src/mo_mpr_data_array_upscale.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/mo_mpr_data_array_upscale.F90 b/src/mo_mpr_data_array_upscale.F90 index 5b115f6..4687aec 100644 --- a/src/mo_mpr_data_array_upscale.F90 +++ b/src/mo_mpr_data_array_upscale.F90 @@ -1545,7 +1545,7 @@ contains polygonNodes(1:nNodes,1) = targetCoord%cornersCoord1(1:nNodes, currentPolygon) polygonNodes(1:nNodes,2) = targetCoord%cornersCoord2(1:nNodes, currentPolygon) orientation = orientpoly(polygonNodes(1:nNodes,:)) - log_trace(*) 'compute_weights_poly: got orientation ', orientation, ' for polygon ', currentPolygon, & + log_info(*) 'compute_weights_poly: got orientation ', orientation, ' for polygon ', currentPolygon, & ' with nodes ', polygonNodes(1:nNodes,:), ' for point ', point if (orientation == 0_i4) then @@ -1567,6 +1567,7 @@ contains exit end if end do + stop 1 lastiPolygon = currentPolygon end do maxSubcells = maxval(self%subcells) -- GitLab From fc4fcdce29bfd6f9d0a570ff2d16c1314dbe1c75 Mon Sep 17 00:00:00 2001 From: Robert Schweppe Date: Wed, 28 Sep 2022 15:21:11 +0200 Subject: [PATCH 4/9] - set debug settings, as there are still problems --- src/mo_mpr_data_array_upscale.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mo_mpr_data_array_upscale.F90 b/src/mo_mpr_data_array_upscale.F90 index 4687aec..59ac556 100644 --- a/src/mo_mpr_data_array_upscale.F90 +++ b/src/mo_mpr_data_array_upscale.F90 @@ -1558,6 +1558,8 @@ contains point(1) = mod_shift(point(1)) end if end if + print*, point + print*, polygonNodes(1:nNodes,:) call inpoly(point, polygonNodes(1:nNodes,:), result) if (result > 0_i4) then log_trace(*) 'compute_weights_poly: subcell ', iSubcell, ' in target polygon ', currentPolygon -- GitLab From bd6507f9db12781e26e48df290bd99a3b4648a36 Mon Sep 17 00:00:00 2001 From: Robert Schweppe Date: Wed, 28 Sep 2022 15:43:14 +0200 Subject: [PATCH 5/9] - set debug settings, as there are still problems --- src/mo_mpr_data_array_upscale.F90 | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/mo_mpr_data_array_upscale.F90 b/src/mo_mpr_data_array_upscale.F90 index 59ac556..9be47b5 100644 --- a/src/mo_mpr_data_array_upscale.F90 +++ b/src/mo_mpr_data_array_upscale.F90 @@ -1545,7 +1545,7 @@ contains polygonNodes(1:nNodes,1) = targetCoord%cornersCoord1(1:nNodes, currentPolygon) polygonNodes(1:nNodes,2) = targetCoord%cornersCoord2(1:nNodes, currentPolygon) orientation = orientpoly(polygonNodes(1:nNodes,:)) - log_info(*) 'compute_weights_poly: got orientation ', orientation, ' for polygon ', currentPolygon, & + log_trace(*) 'compute_weights_poly: got orientation ', orientation, ' for polygon ', currentPolygon, & ' with nodes ', polygonNodes(1:nNodes,:), ' for point ', point if (orientation == 0_i4) then @@ -1558,18 +1558,21 @@ contains point(1) = mod_shift(point(1)) end if end if - print*, point - print*, polygonNodes(1:nNodes,:) call inpoly(point, polygonNodes(1:nNodes,:), result) if (result > 0_i4) then + log_info(*) 'original point:', [sourceCoord%centersCoord1(iSubcell), sourceCoord%centersCoord2(iSubcell)] + log_info(*) 'original nodes:', [targetCoord%cornersCoord1(1:nNodes, currentPolygon), targetCoord%cornersCoord2(1:nNodes, currentPolygon)] + log_info(*) 'compute_weights_poly: got orientation ', orientation, ' for polygon ', currentPolygon, & + ' with nodes ', polygonNodes(1:nNodes,:), ' for point ', point + log_trace(*) 'compute_weights_poly: subcell ', iSubcell, ' in target polygon ', currentPolygon + stop 1 indices(iSubcell) = currentPolygon ! advance counter self%subcells(currentPolygon) = self%subcells(currentPolygon) + 1_i8 exit end if end do - stop 1 lastiPolygon = currentPolygon end do maxSubcells = maxval(self%subcells) -- GitLab From 16cf2b2479bb6660a6ca856aa6bfaaf4d7ac452a Mon Sep 17 00:00:00 2001 From: Robert Schweppe Date: Wed, 28 Sep 2022 15:57:36 +0200 Subject: [PATCH 6/9] - set debug settings, as there are still problems --- src/mo_mpr_data_array_upscale.F90 | 7 +------ src/mo_mpr_upscale_func.F90 | 6 ++++++ 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/mo_mpr_data_array_upscale.F90 b/src/mo_mpr_data_array_upscale.F90 index 9be47b5..a70f39d 100644 --- a/src/mo_mpr_data_array_upscale.F90 +++ b/src/mo_mpr_data_array_upscale.F90 @@ -1561,12 +1561,7 @@ contains call inpoly(point, polygonNodes(1:nNodes,:), result) if (result > 0_i4) then log_info(*) 'original point:', [sourceCoord%centersCoord1(iSubcell), sourceCoord%centersCoord2(iSubcell)] - log_info(*) 'original nodes:', [targetCoord%cornersCoord1(1:nNodes, currentPolygon), targetCoord%cornersCoord2(1:nNodes, currentPolygon)] - log_info(*) 'compute_weights_poly: got orientation ', orientation, ' for polygon ', currentPolygon, & - ' with nodes ', polygonNodes(1:nNodes,:), ' for point ', point - - log_trace(*) 'compute_weights_poly: subcell ', iSubcell, ' in target polygon ', currentPolygon - stop 1 + log_info(*) 'compute_weights_poly: subcell ', iSubcell, ' in target polygon ', currentPolygon indices(iSubcell) = currentPolygon ! advance counter self%subcells(currentPolygon) = self%subcells(currentPolygon) + 1_i8 diff --git a/src/mo_mpr_upscale_func.F90 b/src/mo_mpr_upscale_func.F90 index f4e60fd..eb85fcd 100644 --- a/src/mo_mpr_upscale_func.F90 +++ b/src/mo_mpr_upscale_func.F90 @@ -118,6 +118,7 @@ contains real(dp), dimension(size(weightsIn)) :: weightsAdapted if (all(.not. maskIn)) then + print*, 'all missing' ! if all values are masked missing, then the result is nan and also masked missing valueOut = nodata_dp maskOut = .false. @@ -133,9 +134,14 @@ contains else weightsAdapted = weightsIn / sum(weightsIn, mask=maskIn) end if + print*, 'sliceIn', sliceIn + print*, 'maskIn', maskIn + print*, 'size', size(sliceIn) maskOut = .true. ! the non-missing values are send to func together with their weights valueOut = func(pack(sliceIn, mask=maskIn), pack(weightsAdapted, mask=maskIn), p) + print*, 'valueOut', valueOut + stop 1 end if end subroutine wrap_weighted_upscale -- GitLab From 4acf53d59902ef3b64cf2e813c46d9e3de72a588 Mon Sep 17 00:00:00 2001 From: Robert Schweppe Date: Wed, 28 Sep 2022 16:03:15 +0200 Subject: [PATCH 7/9] - set debug settings, as there are still problems --- src/mo_mpr_data_array_upscale.F90 | 1 + src/mo_mpr_upscale_func.F90 | 6 ------ 2 files changed, 1 insertion(+), 6 deletions(-) diff --git a/src/mo_mpr_data_array_upscale.F90 b/src/mo_mpr_data_array_upscale.F90 index a70f39d..e8ed2d4 100644 --- a/src/mo_mpr_data_array_upscale.F90 +++ b/src/mo_mpr_data_array_upscale.F90 @@ -1570,6 +1570,7 @@ contains end do lastiPolygon = currentPolygon end do + stop 1 maxSubcells = maxval(self%subcells) log_debug(*) 'compute_weights_poly: setting ids and weights for Upscaler with dimensions (', & maxSubcells, 'x', targetCoord%count, ')' diff --git a/src/mo_mpr_upscale_func.F90 b/src/mo_mpr_upscale_func.F90 index eb85fcd..f4e60fd 100644 --- a/src/mo_mpr_upscale_func.F90 +++ b/src/mo_mpr_upscale_func.F90 @@ -118,7 +118,6 @@ contains real(dp), dimension(size(weightsIn)) :: weightsAdapted if (all(.not. maskIn)) then - print*, 'all missing' ! if all values are masked missing, then the result is nan and also masked missing valueOut = nodata_dp maskOut = .false. @@ -134,14 +133,9 @@ contains else weightsAdapted = weightsIn / sum(weightsIn, mask=maskIn) end if - print*, 'sliceIn', sliceIn - print*, 'maskIn', maskIn - print*, 'size', size(sliceIn) maskOut = .true. ! the non-missing values are send to func together with their weights valueOut = func(pack(sliceIn, mask=maskIn), pack(weightsAdapted, mask=maskIn), p) - print*, 'valueOut', valueOut - stop 1 end if end subroutine wrap_weighted_upscale -- GitLab From 59b3420fbb203165def8dbd2424a9386fdc2b7e9 Mon Sep 17 00:00:00 2001 From: Robert Schweppe Date: Wed, 28 Sep 2022 16:18:59 +0200 Subject: [PATCH 8/9] - finally fixed bug --- src/mo_mpr_data_array_upscale.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mo_mpr_data_array_upscale.F90 b/src/mo_mpr_data_array_upscale.F90 index e8ed2d4..a18db71 100644 --- a/src/mo_mpr_data_array_upscale.F90 +++ b/src/mo_mpr_data_array_upscale.F90 @@ -1555,13 +1555,14 @@ 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_info(*) 'original point:', [sourceCoord%centersCoord1(iSubcell), sourceCoord%centersCoord2(iSubcell)] - log_info(*) 'compute_weights_poly: subcell ', iSubcell, ' in target polygon ', currentPolygon + log_trace(*) 'compute_weights_poly: subcell ', iSubcell, ' in target polygon ', currentPolygon indices(iSubcell) = currentPolygon ! advance counter self%subcells(currentPolygon) = self%subcells(currentPolygon) + 1_i8 @@ -1570,7 +1571,6 @@ contains end do lastiPolygon = currentPolygon end do - stop 1 maxSubcells = maxval(self%subcells) log_debug(*) 'compute_weights_poly: setting ids and weights for Upscaler with dimensions (', & maxSubcells, 'x', targetCoord%count, ')' -- GitLab From 21d2f55aed64578a2079d0bbf8a736570e8434cf Mon Sep 17 00:00:00 2001 From: Robert Schweppe Date: Wed, 28 Sep 2022 16:27:45 +0200 Subject: [PATCH 9/9] - reset orientation init flag - updated convert_shp_to_SCRIP_nc.py script --- src/mo_mpr_data_array_upscale.F90 | 2 +- src_python/pre_proc/convert_shp_to_SCRIP_nc.py | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/mo_mpr_data_array_upscale.F90 b/src/mo_mpr_data_array_upscale.F90 index a18db71..41cfff7 100644 --- a/src/mo_mpr_data_array_upscale.F90 +++ b/src/mo_mpr_data_array_upscale.F90 @@ -1515,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 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 c6837ff..3eeb9c6 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'} -- GitLab