libsim Versione 7.2.6
|
◆ grid_transform_v7d_grid_compute()
Compute the output data array from input data array according to the defined transformation. The grid_transform object this must have been properly initialised, so that it contains all the information needed for computing the transformation. This is the sparse points-to-grid and sparse points-to-sparse points version.
Definizione alla linea 4326 del file grid_transform_class.F90. 4328
4329z5 = (z4-z1)*p2+z1
4330z6 = (z3-z2)*p2+z2
4331
4332zp = (z6-z5)*(p1)+z5
4333
4334END FUNCTION hbilin
4335
4336
4337! Shapiro filter of order 2
4338FUNCTION shapiro (z,zp) RESULT(zs)
4339!! z_smoothed(i,j) = z(i,j) * (1-S) + S * sum(z_vicini)/N
4340!! = z(i,j) - S/N (N*z(i,j) - sum(z_vicini))
4341REAL,INTENT(in) :: z(:) ! Z values on the four surrounding points
4342REAL,INTENT(in) :: zp ! Z values on the central point
4343REAL :: zs ! Z smoothed value on the central point
4344INTEGER:: N
4345
4346IF(c_e(zp))THEN
4347 n = count(c_e(z))
4348 zs = zp - 0.5* ( n*zp - sum(z, c_e(z)) )/n
4349ELSE
4350 zs = rmiss
4351END IF
4352
4353END FUNCTION shapiro
4354
4355
4356! Locate index of requested point
4357SUBROUTINE basic_find_index(this, near, nx, ny, xmin, xmax, ymin, ymax, &
4358 lon, lat, extrap, index_x, index_y)
4359TYPE(griddim_def),INTENT(in) :: this ! griddim object (from grid)
4360LOGICAL,INTENT(in) :: near ! near or bilin interpolation (determine wich point is requested)
4361INTEGER,INTENT(in) :: nx,ny ! dimension (to grid)
4362DOUBLE PRECISION,INTENT(in) :: xmin, xmax, ymin, ymax ! extreme coordinate (to grid)
4363DOUBLE PRECISION,INTENT(in) :: lon(:,:),lat(:,:) ! target coordinate
4364LOGICAL,INTENT(in) :: extrap ! extrapolate
4365INTEGER,INTENT(out) :: index_x(:,:),index_y(:,:) ! index of point requested
4366
4367INTEGER :: lnx, lny
4368DOUBLE PRECISION :: x(SIZE(lon,1),SIZE(lon,2)),y(SIZE(lon,1),SIZE(lon,2))
4369
4370IF (near) THEN
4371 CALL proj(this,lon,lat,x,y)
4372 index_x = nint((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4373 index_y = nint((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4374 lnx = nx
4375 lny = ny
4376ELSE
4377 CALL proj(this,lon,lat,x,y)
4378 index_x = floor((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4379 index_y = floor((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4380 lnx = nx-1
4381 lny = ny-1
4382ENDIF
4383
4384IF (extrap) THEN ! trim indices outside grid for extrapolation
4385 index_x = max(index_x, 1)
4386 index_y = max(index_y, 1)
4387 index_x = min(index_x, lnx)
4388 index_y = min(index_y, lny)
4389ELSE ! nullify indices outside grid
4390 WHERE(index_x < 1 .OR. index_x > lnx .OR. index_y < 1 .OR. index_y > lny)
4391 index_x = imiss
4392 index_y = imiss
4393 END WHERE
4394ENDIF
4395
4396END SUBROUTINE basic_find_index
4397
Module for defining transformations between rectangular georeferenced grids and between grids and spa... Definition grid_transform_class.F90:404 |