Commit 2308d41b authored by Robert Schweppe's avatar Robert Schweppe
Browse files

Merge branch 'hotfix_polygon_resize' into 'develop'

- merged develop in reactivated branch hotfix_polygon_resize

See merge request !87
parents 68aa83bb e1f9be2c
Pipeline #95864 failed with stages
in 18 minutes and 31 seconds
......@@ -807,7 +807,62 @@ contains
! read Dataset and store the variable of interest
nc = NcDataset(fileName, "r")
if (nc%hasVariable(trim(self%name))) then
if (nc%hasAttribute('title')) then
! SCRIP format
call nc%getAttribute('title', ncTitle)
self%name = trim(ncTitle)
! found the title, now read all the different field information
! start with the coordinates
call read_coordinate_from_file(nc, 'grid_size', dimSize)
call self%set_count(int(dimSize, kind=i8))
! assuming a rank of 2
! SCRIP uses rank 1 for unstructured grids, we use it differently
! call read_coordinate_from_file(nc, 'grid_rank', self%rank)
self%rank = 2_i4
self%subDims =['lon', 'lat']
call read_coordinate_from_file(nc, 'grid_corners', self%corners)
allocate(self%cornersCoord1(self%corners, self%count))
allocate(self%cornersCoord2(self%corners, self%count))
allocate(self%centersCoord1(self%count))
allocate(self%centersCoord2(self%count))
allocate(self%nodes(self%count))
call read_variable_from_file(nc, 'grid_center_lon', self%centersCoord1)
call read_variable_from_file(nc, 'grid_center_lat', self%centersCoord2)
call read_variable_from_file(nc, 'grid_corner_lon', self%cornersCoord1)
! reading also the mask to infer the number of corners/nodes per grid cell
call read_variable_from_file(nc, 'grid_corner_lat', self%cornersCoord2, ncMask)
if (.not. all(ncMask)) then
do i = 1, self%count
self%nodes(i) = self%corners - count(.not. ncMask(:, i))
end do
else
! set the number of nodes to the maximum number of corners
self%nodes = self%corners
end if
deallocate(ncMask)
! for now ignore the mask as this applied later on with the data array only
! call self%read_mask_from_file(nc, 'grid_imask', self%mask)
call read_variable_from_file(nc, 'grid_dims', dimSizes)
self%subDimSizes = int(dimSizes, kind=i8)
ncVar = nc%getVariable('grid_center_lon')
! this is mimicking the retrieval of grids.f in SCRIP -> assume unit to be similar over all variables
call ncVar%getAttribute('units', self%unit)
log_subtrace(*) 'Reading from SCRIP dimension:'
log_subtrace(*) 'grid_size ', self%count
log_subtrace(*) 'grid_rank ', self%rank
log_subtrace(*) 'grid_corners ', self%corners
log_subtrace(*) 'grid_corner_lon ', self%cornersCoord1
log_subtrace(*) 'grid_corner_lat ', self%cornersCoord2
log_subtrace(*) 'grid_center_lon ', self%centersCoord1
log_subtrace(*) 'grid_center_lat ', self%centersCoord2
log_subtrace(*) 'grid_dims ', self%subDimSizes
log_subtrace(*) 'unit ', trim(self%unit)
else if (nc%hasVariable(trim(self%name))) then
! the coordinate is actually a variable in the NcDataset
ncVar = nc%getVariable(trim(self%name))
varShape = ncVar%getShape()
......@@ -918,62 +973,6 @@ contains
if (ncVar%hasAttribute('units')) then
call ncVar%getAttribute('units', self%unit)
end if
else if (nc%hasAttribute('title')) then
! SCRIP format
call nc%getAttribute('title', ncTitle)
self%name = trim(ncTitle)
! found the title, now read all the different field information
! start with the coordinates
call read_coordinate_from_file(nc, 'grid_size', dimSize)
call self%set_count(int(dimSize, kind=i8))
! assuming a rank of 2
! SCRIP uses rank 1 for unstructured grids, we use it differently
! call read_coordinate_from_file(nc, 'grid_rank', self%rank)
self%rank = 2_i4
self%subDims =['lon', 'lat']
call read_coordinate_from_file(nc, 'grid_corners', self%corners)
allocate(self%cornersCoord1(self%corners, self%count))
allocate(self%cornersCoord2(self%corners, self%count))
allocate(self%centersCoord1(self%count))
allocate(self%centersCoord2(self%count))
allocate(self%nodes(self%count))
call read_variable_from_file(nc, 'grid_center_lon', self%centersCoord1)
call read_variable_from_file(nc, 'grid_center_lat', self%centersCoord2)
call read_variable_from_file(nc, 'grid_corner_lon', self%cornersCoord1)
! reading also the mask to infer the number of corners/nodes per grid cell
call read_variable_from_file(nc, 'grid_corner_lat', self%cornersCoord2, ncMask)
if (.not. all(ncMask)) then
do i = 1, self%count
self%nodes(i) = self%corners - count(.not. ncMask(:, i))
end do
else
! set the number of nodes to the maximum number of corners
self%nodes = self%corners
end if
deallocate(ncMask)
! for now ignore the mask as this applied later on with the data array only
! call self%read_mask_from_file(nc, 'grid_imask', self%mask)
call read_variable_from_file(nc, 'grid_dims', dimSizes)
self%subDimSizes = int(dimSizes, kind=i8)
ncVar = nc%getVariable('grid_center_lon')
! this is mimicking the retrieval of grids.f in SCRIP -> assume unit to be similar over all variables
call ncVar%getAttribute('units', self%unit)
log_subtrace(*) 'Reading from SCRIP dimension:'
log_subtrace(*) 'grid_size ', self%count
log_subtrace(*) 'grid_rank ', self%rank
log_subtrace(*) 'grid_corners ', self%corners
log_subtrace(*) 'grid_corner_lon ', self%cornersCoord1
log_subtrace(*) 'grid_corner_lat ', self%cornersCoord2
log_subtrace(*) 'grid_center_lon ', self%centersCoord1
log_subtrace(*) 'grid_center_lat ', self%centersCoord2
log_subtrace(*) 'grid_dims ', self%subDimSizes
log_subtrace(*) 'unit ', trim(self%unit)
else
log_debug(*) "Cannot properly initialize coordinate '", trim(self%name), "'."
call nc%close()
......@@ -1622,7 +1621,7 @@ contains
else if (self%is_polygon()) then
is_finalized = .true.
! check if the polygon based coordinate is initialized
if (size(self%centersCoord1) <= 2_i4 .or. size(self%centersCoord2) <= 2_i4) then
if (size(self%centersCoord1) == 0_i4 .or. size(self%centersCoord2) == 0_i4) then
is_finalized = .false.
return
end if
......
......@@ -1469,32 +1469,48 @@ contains
real(dp), dimension(:,:), allocatable :: polygonNodes
integer(i4) :: orientation
integer(i8), dimension(:), allocatable :: indices, insertIndices
logical :: useBoundingBox
log_debug(*) 'compute_weights_poly: computing weights for 2d-coordinate to polygon coordinate variables ', &
trim(sourceCoord%name), ' to ', trim(targetCoord%name)
self%mapMethod = 'Nearest destination to source'
useBoundingBox = .true.
allocate(polygonNodes(maxval(targetCoord%nodes), 2))
do iPolygon=1, targetCoord%count
nNodes = targetCoord%nodes(iPolygon)
! select correct polygonNodes
polygonNodes(1:nNodes,1) = targetCoord%cornersCoord1(1:nNodes, iPolygon)
polygonNodes(1:nNodes,2) = targetCoord%cornersCoord2(1:nNodes, iPolygon)
orientation = orientpoly(polygonNodes(1:nNodes,:))
if (orientation /= -1_i4) then
useBoundingBox = .false.
exit
end if
end do
! restrain points to be considered for search by min and max for x and y
if (all(targetCoord%nodes == targetCoord%nodes(1))) then
! all polygons have same number of corners
boundingBox(1) = minval(targetCoord%cornersCoord1(1:targetCoord%nodes(1), :)) - eps_dp
boundingBox(2) = maxval(targetCoord%cornersCoord1(1:targetCoord%nodes(1), :)) + eps_dp
boundingBox(3) = minval(targetCoord%cornersCoord2(1:targetCoord%nodes(1), :)) - eps_dp
boundingBox(4) = maxval(targetCoord%cornersCoord2(1:targetCoord%nodes(1), :)) + eps_dp
allocate(polygonNodes(targetCoord%nodes(1), 2))
log_trace(*) 'compute_weights_poly: all polygons have same number of nodes: ', targetCoord%nodes(1)
else
! polygons have different numbers of corners
boundingBox = [targetCoord%cornersCoord1(1, 1), targetCoord%cornersCoord1(1, 1), &
targetCoord%cornersCoord2(1, 1), targetCoord%cornersCoord2(1, 1)]
do iPolygon=1, targetCoord%count
boundingBox(1) = minval([boundingBox(1), minval(targetCoord%cornersCoord1(1:targetCoord%nodes(iPolygon), iPolygon)) - eps_dp])
boundingBox(2) = maxval([boundingBox(2), maxval(targetCoord%cornersCoord1(1:targetCoord%nodes(iPolygon), iPolygon)) + eps_dp])
boundingBox(3) = minval([boundingBox(3), minval(targetCoord%cornersCoord2(1:targetCoord%nodes(iPolygon), iPolygon)) - eps_dp])
boundingBox(4) = maxval([boundingBox(4), maxval(targetCoord%cornersCoord2(1:targetCoord%nodes(iPolygon), iPolygon)) + eps_dp])
end do
allocate(polygonNodes(maxval(targetCoord%nodes), 2))
log_trace(*) 'compute_weights_poly: polygons have different number of nodes: ', maxval(targetCoord%nodes)
if (useBoundingBox) then
if (all(targetCoord%nodes == targetCoord%nodes(1))) then
! all polygons have same number of corners
boundingBox(1) = minval(targetCoord%cornersCoord1(1:targetCoord%nodes(1), :)) - eps_dp
boundingBox(2) = maxval(targetCoord%cornersCoord1(1:targetCoord%nodes(1), :)) + eps_dp
boundingBox(3) = minval(targetCoord%cornersCoord2(1:targetCoord%nodes(1), :)) - eps_dp
boundingBox(4) = maxval(targetCoord%cornersCoord2(1:targetCoord%nodes(1), :)) + eps_dp
log_trace(*) 'compute_weights_poly: all polygons have same number of nodes: ', targetCoord%nodes(1)
else
! polygons have different numbers of corners
boundingBox = [targetCoord%cornersCoord1(1, 1), targetCoord%cornersCoord1(1, 1), &
targetCoord%cornersCoord2(1, 1), targetCoord%cornersCoord2(1, 1)]
do iPolygon=1, targetCoord%count
boundingBox(1) = minval([boundingBox(1), minval(targetCoord%cornersCoord1(1:targetCoord%nodes(iPolygon), iPolygon)) - eps_dp])
boundingBox(2) = maxval([boundingBox(2), maxval(targetCoord%cornersCoord1(1:targetCoord%nodes(iPolygon), iPolygon)) + eps_dp])
boundingBox(3) = minval([boundingBox(3), minval(targetCoord%cornersCoord2(1:targetCoord%nodes(iPolygon), iPolygon)) - eps_dp])
boundingBox(4) = maxval([boundingBox(4), maxval(targetCoord%cornersCoord2(1:targetCoord%nodes(iPolygon), iPolygon)) + eps_dp])
end do
log_trace(*) 'compute_weights_poly: polygons have different number of nodes: ', maxval(targetCoord%nodes)
end if
end if
allocate(indices(sourceCoord%count))
......@@ -1507,12 +1523,14 @@ contains
! loop over source cells
do iSubcell=1, sourceCoord%count
! check if within bounds
if ((sourceCoord%centersCoord1(iSubcell) < boundingBox(1)) &
.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
cycle
if (useBoundingBox) then
if ((sourceCoord%centersCoord1(iSubcell) < boundingBox(1)) &
.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
cycle
end if
end if
point = [sourceCoord%centersCoord1(iSubcell), sourceCoord%centersCoord2(iSubcell)]
......@@ -1526,7 +1544,8 @@ 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_trace(*) 'compute_weights_poly: got orientation ', orientation, ' for polygon ', currentPolygon, &
' with nodes ', polygonNodes(1:nNodes,:), ' for point ', point
if (orientation == 0_i4) then
! the polygon covers a pole, so we need to modify the coordinates
......@@ -1538,13 +1557,13 @@ contains
point(1) = mod_shift(point(1))
end if
end if
lastiPolygon = currentPolygon
call inpoly(point, polygonNodes(1:nNodes,:), result)
if (result > 0_i4) then
log_trace(*) 'compute_weights_poly: subcell ', iSubcell, ' in target polygon ', currentPolygon
indices(iSubcell) = currentPolygon
! advance counter
self%subcells(currentPolygon) = self%subcells(currentPolygon) + 1_i8
lastiPolygon = currentPolygon
exit
end if
end do
......
......@@ -274,14 +274,14 @@ if __name__ == '__main__':
gdf = handle_latlon_format(ds)
else:
# do we have special MPI-ICON names?
if all(mpi_icon_format):
if mpi_icon_format:
# select and rename
ds = ds[list(RENAME_VAR_DICT.keys()) + SELECT_VARS].rename({**RENAME_VAR_DICT, **RENAME_DIM_DICT})
# set the grid_imask property based on sea_land_mask or to default 1
if 'cell_sea_land_mask' in ds:
ds['grid_imask'] = xr.where(ds['cell_sea_land_mask'] > 1, 1, 0)
else:
ds['grid_imask'] = (('grid_size',), np.ones_like(ds['grid_center_lon'].values, dtype=int))
if 'cell_sea_land_mask' in ds:
ds['grid_imask'] = xr.where(ds['cell_sea_land_mask'] > 1, 1, 0)
else:
ds['grid_imask'] = (('grid_size',), np.ones_like(ds['grid_center_lon'].values, dtype=int))
# select all data_vars that have dimension grid_size
shp_vars = [data_var for data_var in ds.data_vars if
'grid_size' in ds[data_var].dims and data_var not in EXCLUDE_DIMS]
......
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