|
◆ grid_transform_init()
subroutine grid_transform_init |
( |
type(grid_transform), intent(out) |
this, |
|
|
type(transform_def), intent(in) |
trans, |
|
|
type(griddim_def), intent(inout) |
in, |
|
|
type(griddim_def), intent(inout) |
out, |
|
|
real, dimension(:,:), intent(in), optional |
maskgrid, |
|
|
real, dimension(:), intent(in), optional |
maskbounds, |
|
|
character(len=*), intent(in), optional |
categoryappend |
|
) |
| |
Constructor for a grid_transform object, defining a particular grid-to-grid transformation.
It defines an object describing a transformation from one rectangular grid to another; 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), the description of the output grid out (type griddim_def as well). The description of the output grid must be initialized for interpolating type transformations ('inter' and 'boxinter'), while it is generated by this constructor and returned in output for 'zoom', 'boxregrid', 'maskgen' and 'polyinter' transformations.
The generated grid_transform object is specific to the input and output grids involved. 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,out] | in | griddim object to transform |
[in,out] | out | griddim object defining target grid (input or output depending on type of transformation) |
[in] | maskgrid | 2D field to be used for defining valid points, it must have the same shape as the field to be interpolated (for transformation type 'metamorphosis:maskvalid') |
[in] | maskbounds | array of boundary values for defining a subset of valid points where the values of maskgrid are within the first and last value of maskbounds (for transformation type 'metamorphosis:maskvalid/settoinvalid' and others) |
[in] | categoryappend | append this suffix to log4fortran namespace category |
Definizione alla linea 1415 del file grid_transform_class.F90.
1417 this%inter_yp(this%outnx,this%outny))
1420 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1422 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp
1425 CALL this%find_index(in, .true., &
1426 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1427 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1428 this%inter_index_x, this%inter_index_y)
1436 ELSE IF (this%trans%trans_type == 'boxinter') THEN
1438 CALL outgrid_setup()
1439 CALL get_val(in, nx=this%innx, ny=this%inny)
1440 CALL get_val(out, nx=this%outnx, ny=this%outny, &
1441 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1444 IF (.NOT.c_e(this%trans%area_info%boxdx)) &
1445 CALL get_val(out, dx=this%trans%area_info%boxdx)
1446 IF (.NOT.c_e(this%trans%area_info%boxdy)) &
1447 CALL get_val(out, dx=this%trans%area_info%boxdy)
1449 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
1450 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
1452 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1453 this%inter_index_y(this%innx,this%inny))
1459 CALL this%find_index(out, .true., &
1460 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
1461 lin%dim%lon, lin%dim%lat, .false., &
1462 this%inter_index_x, this%inter_index_y)
1467 ELSE IF (this%trans%trans_type == 'stencilinter') THEN
1469 CALL outgrid_setup()
1471 CALL get_val(in, nx=this%innx, ny=this%inny, &
1472 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1473 CALL get_val(out, nx=this%outnx, ny=this%outny)
1475 ALLOCATE (this%inter_index_x(this%outnx,this%outny), &
1476 this%inter_index_y(this%outnx,this%outny))
1477 CALL copy(out, lout)
1479 CALL this%find_index(in, .true., &
1480 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1481 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1482 this%inter_index_x, this%inter_index_y)
1485 nr = int(this%trans%area_info%radius)
1488 r2 = this%trans%area_info%radius**2
1489 ALLOCATE(this%stencil(n,n))
1490 this%stencil(:,:) = .true.
1493 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1500 xnmax = this%innx - nr
1502 ynmax = this%inny - nr
1503 DO iy = 1, this%outny
1504 DO ix = 1, this%outnx
1505 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
1506 this%inter_index_x(ix,iy) > xnmax .OR. &
1507 this%inter_index_y(ix,iy) < ynmin .OR. &
1508 this%inter_index_y(ix,iy) > ynmax) THEN
1509 this%inter_index_x(ix,iy) = imiss
1510 this%inter_index_y(ix,iy) = imiss
1516 CALL l4f_category_log(this%category, l4f_debug, &
1517 'stencilinter: stencil size '//t2c(n*n))
1518 CALL l4f_category_log(this%category, l4f_debug, &
1519 'stencilinter: stencil points '//t2c(count(this%stencil)))
1525 ELSE IF (this%trans%trans_type == 'maskgen') THEN
1527 IF (this%trans%sub_type == 'poly') THEN
1530 CALL get_val(in, nx=this%innx, ny=this%inny)
1531 this%outnx = this%innx
1532 this%outny = this%inny
1535 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1536 this%inter_index_y(this%innx,this%inny))
1537 this%inter_index_x(:,:) = imiss
1538 this%inter_index_y(:,:) = 1
1547 inside_x: DO i = 1, this%innx
1548 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1550 DO n = nprev, this%trans%poly%arraysize
1551 IF (inside(point, this%trans%poly%array(n))) THEN
1552 this%inter_index_x(i,j) = n
1557 DO n = nprev-1, 1, -1
1558 IF (inside(point, this%trans%poly%array(n))) THEN
1559 this%inter_index_x(i,j) = n
1570 ELSE IF (this%trans%sub_type == 'grid') THEN
1573 CALL copy(out, lout)
1576 CALL get_val(in, nx=this%innx, ny=this%inny)
1577 this%outnx = this%innx
1578 this%outny = this%inny
1579 CALL get_val(lout, nx=nx, ny=ny, &
1580 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1583 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1584 this%inter_index_y(this%innx,this%inny))
1590 CALL this%find_index(lout, .true., &
1591 nx, ny, xmin, xmax, ymin, ymax, &
1592 out%dim%lon, out%dim%lat, .false., &
1593 this%inter_index_x, this%inter_index_y)
1595 WHERE(c_e(this%inter_index_x(:,:)))
1596 this%inter_index_x(:,:) = this%inter_index_x(:,:) + &
1597 (this%inter_index_y(:,:)-1)*nx
1605 ELSE IF (this%trans%trans_type == 'polyinter') THEN
1611 CALL get_val(in, nx=this%innx, ny=this%inny)
1612 this%outnx = this%innx
1613 this%outny = this%inny
1616 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1617 this%inter_index_y(this%innx,this%inny))
1618 this%inter_index_x(:,:) = imiss
1619 this%inter_index_y(:,:) = 1
1628 inside_x_2: DO i = 1, this%innx
1629 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1631 DO n = nprev, this%trans%poly%arraysize
1632 IF (inside(point, this%trans%poly%array(n))) THEN
1633 this%inter_index_x(i,j) = n
1638 DO n = nprev-1, 1, -1
1639 IF (inside(point, this%trans%poly%array(n))) THEN
1640 this%inter_index_x(i,j) = n
1653 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
1656 CALL get_val(in, nx=this%innx, ny=this%inny)
1657 this%outnx = this%innx
1658 this%outny = this%inny
1660 IF (this%trans%sub_type == 'maskvalid' .OR. this%trans%sub_type == 'maskinvalid'THEN
1662 IF (.NOT. PRESENT(maskgrid)) THEN
1663 CALL l4f_category_log(this%category,l4f_error, &
1664 'grid_transform_init maskgrid argument missing for metamorphosis:'
1665 trim(this%trans%sub_type)// ' transformation')
1670 IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2) THEN
1671 CALL l4f_category_log(this%category,l4f_error, &
1672 'grid_transform_init mask non conformal with input field')
1673 CALL l4f_category_log(this%category,l4f_error, &
1674 'mask: '//t2c( SIZE(maskgrid,1))// 'x'//t2c( SIZE(maskgrid,2))// &
1675 ' input field:'//t2c(this%innx)// 'x'//t2c(this%inny))
1680 ALLOCATE(this%point_mask(this%innx,this%inny))
1682 IF (this%trans%sub_type == 'maskvalid') THEN
1685 IF (.NOT. PRESENT(maskbounds)) THEN
1686 this%point_mask(:,:) = c_e(maskgrid(:,:))
1687 ELSE IF ( SIZE(maskbounds) < 2) THEN
1688 this%point_mask(:,:) = c_e(maskgrid(:,:))
1690 this%point_mask(:,:) = c_e(maskgrid(:,:)) .AND. &
1691 maskgrid(:,:) > maskbounds(1) .AND. &
1692 maskgrid(:,:) <= maskbounds( SIZE(maskbounds))
1695 this%point_mask(:,:) = .NOT.c_e(maskgrid(:,:))
1700 ELSE IF (this%trans%sub_type == 'lemaskinvalid' .OR. &
1701 this%trans%sub_type == 'ltmaskinvalid' .OR. &
1702 this%trans%sub_type == 'gemaskinvalid' .OR. &
1703 this%trans%sub_type == 'gtmaskinvalid') THEN
1706 this%val_mask = maskgrid
1709 ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
1711 IF (.NOT. PRESENT(maskbounds)) THEN
1712 CALL l4f_category_log(this%category,l4f_error, &
1713 'grid_transform_init maskbounds missing for metamorphosis:'// &
1714 trim(this%trans%sub_type)// ' transformation')
1716 ELSE IF ( SIZE(maskbounds) < 1) THEN
1717 CALL l4f_category_log(this%category,l4f_error, &
1718 'grid_transform_init maskbounds empty for metamorphosis:'// &
1719 trim(this%trans%sub_type)// ' transformation')
1722 this%val1 = maskbounds(1)
1724 CALL l4f_category_log(this%category, l4f_debug, &
1725 "grid_transform_init setting invalid data to "//t2c(this%val1))
1731 ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
1733 IF (.NOT. PRESENT(maskbounds)) THEN
1734 CALL l4f_category_log(this%category,l4f_error, &
1735 'grid_transform_init maskbounds missing for metamorphosis:'// &
1736 trim(this%trans%sub_type)// ' transformation')
1739 ELSE IF ( SIZE(maskbounds) < 2) THEN
1740 CALL l4f_category_log(this%category,l4f_error, &
1741 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'
1742 trim(this%trans%sub_type)// ' transformation')
1746 this%val1 = maskbounds(1)
1747 this%val2 = maskbounds( SIZE(maskbounds))
1749 CALL l4f_category_log(this%category, l4f_debug, &
1750 "grid_transform_init setting to invalid interval ]"//t2c(this%val1 ','
1751 t2c(this%val2)// ']')
1765 SUBROUTINE outgrid_setup()
1772 IF (proj_in == proj_out) THEN
1775 CALL set_val(out, component_flag=cf_i)
1780 CALL l4f_category_log(this%category,l4f_warn, &
1781 'trying to interpolate a grid with component flag 1 to a grid on a different ')
1782 CALL l4f_category_log(this%category,L4F_WARN, &
1783 'vector fields will probably be wrong ')
1785 CALL set_val(out, component_flag=cf_i)
1788 ! rotate the input grid by n*360 degrees to bring it closer to the output grid
1791 END SUBROUTINE outgrid_setup
1793 END SUBROUTINE grid_transform_init
1796 !> Constructor for a \a grid_transform object, defining a particular
1797 !! grid-to-sparse points transformation.
1798 !! It defines an object describing a transformation from a rectangular
1799 !! grid to a set of sparse points; the abstract type of transformation
1800 !! is described in the transformation object \a trans (type
1801 !! transform_def) which must have been properly initialised. The
1802 !! additional information required here is the description of the
1803 !! input grid \a in (type griddim_def), and, if required by the
1804 !! transformation type, the information about the target sparse points
1805 !! over which the transformation should take place:
1807 !! - for 'inter ' transformation, this is provided in the form of a
1808 !! vol7d object (\a v7d_out argument, input), which must have been
1809 !! initialized with the coordinates of desired sparse points
1811 !! - for 'polyinter ' transformation, no target point information has
1812 !! to be provided in input (it is calculated on the basis of input
1813 !! grid and \a trans object), and the coordinates of the target
1814 !! points (polygons' centroids) are returned in output in \a
1838 SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1839 maskgrid, maskbounds, find_index, categoryappend)
1840 TYPE(grid_transform), INTENT(out) :: this
1841 TYPE(transform_def), INTENT(in) :: trans
1842 TYPE(griddim_def), INTENT(in) :: in
1843 TYPE(vol7d), INTENT(inout) :: v7d_out
1844 REAL, INTENT(in), OPTIONAL :: maskgrid(:,:)
1845 REAL, INTENT(in), OPTIONAL :: maskbounds(:)
1846 PROCEDURE(basic_find_index), POINTER, OPTIONAL :: find_index
1849 INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax
1851 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
1852 DOUBLE PRECISION, ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
1853 REAL, ALLOCATABLE :: lmaskbounds(:)
1854 TYPE(georef_coord) :: point
1855 TYPE(griddim_def) :: lin
1858 IF ( PRESENT(find_index)) THEN
1859 IF ( ASSOCIATED(find_index)) THEN
1860 this%find_index => find_index
1870 IF (.NOT. c_e(time_definition)) THEN
1874 IF (this%trans%trans_type == 'inter') THEN
1876 IF (this%trans%sub_type == 'near' .OR. this%trans%sub_type == 'bilin'
1877 .OR. this%trans%sub_type == 'shapiro_near') THEN
1880 CALL get_val(lin, nx=this%innx, ny=this%inny)
1881 this%outnx = SIZE(v7d_out%ana)
1884 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1885 this%inter_index_y(this%outnx,this%outny))
1886 ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
1888 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
1890 lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,
1891 CALL griddim_set_central_lon(lin, lonref)
1892 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1894 IF (this%trans%sub_type == 'bilin') THEN
1895 CALL this%find_index(lin, .false., &
1896 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1897 lon, lat, this%trans%extrap, &
1898 this%inter_index_x, this%inter_index_y)
1900 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny
1901 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx
1903 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1904 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1907 CALL this%find_index(lin, .true., &
1908 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1909 lon, lat, this%trans%extrap, &
1910 this%inter_index_x, this%inter_index_y)
1921 ELSE IF (this%trans%trans_type == 'polyinter') THEN
1923 CALL get_val(in, nx=this%innx, ny=this%inny)
1925 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1926 this%inter_index_y(this%innx,this%inny))
1927 this%inter_index_x(:,:) = imiss
1928 this%inter_index_y(:,:) = 1
1934 this%outnx = this%trans%poly%arraysize
1936 CALL delete(v7d_out)
1937 CALL init(v7d_out, time_definition=time_definition)
1938 CALL vol7d_alloc(v7d_out, nana=this%outnx)
1941 ALLOCATE(lon(this%outnx,1))
1942 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1))
1943 lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1)
1944 CALL griddim_set_central_lon(lin, lonref)
1949 CALL poly_to_coordinates(this%trans%poly, v7d_out)
1954 DO iy = 1, this%inny
1955 inside_x: DO ix = 1, this%innx
1956 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy
1958 DO n = nprev, this%trans%poly%arraysize
1959 IF (inside(point, this%trans%poly%array(n))) THEN
1960 this%inter_index_x(ix,iy) = n
1965 DO n = nprev-1, 1, -1
1966 IF (inside(point, this%trans%poly%array(n))) THEN
1967 this%inter_index_x(ix,iy) = n
1979 DO n = 1, this%trans%poly%arraysize
1980 CALL l4f_category_log(this%category, l4f_debug, &
1981 'polygon: '//t2c(n)//' grid points: '// &
1982 t2c(COUNT(this%inter_index_x(:,:) == n)))
|