|
◆ 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 2032 del file grid_transform_class.F90.
2034 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
2035 this%inter_index_x(ix,iy) > xnmax .OR. &
2036 this%inter_index_y(ix,iy) < ynmin .OR. &
2037 this%inter_index_y(ix,iy) > ynmax) THEN
2038 this%inter_index_x(ix,iy) = imiss
2039 this%inter_index_y(ix,iy) = imiss
2045 CALL l4f_category_log(this%category, l4f_debug, &
2046 'stencilinter: stencil size '//t2c(n*n))
2047 CALL l4f_category_log(this%category, l4f_debug, &
2048 'stencilinter: stencil points '//t2c(count(this%stencil)))
2056 ELSE IF (this%trans%trans_type == 'maskinter') THEN
2058 IF (.NOT. PRESENT(maskgrid)) THEN
2059 CALL l4f_category_log(this%category,l4f_error, &
2060 'grid_transform_init maskgrid argument missing for maskinter transformation ')
2061 CALL raise_fatal_error()
2064 CALL get_val(in, nx=this%innx, ny=this%inny)
2065 .OR. IF (this%innx /= SIZE(maskgrid,1) this%inny /= SIZE(maskgrid,2)) THEN
2066 CALL l4f_category_log(this%category,L4F_ERROR, &
2067 'grid_transform_init mask non conformal with input field ')
2068 CALL l4f_category_log(this%category,L4F_ERROR, &
2069 'mask: '//t2c(SIZE(maskgrid,1))//'x '//t2c(SIZE(maskgrid,2))// &
2070 ' input field: '//t2c(this%innx)//'x '//t2c(this%inny))
2071 CALL raise_fatal_error()
2073 ! unlike before, here index arrays must have the shape of input grid
2074 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2075 this%inter_index_y(this%innx,this%inny))
2076 this%inter_index_x(:,:) = imiss
2077 this%inter_index_y(:,:) = 1
2079 ! generate the subarea boundaries according to maskgrid and maskbounds
2080 CALL gen_mask_class()
2082 ! compute coordinates of input grid in geo system
2086 !$OMP PARALLEL DEFAULT(SHARED)
2087 !$OMP DO PRIVATE(iy, ix, n)
2088 DO iy = 1, this%inny
2089 DO ix = 1, this%innx
2090 IF (c_e(maskgrid(ix,iy))) THEN
2091 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN
2092 DO n = nmaskarea, 1, -1
2093 IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN
2094 this%inter_index_x(ix,iy) = n
2104 this%outnx = nmaskarea
2106 CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
2107 CALL init(v7d_out, time_definition=time_definition)
2108 CALL vol7d_alloc(v7d_out, nana=nmaskarea)
2110 ! setup output point list, equal to average of polygon points
2111 ! warning, in case of concave areas points may coincide!
2113 CALL init(v7d_out%ana(n), &
2114 lon=stat_average(PACK(lin%dim%lon(:,:), &
2115 mask=(this%inter_index_x(:,:) == n))), &
2116 lat=stat_average(PACK(lin%dim%lat(:,:), &
2117 mask=(this%inter_index_x(:,:) == n))))
2121 this%valid = .TRUE. ! warning, no check of subtype
2123 ELSE IF (this%trans%trans_type == 'metamorphosis ') THEN
2125 ! common to all metamorphosis subtypes
2126 ! compute coordinates of input grid in geo system
2130 CALL get_val(in, nx=this%innx, ny=this%inny)
2131 ! allocate index array
2132 ALLOCATE(this%point_index(this%innx,this%inny))
2133 this%point_index(:,:) = imiss
2134 ! setup output coordinates
2135 CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
2136 CALL init(v7d_out, time_definition=time_definition)
2138 IF (this%trans%sub_type == 'all ' ) THEN
2140 this%outnx = this%innx*this%inny
2142 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2147 CALL init(v7d_out%ana((iy-1)*this%innx+ix), &
2148 lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2150 this%point_index(ix,iy) = n
2156 ELSE IF (this%trans%sub_type == 'coordbb ' ) THEN
2158 ! count and mark points falling into requested bounding-box
2161 DO iy = 1, this%inny
2162 DO ix = 1, this%innx
2163 ! IF (geo_coord_inside_rectang()
2164 .AND. IF (lin%dim%lon(ix,iy) > this%trans%rect_coo%ilon &
2165 .AND. lin%dim%lon(ix,iy) < this%trans%rect_coo%flon &
2166 .AND. lin%dim%lat(ix,iy) > this%trans%rect_coo%ilat &
2167 lin%dim%lat(ix,iy) < this%trans%rect_coo%flat) THEN ! improve!
2168 this%outnx = this%outnx + 1
2169 this%point_index(ix,iy) = this%outnx
2174 IF (this%outnx <= 0) THEN
2175 CALL l4f_category_log(this%category,L4F_WARN, &
2176 "metamorphosis:coordbb: no points inside bounding box "//&
2177 TRIM(to_char(this%trans%rect_coo%ilon))//","// &
2178 TRIM(to_char(this%trans%rect_coo%flon))//","// &
2179 TRIM(to_char(this%trans%rect_coo%ilat))//","// &
2180 TRIM(to_char(this%trans%rect_coo%flat)))
2183 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2185 ! collect coordinates of points falling into requested bounding-box
2187 DO iy = 1, this%inny
2188 DO ix = 1, this%innx
2189 IF (c_e(this%point_index(ix,iy))) THEN
2191 CALL init(v7d_out%ana(n), &
2192 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2199 ELSE IF (this%trans%sub_type == 'poly ' ) THEN
2201 ! count and mark points falling into requested polygon
2205 ! this OMP block has to be checked
2206 !$OMP PARALLEL DEFAULT(SHARED)
2207 !$OMP DO PRIVATE(iy, ix, point, n) REDUCTION(+:this%outnx)
2208 DO iy = 1, this%inny
2209 DO ix = 1, this%innx
2210 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
2211 DO n = 1, this%trans%poly%arraysize
2212 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
2213 this%outnx = this%outnx + 1
2214 this%point_index(ix,iy) = n
2218 ! CALL delete(point) ! speedup
2223 IF (this%outnx <= 0) THEN
2224 CALL l4f_category_log(this%category,L4F_WARN, &
2225 "metamorphosis:poly: no points inside polygons")
2228 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2229 ! collect coordinates of points falling into requested polygon
2231 DO iy = 1, this%inny
2232 DO ix = 1, this%innx
2233 IF (c_e(this%point_index(ix,iy))) THEN
2235 CALL init(v7d_out%ana(n), &
2236 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2243 ELSE IF (this%trans%sub_type == 'mask ' ) THEN
2245 .NOT. IF (PRESENT(maskgrid)) THEN
2246 CALL l4f_category_log(this%category,L4F_ERROR, &
2247 'grid_transform_init maskgrid argument missing for metamorphosis:mask ')
2252 .OR. IF (this%innx /= SIZE(maskgrid,1) this%inny /= SIZE(maskgrid,2)) THEN
2253 CALL l4f_category_log(this%category,L4F_ERROR, &
2254 'grid_transform_init mask non conformal with input field ')
2255 CALL l4f_category_log(this%category,L4F_ERROR, &
2256 'mask: '//t2c(SIZE(maskgrid,1))//'x '//t2c(SIZE(maskgrid,2))// &
2257 ' input field: '//t2c(this%innx)//'x '//t2c(this%inny))
2262 ! generate the subarea boundaries according to maskgrid and maskbounds
2263 CALL gen_mask_class()
2268 ! this OMP block has to be checked
2269 !$OMP PARALLEL DEFAULT(SHARED)
2270 !$OMP DO PRIVATE(iy, ix) REDUCTION(+:this%outnx)
2271 DO iy = 1, this%inny
2272 DO ix = 1, this%innx
2273 IF (c_e(maskgrid(ix,iy))) THEN
2274 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN
2275 DO n = nmaskarea, 1, -1
2276 IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN
2277 this%outnx = this%outnx + 1
2278 this%point_index(ix,iy) = n
2288 IF (this%outnx <= 0) THEN
2289 CALL l4f_category_log(this%category,L4F_WARN, &
2290 "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2294 CALL l4f_category_log(this%category,L4F_INFO, &
2295 "points in subarea "//t2c(n)//": "// &
2296 t2c(COUNT(this%point_index(:,:) == n)))
2299 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2300 ! collect coordinates of points falling into requested polygon
2302 DO iy = 1, this%inny
2303 DO ix = 1, this%innx
2304 IF (c_e(this%point_index(ix,iy))) THEN
2306 CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2319 SUBROUTINE gen_mask_class()
2320 REAL :: startmaskclass, mmin, mmax
2323 IF (PRESENT(maskbounds)) THEN
2324 nmaskarea = SIZE(maskbounds) - 1
2325 IF (nmaskarea > 0) THEN
2326 lmaskbounds = maskbounds ! automatic f2003 allocation
2328 ELSE IF (nmaskarea == 0) THEN
2329 CALL l4f_category_log(this%category,L4F_WARN, &
2330 'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2331 //', it will be ignored ')
2332 CALL l4f_category_log(this%category,L4F_WARN, &
2333 'at least 2 values are required for maskbounds ')
2337 mmin = MINVAL(maskgrid, mask=c_e(maskgrid))
2338 mmax = MAXVAL(maskgrid, mask=c_e(maskgrid))
2340 nmaskarea = INT(mmax - mmin + 1.5)
2341 startmaskclass = mmin - 0.5
2342 ! assign limits for each class
2343 ALLOCATE(lmaskbounds(nmaskarea+1))
2344 lmaskbounds(:) = (/(startmaskclass+REAL(i),i=0,nmaskarea)/)
2347 'maskinter, '//t2c(nmaskarea)//' subareas defined, with boundaries: ')
2348 DO i = 1, SIZE(lmaskbounds)
2349 CALL l4f_category_log(this%category,L4F_DEBUG, &
2350 'maskinter '//t2c(i)//' '//t2c(lmaskbounds(i)))
2354 END SUBROUTINE gen_mask_class
2356 END SUBROUTINE grid_transform_grid_vol7d_init
2359 !> Constructor for a \a grid_transform object, defining a particular
2360 !! sparse points-to-grid transformation.
2361 !! It defines an object describing a transformation from a set of
2362 !! sparse points to a rectangular grid; the abstract type of
2363 !! transformation is described in the transformation object \a trans
2364 !! (type transform_def) which must have been properly initialised. The
2365 !! additional information required here is the list of the input
2366 !! sparse points in the form of a \a vol7d object (parameter \a v7d_in),
2367 !! which can be the same volume that will be successively used for
2368 !! interpolation, or a volume with just the same coordinate data, and
2369 !! the description of the output grid \a griddim (a \a griddim_def
2372 !! The generated \a grid_transform object is specific to the sparse
2373 !! point list and grid provided. The function \a c_e can be used in
2374 !! order to check whether the object has been successfully
2375 !! initialised, if the result is \a .FALSE., it should not be used
2377 SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2378 TYPE(grid_transform),INTENT(out) :: this !< grid transformation object
2379 TYPE(transform_def),INTENT(in) :: trans !< transformation object
2380 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)
2381 TYPE(griddim_def),INTENT(in) :: out !< griddim object defining target grid
2385 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2386 DOUBLE PRECISION,ALLOCATABLE :: lon(:,:),lat(:,:)
2387 TYPE(griddim_def) :: lout
2395 IF (this%trans%trans_type == 'inter ') THEN
2397 IF ( this%trans%sub_type == 'linear ' ) THEN
2399 this%innx=SIZE(v7d_in%ana)
2401 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2402 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2403 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2405 CALL copy (out, lout)
2406 ! equalize in/out coordinates
2407 lonref = 0.5D0*(MAXVAL(lon(:,1), mask=c_e(lon(:,1))) + MINVAL(lon(:,1), mask=c_e(lon(:,1))))
2408 CALL griddim_set_central_lon(lout, lonref)
2410 CALL get_val(lout, nx=nx, ny=ny)
2413 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2415 CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2416 CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
2417 CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)
2426 ELSE IF (this%trans%trans_type == 'boxinter ') THEN
2428 this%innx=SIZE(v7d_in%ana)
2430 ! index arrays must have the shape of input grid
2431 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2432 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2433 this%inter_index_y(this%innx,this%inny))
2434 ! get coordinates of input grid in geo system
2435 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2437 CALL copy (out, lout)
2438 ! equalize in/out coordinates
2439 lonref = 0.5D0*(MAXVAL(lon(:,1), mask=c_e(lon(:,1))) + MINVAL(lon(:,1), mask=c_e(lon(:,1))))
2440 CALL griddim_set_central_lon(lout, lonref)
2442 CALL get_val(lout, nx=this%outnx, ny=this%outny, &
2443 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2444 ! TODO now box size is ignored
2445 ! if box size not provided, use the actual grid step
2446 .NOT. IF (c_e(this%trans%area_info%boxdx)) &
2447 CALL get_val(out, dx=this%trans%area_info%boxdx)
2448 .NOT. IF (c_e(this%trans%area_info%boxdy)) &
2449 CALL get_val(out, dx=this%trans%area_info%boxdy)
2450 ! half size is actually needed
2451 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5D0
2452 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5D0
2454 ! use find_index in the opposite way, here extrap does not make sense
2455 CALL this%find_index(lout, .TRUE., &
2456 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2457 lon, lat, .FALSE., &
2458 this%inter_index_x, this%inter_index_y)
2463 this%valid = .TRUE. ! warning, no check of subtype
2467 END SUBROUTINE grid_transform_vol7d_grid_init
2470 !> Constructor for a \a grid_transform object, defining a particular
2471 !! sparse points-to-sparse points transformation.
2472 !! It defines an object describing a transformation from a set of
2473 !! sparse points to a set of sparse points; the abstract type of
2474 !! transformation is described in the transformation object \a trans
2475 !! (type transform_def) which must have been properly initialised. The
2476 !! additional information required here is the list of the input
2477 !! sparse points in the form of a \a vol7d object (parameter \a
2478 !! v7d_in), which can be the same volume that will be successively
2479 !! used for interpolation, or a volume with just the same coordinate
2480 !! data, and, if required by the transformation type, the information
2481 !! about the target sparse points over which the transformation should
2484 !! - for 'inter ' transformation, this is provided in the form of a
2485 !! vol7d object (\a v7d_out argument, input), which must have been
2486 !! initialized with the coordinates of desired sparse points
2488 !! - for 'polyinter ' transformation, no target point information has
2489 !! to be provided in input (it is calculated on the basis of input
2490 !! grid and \a trans object), and the coordinates of the target
2491 !! points (polygons' centroids) are returned in output in \a
2504 SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2505 maskbounds, categoryappend)
2506 TYPE(grid_transform), INTENT(out) :: this
2507 TYPE(transform_def), INTENT(in) :: trans
2508 TYPE(vol7d), INTENT(in) :: v7d_in
2509 TYPE(vol7d), INTENT(inout) :: v7d_out
2510 REAL, INTENT(in), OPTIONAL :: maskbounds(:)
2514 DOUBLE PRECISION, ALLOCATABLE :: lon(:), lat(:)
2516 DOUBLE PRECISION :: lon1, lat1
2517 TYPE(georef_coord) :: point
2525 IF (this%trans%trans_type == 'inter') THEN
2527 IF ( this%trans%sub_type == 'linear' ) THEN
2529 CALL vol7d_alloc_vol(v7d_out)
2530 this%outnx= SIZE(v7d_out%ana)
2533 this%innx= SIZE(v7d_in%ana)
2536 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny
2537 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx
2539 CALL getval(v7d_in%ana(:)%coord,lon=this%inter_xp(:,1),lat=this%inter_yp
2540 CALL getval(v7d_out%ana(:)%coord,lon=this%inter_x(:,1),lat=this%inter_y
|