|
◆ grid_transform_grid_vol7d_init()
subroutine grid_transform_class::grid_transform_grid_vol7d_init |
( |
type(grid_transform), intent(out) |
this, |
|
|
type(transform_def), intent(in) |
trans, |
|
|
type(griddim_def), intent(in) |
in, |
|
|
type(vol7d), intent(inout) |
v7d_out, |
|
|
real, dimension(:,:), intent(in), optional |
maskgrid, |
|
|
real, dimension(:), intent(in), optional |
maskbounds, |
|
|
procedure(basic_find_index), optional, pointer |
find_index, |
|
|
character(len=*), intent(in), optional |
categoryappend |
|
) |
| |
|
private |
Constructor for a grid_transform object, defining a particular grid-to-sparse points transformation.
It defines an object describing a transformation from a rectangular grid to a set of sparse points; the abstract type of transformation is described in the transformation object trans (type transform_def) which must have been properly initialised. The additional information required here is the description of the input grid in (type griddim_def), and, if required by the transformation type, the information about the target sparse points over which the transformation should take place:
- for 'inter' transformation, this is provided in the form of a vol7d object (v7d_out argument, input), which must have been initialized with the coordinates of desired sparse points
- for 'polyinter' transformation, no target point information has to be provided in input (it is calculated on the basis of input grid and trans object), and the coordinates of the target points (polygons' centroids) are returned in output in v7d_out argument
- for 'maskinter' transformation, this is a two dimensional real field (maskgrid argument), which, together with the maskbounds argument (optional with default), divides the input grid in a number of subareas according to the values of maskinter, and, in this case, v7d_out is an output argument which returns the coordinates of the target points (subareas' centroids)
- for 'metamorphosis' transformation, no target point information has to be provided in input (it is calculated on the basis of input grid and trans object), except for 'mask' subtype, for which the same information as for 'maskinter' transformation has to be provided; in all the cases, as for 'polyinter', the information about target points is returned in output in v7d_out argument.
The generated grid_transform object is specific to the grid and sparse point list provided or computed. The function c_e can be used in order to check whether the object has been successfully initialised, if the result is .FALSE., it should not be used further on. - Parametri
-
[out] | this | grid_transformation object |
[in] | trans | transformation object |
[in] | in | griddim object to transform |
[in,out] | v7d_out | vol7d object with the coordinates of the sparse points to be used as transformation target (input or output depending on type of transformation) |
[in] | maskgrid | 2D field to be used for defining subareas according to its values, it must have the same shape as the field to be interpolated (for transformation type 'maskinter' and 'metamorphosis:mask') |
[in] | maskbounds | array of boundary values for defining subareas from the values of maskgrid, the number of subareas is SIZE(maskbounds) - 1, if not provided a default based on extreme values of maskgrid is used |
[in] | categoryappend | append this suffix to log4fortran namespace category |
Definizione alla linea 2016 del file grid_transform_class.F90.
2018 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
2019 this%inter_index_x(ix,iy) > xnmax .OR. &
2020 this%inter_index_y(ix,iy) < ynmin .OR. &
2021 this%inter_index_y(ix,iy) > ynmax) THEN
2022 this%inter_index_x(ix,iy) = imiss
2023 this%inter_index_y(ix,iy) = imiss
2029 CALL l4f_category_log(this%category, l4f_debug, &
2030 'stencilinter: stencil size '//t2c(n*n))
2031 CALL l4f_category_log(this%category, l4f_debug, &
2032 'stencilinter: stencil points '//t2c(count(this%stencil)))
2040 ELSE IF (this%trans%trans_type == 'maskinter') THEN
2042 IF (.NOT. PRESENT(maskgrid)) THEN
2043 CALL l4f_category_log(this%category,l4f_error, &
2044 'grid_transform_init maskgrid argument missing for maskinter transformation ')
2045 CALL raise_fatal_error()
2048 CALL get_val(in, nx=this%innx, ny=this%inny)
2049 .OR. IF (this%innx /= SIZE(maskgrid,1) this%inny /= SIZE(maskgrid,2)) THEN
2050 CALL l4f_category_log(this%category,L4F_ERROR, &
2051 'grid_transform_init mask non conformal with input field ')
2052 CALL l4f_category_log(this%category,L4F_ERROR, &
2053 'mask: '//t2c(SIZE(maskgrid,1))//'x '//t2c(SIZE(maskgrid,2))// &
2054 ' input field: '//t2c(this%innx)//'x '//t2c(this%inny))
2055 CALL raise_fatal_error()
2057 ! unlike before, here index arrays must have the shape of input grid
2058 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2059 this%inter_index_y(this%innx,this%inny))
2060 this%inter_index_x(:,:) = imiss
2061 this%inter_index_y(:,:) = 1
2063 ! generate the subarea boundaries according to maskgrid and maskbounds
2064 CALL gen_mask_class()
2066 ! compute coordinates of input grid in geo system
2070 !$OMP PARALLEL DEFAULT(SHARED)
2071 !$OMP DO PRIVATE(iy, ix, n)
2072 DO iy = 1, this%inny
2073 DO ix = 1, this%innx
2074 IF (c_e(maskgrid(ix,iy))) THEN
2075 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN
2076 DO n = nmaskarea, 1, -1
2077 IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN
2078 this%inter_index_x(ix,iy) = n
2088 this%outnx = nmaskarea
2090 CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
2091 CALL init(v7d_out, time_definition=time_definition)
2092 CALL vol7d_alloc(v7d_out, nana=nmaskarea)
2094 ! setup output point list, equal to average of polygon points
2095 ! warning, in case of concave areas points may coincide!
2097 CALL init(v7d_out%ana(n), &
2098 lon=stat_average(PACK(lin%dim%lon(:,:), &
2099 mask=(this%inter_index_x(:,:) == n))), &
2100 lat=stat_average(PACK(lin%dim%lat(:,:), &
2101 mask=(this%inter_index_x(:,:) == n))))
2105 this%valid = .TRUE. ! warning, no check of subtype
2107 ELSE IF (this%trans%trans_type == 'metamorphosis ') THEN
2109 ! common to all metamorphosis subtypes
2110 ! compute coordinates of input grid in geo system
2114 CALL get_val(in, nx=this%innx, ny=this%inny)
2115 ! allocate index array
2116 ALLOCATE(this%point_index(this%innx,this%inny))
2117 this%point_index(:,:) = imiss
2118 ! setup output coordinates
2119 CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
2120 CALL init(v7d_out, time_definition=time_definition)
2122 IF (this%trans%sub_type == 'all ' ) THEN
2124 this%outnx = this%innx*this%inny
2126 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2131 CALL init(v7d_out%ana((iy-1)*this%innx+ix), &
2132 lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2134 this%point_index(ix,iy) = n
2140 ELSE IF (this%trans%sub_type == 'coordbb ' ) THEN
2142 ! count and mark points falling into requested bounding-box
2145 DO iy = 1, this%inny
2146 DO ix = 1, this%innx
2147 ! IF (geo_coord_inside_rectang()
2148 .AND. IF (lin%dim%lon(ix,iy) > this%trans%rect_coo%ilon &
2149 .AND. lin%dim%lon(ix,iy) < this%trans%rect_coo%flon &
2150 .AND. lin%dim%lat(ix,iy) > this%trans%rect_coo%ilat &
2151 lin%dim%lat(ix,iy) < this%trans%rect_coo%flat) THEN ! improve!
2152 this%outnx = this%outnx + 1
2153 this%point_index(ix,iy) = this%outnx
2158 IF (this%outnx <= 0) THEN
2159 CALL l4f_category_log(this%category,L4F_WARN, &
2160 "metamorphosis:coordbb: no points inside bounding box "//&
2161 TRIM(to_char(this%trans%rect_coo%ilon))//","// &
2162 TRIM(to_char(this%trans%rect_coo%flon))//","// &
2163 TRIM(to_char(this%trans%rect_coo%ilat))//","// &
2164 TRIM(to_char(this%trans%rect_coo%flat)))
2167 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2169 ! collect coordinates of points falling into requested bounding-box
2171 DO iy = 1, this%inny
2172 DO ix = 1, this%innx
2173 IF (c_e(this%point_index(ix,iy))) THEN
2175 CALL init(v7d_out%ana(n), &
2176 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2183 ELSE IF (this%trans%sub_type == 'poly ' ) THEN
2185 ! count and mark points falling into requested polygon
2189 ! this OMP block has to be checked
2190 !$OMP PARALLEL DEFAULT(SHARED)
2191 !$OMP DO PRIVATE(iy, ix, point, n) REDUCTION(+:this%outnx)
2192 DO iy = 1, this%inny
2193 DO ix = 1, this%innx
2194 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
2195 DO n = 1, this%trans%poly%arraysize
2196 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
2197 this%outnx = this%outnx + 1
2198 this%point_index(ix,iy) = n
2202 ! CALL delete(point) ! speedup
2207 IF (this%outnx <= 0) THEN
2208 CALL l4f_category_log(this%category,L4F_WARN, &
2209 "metamorphosis:poly: no points inside polygons")
2212 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2213 ! collect coordinates of points falling into requested polygon
2215 DO iy = 1, this%inny
2216 DO ix = 1, this%innx
2217 IF (c_e(this%point_index(ix,iy))) THEN
2219 CALL init(v7d_out%ana(n), &
2220 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2227 ELSE IF (this%trans%sub_type == 'mask ' ) THEN
2229 .NOT. IF (PRESENT(maskgrid)) THEN
2230 CALL l4f_category_log(this%category,L4F_ERROR, &
2231 'grid_transform_init maskgrid argument missing for metamorphosis:mask ')
2236 .OR. IF (this%innx /= SIZE(maskgrid,1) this%inny /= SIZE(maskgrid,2)) THEN
2237 CALL l4f_category_log(this%category,L4F_ERROR, &
2238 'grid_transform_init mask non conformal with input field ')
2239 CALL l4f_category_log(this%category,L4F_ERROR, &
2240 'mask: '//t2c(SIZE(maskgrid,1))//'x '//t2c(SIZE(maskgrid,2))// &
2241 ' input field: '//t2c(this%innx)//'x '//t2c(this%inny))
2246 ! generate the subarea boundaries according to maskgrid and maskbounds
2247 CALL gen_mask_class()
2252 ! this OMP block has to be checked
2253 !$OMP PARALLEL DEFAULT(SHARED)
2254 !$OMP DO PRIVATE(iy, ix) REDUCTION(+:this%outnx)
2255 DO iy = 1, this%inny
2256 DO ix = 1, this%innx
2257 IF (c_e(maskgrid(ix,iy))) THEN
2258 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN
2259 DO n = nmaskarea, 1, -1
2260 IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN
2261 this%outnx = this%outnx + 1
2262 this%point_index(ix,iy) = n
2272 IF (this%outnx <= 0) THEN
2273 CALL l4f_category_log(this%category,L4F_WARN, &
2274 "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2278 CALL l4f_category_log(this%category,L4F_INFO, &
2279 "points in subarea "//t2c(n)//": "// &
2280 t2c(COUNT(this%point_index(:,:) == n)))
2283 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2284 ! collect coordinates of points falling into requested polygon
2286 DO iy = 1, this%inny
2287 DO ix = 1, this%innx
2288 IF (c_e(this%point_index(ix,iy))) THEN
2290 CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2303 SUBROUTINE gen_mask_class()
2304 REAL :: startmaskclass, mmin, mmax
2307 IF (PRESENT(maskbounds)) THEN
2308 nmaskarea = SIZE(maskbounds) - 1
2309 IF (nmaskarea > 0) THEN
2310 lmaskbounds = maskbounds ! automatic f2003 allocation
2312 ELSE IF (nmaskarea == 0) THEN
2313 CALL l4f_category_log(this%category,L4F_WARN, &
2314 'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2315 //', it will be ignored ')
2316 CALL l4f_category_log(this%category,L4F_WARN, &
2317 'at least 2 values are required for maskbounds ')
2321 mmin = MINVAL(maskgrid, mask=c_e(maskgrid))
2322 mmax = MAXVAL(maskgrid, mask=c_e(maskgrid))
2324 nmaskarea = INT(mmax - mmin + 1.5)
2325 startmaskclass = mmin - 0.5
2326 ! assign limits for each class
2327 ALLOCATE(lmaskbounds(nmaskarea+1))
2328 lmaskbounds(:) = (/(startmaskclass+REAL(i),i=0,nmaskarea)/)
2331 'maskinter, '//t2c(nmaskarea)//' subareas defined, with boundaries: ')
2332 DO i = 1, SIZE(lmaskbounds)
2333 CALL l4f_category_log(this%category,L4F_DEBUG, &
2334 'maskinter '//t2c(i)//' '//t2c(lmaskbounds(i)))
2338 END SUBROUTINE gen_mask_class
2340 END SUBROUTINE grid_transform_grid_vol7d_init
2343 !> Constructor for a \a grid_transform object, defining a particular
2344 !! sparse points-to-grid transformation.
2345 !! It defines an object describing a transformation from a set of
2346 !! sparse points to a rectangular grid; the abstract type of
2347 !! transformation is described in the transformation object \a trans
2348 !! (type transform_def) which must have been properly initialised. The
2349 !! additional information required here is the list of the input
2350 !! sparse points in the form of a \a vol7d object (parameter \a v7d_in),
2351 !! which can be the same volume that will be successively used for
2352 !! interpolation, or a volume with just the same coordinate data, and
2353 !! the description of the output grid \a griddim (a \a griddim_def
2356 !! The generated \a grid_transform object is specific to the sparse
2357 !! point list and grid provided. The function \a c_e can be used in
2358 !! order to check whether the object has been successfully
2359 !! initialised, if the result is \a .FALSE., it should not be used
2361 SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2362 TYPE(grid_transform),INTENT(out) :: this !< grid transformation object
2363 TYPE(transform_def),INTENT(in) :: trans !< transformation object
2364 TYPE(vol7d),INTENT(in) :: v7d_in !< vol7d object with the coordinates of the sparse point to be used as input (only information about coordinates is used)
2365 TYPE(griddim_def),INTENT(in) :: out !< griddim object defining target grid
2369 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2370 DOUBLE PRECISION,ALLOCATABLE :: lon(:,:),lat(:,:)
2371 TYPE(griddim_def) :: lout
2379 IF (this%trans%trans_type == 'inter ') THEN
2381 IF ( this%trans%sub_type == 'linear ' ) THEN
2383 this%innx=SIZE(v7d_in%ana)
2385 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2386 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2387 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2389 CALL copy (out, lout)
2390 ! equalize in/out coordinates
2391 lonref = 0.5D0*(MAXVAL(lon(:,1), mask=c_e(lon(:,1))) + MINVAL(lon(:,1), mask=c_e(lon(:,1))))
2392 CALL griddim_set_central_lon(lout, lonref)
2394 CALL get_val(lout, nx=nx, ny=ny)
2397 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2399 CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2400 CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
2401 CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)
2410 ELSE IF (this%trans%trans_type == 'boxinter ') THEN
2412 this%innx=SIZE(v7d_in%ana)
2414 ! index arrays must have the shape of input grid
2415 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2416 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2417 this%inter_index_y(this%innx,this%inny))
2418 ! get coordinates of input grid in geo system
2419 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2421 CALL copy (out, lout)
2422 ! equalize in/out coordinates
2423 lonref = 0.5D0*(MAXVAL(lon(:,1), mask=c_e(lon(:,1))) + MINVAL(lon(:,1), mask=c_e(lon(:,1))))
2424 CALL griddim_set_central_lon(lout, lonref)
2426 CALL get_val(lout, nx=this%outnx, ny=this%outny, &
2427 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2428 ! TODO now box size is ignored
2429 ! if box size not provided, use the actual grid step
2430 .NOT. IF (c_e(this%trans%area_info%boxdx)) &
2431 CALL get_val(out, dx=this%trans%area_info%boxdx)
2432 .NOT. IF (c_e(this%trans%area_info%boxdy)) &
2433 CALL get_val(out, dx=this%trans%area_info%boxdy)
2434 ! half size is actually needed
2435 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5D0
2436 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5D0
2438 ! use find_index in the opposite way, here extrap does not make sense
2439 CALL this%find_index(lout, .TRUE., &
2440 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2441 lon, lat, .FALSE., &
2442 this%inter_index_x, this%inter_index_y)
2447 this%valid = .TRUE. ! warning, no check of subtype
2451 END SUBROUTINE grid_transform_vol7d_grid_init
2454 !> Constructor for a \a grid_transform object, defining a particular
2455 !! sparse points-to-sparse points transformation.
2456 !! It defines an object describing a transformation from a set of
2457 !! sparse points to a set of sparse points; the abstract type of
2458 !! transformation is described in the transformation object \a trans
2459 !! (type transform_def) which must have been properly initialised. The
2460 !! additional information required here is the list of the input
2461 !! sparse points in the form of a \a vol7d object (parameter \a
2462 !! v7d_in), which can be the same volume that will be successively
2463 !! used for interpolation, or a volume with just the same coordinate
2464 !! data, and, if required by the transformation type, the information
2465 !! about the target sparse points over which the transformation should
2468 !! - for 'inter ' transformation, this is provided in the form of a
2469 !! vol7d object (\a v7d_out argument, input), which must have been
2470 !! initialized with the coordinates of desired sparse points
2472 !! - for 'polyinter ' transformation, no target point information has
2473 !! to be provided in input (it is calculated on the basis of input
2474 !! grid and \a trans object), and the coordinates of the target
2475 !! points (polygons' centroids) are returned in output in \a
2488 SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2489 maskbounds, categoryappend)
2490 TYPE(grid_transform), INTENT(out) :: this
2491 TYPE(transform_def), INTENT(in) :: trans
2492 TYPE(vol7d), INTENT(in) :: v7d_in
2493 TYPE(vol7d), INTENT(inout) :: v7d_out
2494 REAL, INTENT(in), OPTIONAL :: maskbounds(:)
2498 DOUBLE PRECISION, ALLOCATABLE :: lon(:), lat(:)
2500 DOUBLE PRECISION :: lon1, lat1
2501 TYPE(georef_coord) :: point
2509 IF (this%trans%trans_type == 'inter') THEN
2511 IF ( this%trans%sub_type == 'linear' ) THEN
2513 CALL vol7d_alloc_vol(v7d_out)
2514 this%outnx= SIZE(v7d_out%ana)
2517 this%innx= SIZE(v7d_in%ana)
2520 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny
2521 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx
2523 CALL getval(v7d_in%ana(:)%coord,lon=this%inter_xp(:,1),lat=this%inter_yp
2524 CALL getval(v7d_out%ana(:)%coord,lon=this%inter_x(:,1),lat=this%inter_y
|