Soil Water Balance (SWB2)
Loading...
Searching...
No Matches
grid.F90
Go to the documentation of this file.
1!> @file
2!> Contains a single module, grid, which
3!> provides support for gridded ASCII data file and data structure operations
4!> @ingroup grid
5
6!> Provides support for input and output of gridded ASCII data,
7!> as well as for creation and destruction of grid data structures (defined types).
8module grid
9
10 use iso_c_binding
11 use constants_and_conversions
12 use exceptions
13 use fstring
14 use logfiles
15
16 implicit none
17
18 private
19
20 integer (c_int), public, parameter :: grid_nodata_int = -9999_c_int
21 real (c_float), public, parameter :: grid_nodata_real = -9999._c_float
22 real (c_double), public, parameter :: grid_nodata_double = -9999._c_double
23
24 ! @TODO: why redefine these when types have been defined in 'constants_and_conversions'?
25 integer (c_int), public, parameter :: grid_datatype_int = datatype_int
26 integer (c_int), public, parameter :: grid_datatype_real = datatype_real
27 integer (c_int), public, parameter :: grid_datatype_double = datatype_double
28
29 integer (c_int), public, parameter :: output_surfer = 0
30 integer (c_int), public, parameter :: output_arc = 1
31
32 integer (c_int), public, parameter :: grid_active_cell = 1
33 integer(c_int), parameter :: nc_fill_int = grid_nodata_int
34 real(c_float), parameter :: nc_fill_float = grid_nodata_real
35 real(c_double), parameter :: nc_fill_double = grid_nodata_double
36
37 !> interface to C code that provides a simplified entry point to PROJ4
38 !> capabilities: it has been modified so that all C pointers are kept within the
39 !> C code; no pointers are returned to fortran
40 public :: pj_init_and_transform
41 interface
42 function pj_init_and_transform(from_projection, to_projection, caller_name, &
43 caller_linenum, point_count, x, y) bind(c,name='pj_init_and_transform')
44 import
45 character(c_char) :: from_projection(*)
46 character(c_char) :: to_projection(*)
47 character(c_char) :: caller_name(*)
48 integer(c_int),value :: caller_linenum
49 integer(c_long),value :: point_count
50 real(c_double) :: x(*)
51 real(c_double) :: y(*)
52 integer(c_int) :: pj_init_and_transform
53 end function pj_init_and_transform
54 end interface
55
56
57 type, public :: general_grid_t
58 integer (c_int) :: inx ! Number of cells in the x-direction
59 integer (c_int) :: iny ! Number of cells in the y-direction
60 integer (c_int) :: inumgridcells ! Total number of grid cells
61 integer (c_int) :: idatatype ! Data type contained in the grid (integer, real, SWB cell)
62 character (len=:), allocatable :: sproj4_string ! proj4 string defining coordinate system of grid
63 character (len=:), allocatable :: sfilename ! original file name that the data was read from
64 real (c_double) :: rgridcellsize ! size of one side of a grid cell
65 integer (c_int) :: ilengthunits= -99999 ! length units code
66 real (c_double) :: rx0, rx1 ! World-coordinate range in X
67 real (c_double) :: ry0, ry1 ! World-coordinate range in Y
68
69 integer (c_int), dimension(:,:), allocatable :: idata ! Integer data
70 real (c_float), dimension(:,:), allocatable :: rdata ! Real data
71 real (c_float), dimension(:,:), allocatable :: fdata ! Float data
72 real (c_double), dimension(:,:), allocatable :: dpdata ! Douple-precision data
73 real (c_double), dimension(:,:), allocatable :: rx ! x coordinate associated with data
74 real (c_double), dimension(:,:), allocatable :: ry ! y coordinate associated with data
75 integer (c_int), dimension(:,:), allocatable :: imask ! Mask for processing results
76
77 integer (c_int) :: inodatavalue = nc_fill_int
78 real (c_float) :: rnodatavalue = nc_fill_float
79 real (c_double) :: dpnodatavalue = nc_fill_double
80! type (T_CELL), dimension(:,:), pointer :: Cells ! T_CELL objects
81 end type general_grid_t
82
83
84 type, public :: grid_bounds_t
85 real (c_double) :: rxll, ryll
86 real (c_double) :: rxul, ryul
87 real (c_double) :: rxlr, rylr
88 real (c_double) :: rxur, ryur
89 character (len=256) :: sproj4_string
90 end type grid_bounds_t
91
92
93 type error_message_t
94 character(len=45) :: cerrormessagetext
95 integer :: ierrornumber
96 end type error_message_t
97
98 ! error messages adapted from file "pj_strerrno.c" in PROJ.4.8.0
99 type(error_message_t), dimension(49) :: terrormessage = (/ &
100 error_message_t("no arguments in initialization list ", -1), &
101 error_message_t("no options found in 'init' file ", -2), &
102 error_message_t("no colon in init= string ", -3), &
103 error_message_t("projection not named ", -4), &
104 error_message_t("unknown projection id ", -5), &
105 error_message_t("effective eccentricity = 1. ", -6), &
106 error_message_t("unknown unit conversion id ", -7), &
107 error_message_t("invalid boolean param argument ", -8), &
108 error_message_t("unknown elliptical parameter name ", -9), &
109 error_message_t("reciprocal flattening (1/f) = 0 ", -10), &
110 error_message_t("|radius reference latitude| > 90 ", -11), &
111 error_message_t("squared eccentricity < 0 ", -12), &
112 error_message_t("major axis or radius = 0 or not given ", -13), &
113 error_message_t("latitude or longitude exceeded limits ", -14), &
114 error_message_t("invalid x or y ", -15), &
115 error_message_t("improperly formed DMS value ", -16), &
116 error_message_t("non-convergent inverse meridional dist ", -17), &
117 error_message_t("non-convergent inverse phi2 ", -18), &
118 error_message_t("acos/asin: |arg| >1.+1e-14 ", -19), &
119 error_message_t("tolerance condition error ", -20), &
120 error_message_t("conic lat_1 = -lat_2 ", -21), &
121 error_message_t("lat_1 >= 90 ", -22), &
122 error_message_t("lat_1 = 0 ", -23), &
123 error_message_t("lat_ts >= 90 ", -24), &
124 error_message_t("no distance between control points ", -25), &
125 error_message_t("projection not selected to be rotated ", -26), &
126 error_message_t("W <= 0 or M <= 0 ", -27), &
127 error_message_t("lsat not in 1-5 range ", -28), &
128 error_message_t("path not in range ", -29), &
129 error_message_t("h <= 0 ", -30), &
130 error_message_t("k <= 0 ", -31), &
131 error_message_t("lat_0 = 0 or 90 or alpha = 90 ", -32), &
132 error_message_t("lat_1=lat_2 or lat_1=0 or lat_2=90 ", -33), &
133 error_message_t("elliptical usage required ", -34), &
134 error_message_t("invalid UTM zone number ", -35), &
135 error_message_t("arg(s) out of range for Tcheby eval ", -36), &
136 error_message_t("failed to find projection to be rotated ", -37), &
137 error_message_t("failed to load datum shift file ", -38), &
138 error_message_t("both n ), & m must be spec'd and > 0 ", -39), &
139 error_message_t("n <= 0, n > 1 or not specified ", -40), &
140 error_message_t("lat_1 or lat_2 not specified ", -41), &
141 error_message_t("|lat_1| == |lat_2| ", -42), &
142 error_message_t("lat_0 is pi/2 from mean lat ", -43), &
143 error_message_t("unparseable coordinate system definition ", -44), &
144 error_message_t("geocentric transformation missing z or ellps ", -45), &
145 error_message_t("unknown prime meridian conversion id ", -46), &
146 error_message_t("illegal axis orientation combination ", -47), &
147 error_message_t("point not within available datum shift grids ", -48), &
148 error_message_t("invalid sweep axis, choose x or y ", -49) /)
149
150 public :: grid_gridtogrid_int, grid_gridtogrid_sgl, grid_create, grid_destroy, grid_read, &
151 grid_checkforproj4error, grid_completelycover, &
152 grid_readexisting, grid_dumpgridextent, &
153 grid_conform, grid_transform, grid_interpolate, grid_populatexy, &
154 grid_getgridx, grid_getgridy, &
155 grid_set_nodata_value, &
156 grid_createcomplete, grid_createsimple, grid_writegrid, &
157 grid_writearcgrid, grid_writesurfergrid, grid_set_output_directory_name, &
158 grid_getgridrownum, grid_getgridcolnum, &
159 grid_rowcolfallsinsidegrid, grid_coordinatesfallinsidegrid
160
161 interface grid_create
162 module procedure grid_createsimple
163 module procedure grid_createcomplete
164 end interface grid_create
165
166 ! global parameters for use with the majority filter
167 integer (c_int), parameter :: four_cells = 1
168 integer (c_int), parameter :: eight_cells = 2
169
170 integer (c_int), parameter, private :: column = 1
171 integer (c_int), parameter, private :: row = 2
172
173 !> @todo change these global (module) variables to local variables
174 integer (c_int) :: lu_temp, lu_grid
175
176 character (len=256) :: output_grid_directory_name = ""
177
178contains
179
180
181 subroutine grid_set_output_directory_name( sDirName )
182
183 character (len=*), intent(in) :: sdirname
184
185 output_grid_directory_name = trim(sdirname)
186
187 call logs%write("ASCII grids will be written to subdirectory " &
188 //dquote( output_grid_directory_name ), iloglevel=log_all, &
189 ilinesbefore=1, lecho=true )
190
191 end subroutine grid_set_output_directory_name
192
193
194!--------------------------------------------------------------------------
195
196!> Creates a grid of a specified type.
197!>
198!> Creates a grid pointer object and allocates memory for the data
199!> associated with the grid (REAL, INTEGER, or T_CELL).
200!
201!> @param iNX Number of grid cells in the x direction
202!> @param iNY Number of grid cells in the y direction
203!> @param rX0 X coordinate for the lower left corner of the grid
204!> @param rY0 Y coordinate for the lower left corner of the grid
205!> @param rX1 X coordinate for the upper right corner of the grid
206!> @param rY1 Y coordinate for the upper right corner of the grid
207!> @param iDataType Integer value corresponding to the type of data contained in the grid
208!>
209!> @return pGrd Pointer to a grid object
210function grid_createcomplete ( iNX, iNY, rX0, rY0, rX1, rY1, iDataType ) result ( pGrd )
211
212 ! ARGUMENTS
213 integer (c_int), intent(in) :: inx, iny ! Grid dimensions
214 real (c_double), intent(in) :: rx0, ry0 ! Lower-left corner (world coords)
215 real (c_double), intent(in) :: rx1, ry1 ! Upper-right corner (world coords)
216 integer (c_int), intent(in) :: idatatype ! Data type (DATATYPE_INT, etc.)
217 ! RETURN VALUE
218 type (general_grid_t), pointer :: pgrd
219 ! LOCALS
220 integer (c_int) :: istat
221
222 allocate ( pgrd, stat=istat )
223 call assert ( istat == 0, &
224 "Could not allocate pointer to T_GRID object", __file__,__line__ )
225
226 if (inx <= 0 .or. iny <= 0) then
227 call logs%write("Illegal grid dimensions: ")
228 call logs%write("iNX: "//ascharacter(inx) )
229 call logs%write("iNY: "//ascharacter(iny) )
230 call logs%write("rX0: "//ascharacter(rx0) )
231 call logs%write("rY0: "//ascharacter(ry0) )
232 call logs%write("rX1: "//ascharacter(rx1) )
233 call logs%write("rY1: "//ascharacter(ry1) )
234 call assert ( false, &
235 "INTERNAL PROGRAMMING ERROR? - Illegal grid dimensions specified", __file__,__line__)
236 endif
237
238 select case (idatatype)
239 case ( grid_datatype_int )
240
241 allocate ( pgrd%iData( inx, iny ), stat=istat )
242 call assert (istat == 0, &
243 "Could not allocate integer data", &
244 __file__,__line__)
245 pgrd%iData = 0
246
247 case ( grid_datatype_real )
248 allocate ( pgrd%rData( inx, iny ), stat=istat )
249 call assert (istat == 0, &
250 "Could not allocate real data", &
251 __file__,__line__)
252 pgrd%rData = rzero
253
254 case ( grid_datatype_double )
255 allocate ( pgrd%dpData( inx, iny ), stat=istat )
256 call assert (istat == 0, &
257 "Could not allocate double-precision data", &
258 __file__,__line__)
259 pgrd%dpData = 0.0_c_double
260
261 case default
262 call assert ( false, 'Internal error -- illegal grid data type' )
263 end select
264
265 pgrd%iDataType = idatatype
266 pgrd%iNX = inx
267 pgrd%iNY = iny
268 pgrd%rX0 = rx0
269 pgrd%rX1 = rx1
270 pgrd%rY0 = ry0
271 pgrd%rY1 = ry1
272 pgrd%rGridCellSize = (pgrd%rX1 - pgrd%rX0) / real(pgrd%iNX, c_double)
273 pgrd%iNumGridCells = inx * iny
274
275 allocate(pgrd%iMask(inx, iny))
276 pgrd%iMask = grid_active_cell
277
278end function grid_createcomplete
279
280function grid_createsimple ( iNX, iNY, rX0, rY0, rGridCellSize, iDataType ) result ( pGrd )
281 !! Creates a new iNX-by-iNY T_GRID of data type iDataType, over the extent
282 !! (rX0,rY0)-(rX1,rY1), and returns a pointer.
283 ! ARGUMENTS
284 integer (c_int), intent(in) :: inx, iny ! Grid dimensions
285 real (c_double), intent(in) :: rx0, ry0 ! Lower-left corner (world coords)
286 real (c_double), intent(in) :: rgridcellsize
287 integer (c_int), intent(in) :: idatatype ! Data type (DATATYPE_INT, etc.)
288
289 ! RETURN VALUE
290 type (general_grid_t), pointer :: pgrd
291 ! LOCALS
292 integer (c_int) :: istat
293
294 allocate ( pgrd, stat=istat )
295 call assert ( istat == 0, &
296 "Could not allocate pointer to T_GRID object", __file__,__line__ )
297 call assert ( inx > 0 .and. iny > 0, &
298 "Illegal grid dimensions specified: NX="//ascharacter(inx) &
299 //" NY="//ascharacter(iny), __file__,__line__)
300
301 select case (idatatype)
302 case ( grid_datatype_int )
303
304 allocate ( pgrd%iData( inx, iny ), stat=istat )
305 call assert (istat == 0, &
306 "Could not allocate integer data", &
307 __file__,__line__)
308 pgrd%iData = pgrd%iNoDataValue
309
310 case ( grid_datatype_real )
311 allocate ( pgrd%rData( inx, iny ), stat=istat )
312 call assert (istat == 0, &
313 "Could not allocate real data", &
314 __file__,__line__)
315 pgrd%rData = pgrd%rNoDataValue
316
317 case ( grid_datatype_double )
318 allocate ( pgrd%dpData( inx, iny ), stat=istat )
319 call assert (istat == 0, &
320 "Could not allocate double-precision data", &
321 __file__,__line__)
322 pgrd%dpData = pgrd%dpNoDataValue
323
324 case default
325 call assert ( false, 'Internal error -- illegal grid data type' )
326 end select
327
328 pgrd%iDataType = idatatype
329 pgrd%iNX = inx
330 pgrd%iNY = iny
331 pgrd%rX0 = rx0
332 pgrd%rX1 = rx0 + real(inx, c_double) * rgridcellsize
333 pgrd%rY0 = ry0
334 pgrd%rY1 = ry0 + real(iny, c_double) * rgridcellsize
335 pgrd%rGridCellSize = rgridcellsize
336 pgrd%iNumGridCells = inx * iny
337
338 allocate(pgrd%iMask(inx, iny))
339 pgrd%iMask = grid_active_cell
340
341end function grid_createsimple
342
343!!***
344
345!--------------------------------------------------------------------------
346!!****f* grid/grid_Destroy
347! NAME
348! grid_Destroy - Destroys a grid of a specified type.
349!
350! SYNOPSIS
351! Destroys a grid pointer object and deallocates memory for the data
352! associated with the grid (REAL, INTEGER, or T_CELL).
353!
354! INPUTS
355! pGrd - Pointer to a grid object
356!
357! OUTPUTS
358! None
359!
360! NOTES
361! Code refers to parameters that are set within types.f95.
362!
363! SOURCE
364
365subroutine grid_destroy ( pGrd )
366 !! Destroys the data in the GENERAL_GRID_T pGrd
367 ! ARGUMENTS
368 type ( general_grid_t ), pointer :: pgrd
369 ! LOCALS
370 integer (c_int) :: istat
371
372 if(associated(pgrd) )then
373
374 if ( pgrd%iDataType == grid_datatype_int ) then
375 deallocate ( pgrd%iData, stat=istat )
376 call assert ( istat == 0, "Failed to deallocate integer grid" )
377 else if ( pgrd%iDataType == grid_datatype_real ) then
378 deallocate ( pgrd%rData, stat=istat )
379 call assert ( istat == 0, "Failed to deallocate real grid" )
380 else if ( pgrd%iDataType == grid_datatype_double ) then
381 deallocate ( pgrd%dpData, stat=istat )
382 call assert ( istat == 0, "Failed to deallocate double-precision grid" )
383
384! else if ( pGrd%iDataType == GRID_DATATYPE_CELL_GRID ) then
385! deallocate ( pGrd%Cells, stat=iStat )
386! call assert ( iStat == 0, "Failed to deallocate cell grid" )
387! else if ( pGrd%iDataType == GRID_DATATYPE_ALL ) then
388! deallocate ( pGrd%iData, stat=iStat )
389! call assert ( iStat == 0, "Failed to deallocate integer grid" )
390! deallocate ( pGrd%rData, stat=iStat )
391! call assert ( iStat == 0, "Failed to deallocate real grid" )
392! deallocate ( pGrd%Cells, stat=iStat )
393! call assert ( iStat == 0, "Failed to deallocate cell grid" )
394 else
395 call assert ( false, "Internal error -- unknown grid type", &
396 __file__, __line__)
397 end if
398
399 if( allocated(pgrd%rX) ) then
400 deallocate( pgrd%rX, stat=istat)
401 call assert ( istat == 0, "Failed to deallocate X-coordinate data structure associated with grid", &
402 __file__, __line__ )
403 endif
404
405 if( allocated(pgrd%rY) ) then
406 deallocate( pgrd%rY, stat=istat)
407 call assert ( istat == 0, "Failed to deallocate Y-coordinate data structure associated with grid", &
408 __file__, __line__ )
409 endif
410
411 endif
412
413 pgrd => null()
414
415end subroutine grid_destroy
416!!***
417
418!--------------------------------------------------------------------------
419!!****f* grid/grid_Read
420! NAME
421! grid_Read - Reads a grid of a specified type.
422!
423! SYNOPSIS
424! Reads a grid of the
425! INPUTS
426! sFileName - Character string containing the name of the file to be read
427! sFileType - Character string indicating the type of file to be read
428! iDataType - Integer value corresponding to the type of data contained
429! in the input data file
430!
431! OUTPUTS
432! pGrd - Pointer to a grid object
433!
434! NOTES
435!
436! Current legal file types are
437! "ARC_GRID" ESRI ARC Grid (ASCII)
438! "SURFER" Golden Software SURFER grid (ASCII)
439!
440! Valid data types are (see module 'types'):
441! DATATYPE_INT
442! DATATYPE_REAL
443!
444! SOURCE
445
446function grid_read ( sFilename, sFileType, iDataType ) result ( pGrd )
447
448 ! ARGUMENTS
449 character (len=*), intent(in) :: sfilename ! Name of the grid input file
450 character (len=*), intent(in) :: sfiletype ! File type (see above)
451 integer (c_int), intent(in) :: idatatype ! T_GRID_INT or T_GRID_REAL
452 ! RETURN VALUE
453 type (general_grid_t), pointer :: pgrd
454
455 pgrd => null()
456
457 if ( trim(sfiletype) == "ARC_GRID" ) then
458 pgrd => grid_readarcgrid_fn( sfilename, idatatype )
459 else if ( trim(sfiletype) == "SURFER" ) then
460 pgrd => grid_readsurfergrid_fn( sfilename, idatatype )
461 else
462 call assert( false, "Illegal grid file type requested" )
463 end if
464
465 pgrd%sFilename = trim(sfilename)
466
467end function grid_read
468
469!--------------------------------------------------------------------------
470
471subroutine grid_readexisting ( sFileName, sFileType, pGrd )
472
473 ! ARGUMENTS
474 character (len=*), intent(in) :: sfilename ! Name of the grid input file
475 character (len=*), intent(in) :: sfiletype ! File type (see above)
476 type (general_grid_t), pointer :: pgrd
477
478 if ( trim(sfiletype) == "ARC_GRID" ) then
479 call grid_readarcgrid_sub( sfilename, pgrd )
480 else if ( trim(sfiletype) == "SURFER" ) then
481 call grid_readsurfergrid_sub( sfilename, pgrd )
482 else
483 call assert( false, "Illegal grid file type requested" )
484 end if
485
486 pgrd%sFilename = trim(sfilename)
487
488end subroutine grid_readexisting
489
490!!***
491!--------------------------------------------------------------------------
492!!****f* grid/grid_ReadArcGrid_fn
493! NAME
494! grid_ReadArcGrid_fn - Reads an ARC ASCII grid of a specified type.
495!
496! SYNOPSIS
497! Reads an ARC ASCII grid of specified data type
498!
499! INPUTS
500! sFileName - Character string containing the name of the file to be read
501! iDataType - Integer value corresponding to the type of data contained
502! in the input data file
503!
504! OUTPUTS
505! pGrd - Pointer to a grid object
506!
507! NOTES
508!
509! Valid data types are (see module 'types'):
510! DATATYPE_INT
511! DATATYPE_REAL
512!
513! SOURCE
514
515function grid_readarcgrid_fn ( sFileName, iDataType ) result ( pGrd )
516
517 ! ARGUMENTS
518 character (len=*), intent(in) :: sfilename ! Name of the grid input file
519 integer (c_int), intent(in) :: idatatype ! T_GRID_INT or T_GRID_REAL
520 ! RETURN VALUE
521 type (general_grid_t), pointer :: pgrd
522 ! LOCALS
523 character (len=256) :: sinputrecord ! Text record for input
524 character (len=256) :: sdirective ! Directive for input
525 character (len=256) :: sargument ! Argument for keyword directives
526 character (len=256) :: snodatavalue ! String to hold nodata value
527 integer (c_int) :: istat ! For "iostat="
528 integer (c_int) :: inx,iny ! Grid dimensions
529 integer (c_int) :: ihdrrecs ! Number of records in header
530 integer (c_int) :: icol,irow ! Loop indices for grid reading
531 real (c_double) :: rx0,rx1 ! Limits in X
532 real (c_double) :: ry0,ry1 ! Limits in Y
533 real (c_double) :: rcellsize ! Cell size
534 logical (c_bool) :: lxllcenter, lyllcenter ! Flags XLLCENTER / XLLCORNER
535 logical :: lfileexists
536 logical :: lisopen
537
538 ! Pre-scan for the number of header records and read the header
539 inquire(file=trim(sfilename), exist=lfileexists)
540 call assert( lfileexists, "The Arc ASCII grid file "//dquote(sfilename)// &
541 " could not be found.",__file__,__line__)
542
543 inquire(unit=lu_grid, opened=lisopen )
544 if (lisopen ) close( unit=lu_grid )
545
546 open ( newunit=lu_grid, iostat=istat, file=trim(sfilename) )
547 call assert( istat == 0, &
548 "Could not open input file " // trim(sfilename) )
549
550 ihdrrecs = 0
551 lxllcenter = false
552 lyllcenter = false
553 snodatavalue = ""
554
555 do
556 read ( unit=lu_grid, fmt="(a256)", iostat=istat ) sinputrecord
557 call assert ( istat == 0, &
558 "Could not read input record - file:"//trim(sfilename) )
559 call chomp ( sinputrecord, sdirective )
560 call touppercase ( sdirective )
561 call chomp ( sinputrecord, sargument )
562 if ( sdirective == "NCOLS" ) then
563 read ( unit=sargument, fmt=*, iostat=istat ) inx
564 call assert ( istat == 0, "Could not read NCOLS" )
565 ihdrrecs = ihdrrecs + 1
566 else if ( sdirective == "NROWS" ) then
567 read ( unit=sargument, fmt=*, iostat=istat ) iny
568 call assert ( istat == 0, "Could not read NROWS" )
569 ihdrrecs = ihdrrecs + 1
570 else if ( sdirective == "XLLCENTER" ) then
571 read ( unit=sargument, fmt=*, iostat=istat ) rx0
572 call assert ( istat == 0, "Could not read XLLCENTER" )
573 lxllcenter = true
574 ihdrrecs = ihdrrecs + 1
575 else if ( sdirective == "YLLCENTER" ) then
576 read ( unit=sargument, fmt=*, iostat=istat ) ry0
577 call assert ( istat == 0, "Could not read YLLCENTER" )
578 lxllcenter = true
579 ihdrrecs = ihdrrecs + 1
580 else if ( sdirective == "XLLCORNER" ) then
581 read ( unit=sargument, fmt=*, iostat=istat ) rx0
582 call assert ( istat == 0, "Could not read XLLCORNER" )
583 ihdrrecs = ihdrrecs + 1
584 else if ( sdirective == "YLLCORNER" ) then
585 read ( unit=sargument, fmt=*, iostat=istat ) ry0
586 call assert ( istat == 0, "Could not read YLLCORNER" )
587 ihdrrecs = ihdrrecs + 1
588 else if ( sdirective == "CELLSIZE" ) then
589 read ( unit=sargument, fmt=*, iostat=istat ) rcellsize
590 call assert ( istat == 0, "Could not read CELLSIZE" )
591 ihdrrecs = ihdrrecs + 1
592 else if ( sdirective == "NODATA_VALUE" ) then
593 snodatavalue = trim(sargument)
594 ihdrrecs = ihdrrecs + 1
595 else
596 ! Found the data -- construct the grid
597 if ( lxllcenter ) rx0 = rx0 - 0.5_c_double*rcellsize
598 if ( lyllcenter ) ry0 = ry0 - 0.5_c_double*rcellsize
599 rx1 = rx0 + real(inx, c_double) * rcellsize
600 ry1 = ry0 + real(iny, c_double) * rcellsize
601
602 pgrd => grid_create( inx, iny, rx0, ry0, rx1, ry1, idatatype )
603
604 pgrd%rGridCellSize = rcellsize
605 ! Go back to the top, skip the header...
606 rewind( unit=lu_grid, iostat=istat )
607 call assert ( istat == 0, "Failed to rewind grid file" )
608 do icol=1,ihdrrecs
609 read ( unit=lu_grid, fmt="(a256)", iostat=istat ) sinputrecord
610 call assert ( istat == 0, &
611 "Could not read input record - file: "//trim(sfilename))
612 end do
613 ! ... and read the data.
614 select case ( idatatype )
615
616 case ( datatype_int )
617
618 do irow=1,pgrd%iNY
619 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%iData(:,irow)
620 call assert ( istat == 0, &
621 "Failed to read integer grid data - file: " &
622 //trim(sfilename)//" row num: "//trim( ascharacter(irow)), &
623 __file__,__line__ )
624 end do
625 if(len_trim(snodatavalue) > 0) then
626 read(unit=snodatavalue, fmt=*, iostat=istat) pgrd%iNoDataValue
627 call assert ( istat == 0, &
628 "Failed to read NODATA value in grid data - file: " &
629 //dquote(sfilename), __file__,__line__ )
630 endif
631
632 case ( datatype_real )
633
634 do irow=1,pgrd%iNY
635 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%rData(:,irow)
636 call assert ( istat == 0, &
637 "Failed to read real grid data - file: " &
638 //trim(sfilename)//" row num: "//trim( ascharacter(irow)), &
639 __file__,__line__ )
640 end do
641 if(len_trim(snodatavalue) > 0) then
642 read(unit=snodatavalue, fmt=*, iostat=istat) pgrd%rNoDataValue
643 call assert ( istat == 0, &
644 "Failed to read NODATA value in grid data - file: " &
645 //trim(sfilename), __file__,__line__ )
646 endif
647
648 case ( datatype_double )
649
650 do irow=1,pgrd%iNY
651 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%dpData(:,irow)
652 call assert ( istat == 0, &
653 "Failed to read double-precision grid data - file: " &
654 //trim(sfilename)//" row num: "//trim( ascharacter(irow)), &
655 __file__,__line__ )
656 end do
657 if(len_trim(snodatavalue) > 0) then
658 read(unit=snodatavalue, fmt=*, iostat=istat) pgrd%dpNoDataValue
659 call assert ( istat == 0, &
660 "Failed to read NODATA value in grid data - file: " &
661 //trim(sfilename), __file__,__line__ )
662 endif
663
664 case default
665
666 call assert ( false, &
667 "Internal error -- illegal ARC GRID data type", &
668 __file__,__line__)
669
670 end select
671 exit
672 end if
673 end do
674
675 if ( idatatype == datatype_int ) call grid_checkintegergridvalues(pgrd, sfilename)
676
677 open ( unit=lu_grid, iostat=istat )
678 call assert ( istat == 0, "Failed to close grid file" )
679
680end function grid_readarcgrid_fn
681
682!--------------------------------------------------------------------------
683
684subroutine grid_readarcgrid_sub ( sFileName, pGrd )
685
686 ! ARGUMENTS
687 character (len=*), intent(in) :: sFileName ! Name of the grid input file
688 type (GENERAL_GRID_T), pointer :: pGrd
689
690 ! LOCALS
691 character (len=256) :: sInputRecord ! Text record for input
692 character (len=256) :: sDirective ! Directive for input
693 character (len=256) :: sArgument ! Argument for keyword directives
694 character (len=256) :: sNoDataValue ! String to hold nodata value
695 integer (c_int) :: iStat ! For "iostat="
696 integer (c_int) :: iNX,iNY ! Grid dimensions
697 integer (c_int) :: iHdrRecs ! Number of records in header
698 integer (c_int) :: iCol,iRow ! Loop indices for grid reading
699 real (c_double) :: rX0 ! lower left-hand corner X
700 real (c_double) :: rY0 ! lower left-hand corner Y
701 real (c_double) :: rCellSize ! Cell size
702 logical (c_bool) :: lXLLCenter, lYLLCenter ! Flags XLLCENTER / XLLCORNER
703 logical :: lFileExists
704 logical :: lIsOpen
705
706 ! Pre-scan for the number of header records and read the header
707
708 inquire(file=trim(sfilename), exist=lfileexists, opened=lisopen )
709 call assert( lfileexists, "The Arc ASCII grid file "//dquote(sfilename)// &
710 " could not be found.",__file__,__line__)
711
712 inquire(unit=lu_grid, opened=lisopen )
713 if (lisopen ) close( unit=lu_grid )
714
715 open ( newunit=lu_grid, iostat=istat, file=trim(sfilename) )
716 call assert( istat == 0, &
717 "Could not open input file " // trim(sfilename) )
718
719 ihdrrecs = 0
720 lxllcenter = false
721 lyllcenter = false
722 snodatavalue = ""
723
724 do
725 read ( unit=lu_grid, fmt="(a256)", iostat=istat ) sinputrecord
726 call assert ( istat == 0, &
727 "Could not read input record - file:"//trim(sfilename) )
728 call chomp ( sinputrecord, sdirective )
729 call touppercase ( sdirective )
730 call chomp ( sinputrecord, sargument )
731 if ( sdirective == "NCOLS" ) then
732 read ( unit=sargument, fmt=*, iostat=istat ) inx
733 call assert ( istat == 0, "Could not read NCOLS" )
734 ihdrrecs = ihdrrecs + 1
735 else if ( sdirective == "NROWS" ) then
736 read ( unit=sargument, fmt=*, iostat=istat ) iny
737 call assert ( istat == 0, "Could not read NROWS" )
738 ihdrrecs = ihdrrecs + 1
739 else if ( sdirective == "XLLCENTER" ) then
740 read ( unit=sargument, fmt=*, iostat=istat ) rx0
741 call assert ( istat == 0, "Could not read XLLCENTER" )
742 lxllcenter = true
743 ihdrrecs = ihdrrecs + 1
744 else if ( sdirective == "YLLCENTER" ) then
745 read ( unit=sargument, fmt=*, iostat=istat ) ry0
746 call assert ( istat == 0, "Could not read YLLCENTER" )
747 lxllcenter = true
748 ihdrrecs = ihdrrecs + 1
749 else if ( sdirective == "XLLCORNER" ) then
750 read ( unit=sargument, fmt=*, iostat=istat ) rx0
751 call assert ( istat == 0, "Could not read XLLCORNER" )
752 ihdrrecs = ihdrrecs + 1
753 else if ( sdirective == "YLLCORNER" ) then
754 read ( unit=sargument, fmt=*, iostat=istat ) ry0
755 call assert ( istat == 0, "Could not read YLLCORNER" )
756 ihdrrecs = ihdrrecs + 1
757 else if ( sdirective == "CELLSIZE" ) then
758 read ( unit=sargument, fmt=*, iostat=istat ) rcellsize
759 call assert ( istat == 0, "Could not read CELLSIZE" )
760 ihdrrecs = ihdrrecs + 1
761 else if ( sdirective == "NODATA_VALUE" ) then
762 snodatavalue = trim(sargument)
763 ihdrrecs = ihdrrecs + 1
764 else
765 ! Found the data -- construct the grid
766
767 ! Go back to the top, skip the header...
768 rewind( unit=lu_grid, iostat=istat )
769 call assert ( istat == 0, "Failed to rewind grid file" )
770 do icol=1,ihdrrecs
771 read ( unit=lu_grid, fmt="(a256)", iostat=istat ) sinputrecord
772 call assert ( istat == 0, &
773 "Could not read input record - file: "//trim(sfilename))
774 end do
775 ! ... and read the data.
776 select case ( pgrd%iDataType )
777
778 case ( datatype_int )
779 do irow=1,pgrd%iNY
780 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%iData(:,irow)
781 call assert ( istat == 0, &
782 "Failed to read integer grid data - file: " &
783 //trim(sfilename)//"; row num: "//trim( ascharacter(irow)) &
784 //"; line num: "//trim( ascharacter(irow + ihdrrecs)) &
785 //'; error code: '//as_character(istat), &
786 __file__,__line__ )
787 end do
788 if(len_trim(snodatavalue) > 0) then
789 read(unit=snodatavalue, fmt=*, iostat=istat) pgrd%iNoDataValue
790 call assert ( istat == 0, &
791 "Failed to read NODATA value in grid data - file: " &
792 //trim(sfilename), __file__,__line__ )
793 endif
794
795 case ( datatype_real )
796
797 do irow=1,pgrd%iNY
798 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%rData(:,irow)
799
800 call assert ( istat == 0, &
801 "Failed to read real grid data - file: " &
802 //trim(sfilename)//"; row num: "//trim( ascharacter(irow)) &
803 //"; line num: "//trim( ascharacter(irow + ihdrrecs)) &
804 //'; error code: '//as_character(istat), &
805 __file__,__line__ )
806 end do
807 if(len_trim(snodatavalue) > 0) then
808 read(unit=snodatavalue, fmt=*, iostat=istat) pgrd%rNoDataValue
809 call assert ( istat == 0, &
810 "Failed to read NODATA value in grid data - file: " &
811 //trim(sfilename), __file__,__line__ )
812 endif
813
814 case ( datatype_double )
815 do irow=1,pgrd%iNY
816 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%dpData(:,irow)
817 call assert ( istat == 0, &
818 "Failed to read double-precision grid data - file: " &
819 //trim(sfilename)//"; row num: "//trim( ascharacter(irow)) &
820 //"; line num: "//trim( ascharacter(irow + ihdrrecs)) &
821 //'; error code: '//as_character(istat), &
822 __file__,__line__ )
823 end do
824 if(len_trim(snodatavalue) > 0) then
825 read(unit=snodatavalue, fmt=*, iostat=istat) pgrd%dpNoDataValue
826 call assert ( istat == 0, &
827 "Failed to read NODATA value in grid data - file: " &
828 //trim(sfilename), __file__,__line__ )
829 endif
830
831 case default
832 call assert ( false, &
833 "Internal error -- illegal ARC GRID data type", &
834 __file__,__line__)
835 end select
836 exit
837 end if
838 end do
839
840 if ( pgrd%iDataType == datatype_int ) call grid_checkintegergridvalues(pgrd, sfilename)
841
842 open ( unit=lu_grid, iostat=istat )
843 call assert ( istat == 0, "Failed to close grid file" )
844
845end subroutine grid_readarcgrid_sub
846
847!!***
848!--------------------------------------------------------------------------
849!!****f* grid/grid_ReadSurferGrid_fn
850! NAME
851! grid_ReadSurferGrid_fn - Reads an Surfer ASCII grid of a specified type.
852!
853! SYNOPSIS
854! Reads an Surfer ASCII grid of specified data type
855!
856! INPUTS
857! sFileName - Character string containing the name of the file to be read
858! iDataType - Integer value corresponding to the type of data contained
859! in the input data file
860!
861! OUTPUTS
862! pGrd - Pointer to a grid object
863!
864! NOTES
865!
866! Valid data types are (see module 'types'):
867! DATATYPE_INT
868! DATATYPE_REAL
869!
870! SOURCE
871
872function grid_readsurfergrid_fn ( sFileName, iDataType ) result ( pGrd )
873
874 ! ARGUMENTS
875 character (len=*), intent(in) :: sfilename ! Name of the grid input file
876 integer (c_int), intent(in) :: idatatype ! T_GRID_INT or T_GRID_REAL
877 ! RETURN VALUE
878 type (general_grid_t), pointer :: pgrd
879 ! LOCALS
880 character (len=4) :: ssentinel ! Better be "DSAA" for SURFER!
881 integer (c_int) :: istat ! For "iostat="
882 integer (c_int) :: inx,iny ! Grid dimensions
883 integer (c_int) :: irow ! Loop index for grid reading
884 real (c_double) :: rx0,rx1 ! Limits in X
885 real (c_double) :: ry0,ry1 ! Limits in Y
886 real (c_float) :: rz0,rz1 ! Limits in Z (not used)
887 logical :: lfileexists
888 logical :: lisopen
889
890 inquire(file=trim(sfilename), exist=lfileexists)
891 call assert( lfileexists, "The Surfer ASCII grid file "//dquote(sfilename)// &
892 " could not be found.",__file__,__line__)
893
894 inquire(unit=lu_grid, opened=lisopen )
895 if (lisopen ) close( unit=lu_grid )
896
897 open ( newunit=lu_grid, iostat=istat, file=trim(sfilename) )
898 call assert( istat == 0, &
899 "Could not open input file " // trim(sfilename) )
900
901 read ( unit=lu_grid, fmt=*, iostat=istat ) ssentinel
902 call assert ( istat == 0, &
903 "Could not read first record of SURFER grid" )
904 call assert ( LOGICAL(trim(sSentinel) == "DSAA",c_bool), &
905 "This is not a SURFER grid" )
906 read ( unit=lu_grid, fmt=*, iostat=istat ) inx, iny
907 call assert ( istat == 0, &
908 "Error reading SURFER grid dimensions" )
909 read ( unit=lu_grid, fmt=*, iostat=istat ) rx0, rx1
910 call assert ( istat == 0, &
911 "Error reading SURFER X limits" )
912 read ( unit=lu_grid, fmt=*, iostat=istat ) ry0, ry1
913 call assert ( istat == 0, &
914 "Error reading SURFER y limits" )
915 read ( unit=lu_grid, fmt=*, iostat=istat ) rz0, rz1
916 call assert ( istat == 0, &
917 "Error reading SURFER Z limits" )
918
919 pgrd => grid_create( inx, iny, rx0, ry0, rx1, ry1, idatatype )
920 select case ( idatatype )
921 case ( datatype_int )
922 do irow=1, iny
923 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%iData(:,irow)
924 call assert ( istat == 0, "Failed to read integer grid data" )
925 end do
926 case ( datatype_real )
927 do irow=1, iny
928 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%rData(:,irow)
929 call assert ( istat == 0, "Failed to read real grid data" )
930 end do
931 case ( datatype_double )
932 do irow=1, iny
933 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%dpData(:,irow)
934 call assert ( istat == 0, "Failed to read double-precision grid data" )
935 end do
936 case default
937 call assert ( false, "Internal error -- illegal SURFER grid data type" )
938 end select
939
940 call assert(inx>0, "Must have a non-zero number of columns surfer grid file...")
941 call assert(iny>0, "Must have a non-zero number of rows in a surfer grid file...")
942
943 pgrd%rGridCellSize = (rx1-rx0)/inx
944
945 if ( idatatype == datatype_int ) call grid_checkintegergridvalues(pgrd, sfilename)
946
947 open ( unit=lu_grid, iostat=istat )
948 call assert ( istat == 0, "Failed to close grid file" )
949
950end function grid_readsurfergrid_fn
951
952!------------------------------------------------------------------------------
953
954subroutine grid_readsurfergrid_sub ( sFileName, pGrd )
955
956 ! ARGUMENTS
957 character (len=*), intent(in) :: sFileName ! Name of the grid input file
958 type (GENERAL_GRID_T), pointer :: pGrd
959
960 ! LOCALS
961 character (len=4) :: sSentinel ! Better be "DSAA" for SURFER!
962 integer (c_int) :: iStat ! For "iostat="
963 integer (c_int) :: iNX,iNY ! Grid dimensions
964 integer (c_int) :: iRow ! Loop index for grid reading
965 real (c_double) :: rX0,rX1 ! Limits in X
966 real (c_double) :: rY0,rY1 ! Limits in Y
967 real (c_float) :: rZ0,rZ1 ! Limits in Z (not used)
968 logical :: lFileExists
969 logical :: lIsOpen
970
971 inquire(file=trim(sfilename), exist=lfileexists)
972 call assert( lfileexists, "The Surfer ASCII grid file "//dquote(sfilename)// &
973 " could not be found.",__file__,__line__)
974
975 inquire(unit=lu_grid, opened=lisopen )
976 if (lisopen ) close( unit=lu_grid )
977
978 open ( newunit=lu_grid, iostat=istat, file=trim(sfilename) )
979 call assert( istat == 0, &
980 "Could not open input file " // trim(sfilename) )
981
982 read ( unit=lu_grid, fmt=*, iostat=istat ) ssentinel
983 call assert ( istat == 0, &
984 "Could not read first record of SURFER grid" )
985 call assert ( LOGICAL(trim(sSentinel) == "DSAA",c_bool), &
986 "This is not a SURFER grid" )
987 read ( unit=lu_grid, fmt=*, iostat=istat ) inx, iny
988 call assert ( istat == 0, &
989 "Error reading SURFER grid dimensions" )
990 read ( unit=lu_grid, fmt=*, iostat=istat ) rx0, rx1
991 call assert ( istat == 0, &
992 "Error reading SURFER X limits" )
993 read ( unit=lu_grid, fmt=*, iostat=istat ) ry0, ry1
994 call assert ( istat == 0, &
995 "Error reading SURFER y limits" )
996 read ( unit=lu_grid, fmt=*, iostat=istat ) rz0, rz1
997 call assert ( istat == 0, &
998 "Error reading SURFER Z limits" )
999
1000 select case ( pgrd%iDataType )
1001 case ( datatype_int )
1002 do irow=1, iny
1003 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%iData(:,irow)
1004 call assert ( istat == 0, "Failed to read integer grid data" )
1005 end do
1006 case ( datatype_real )
1007 do irow=1, iny
1008 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%rData(:,irow)
1009 call assert ( istat == 0, "Failed to read real grid data" )
1010 end do
1011 case ( datatype_double )
1012 do irow=1, iny
1013 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%dpData(:,irow)
1014 call assert ( istat == 0, "Failed to read double-precision grid data" )
1015 end do
1016 case default
1017 call assert ( false, "Internal error -- illegal SURFER grid data type" )
1018 end select
1019
1020 call assert(inx > 0,"Must have a non-zero number of grid cells in a surfer grid file...")
1021 call assert(iny > 0,"Must have a non-zero number of grid cells in a surfer grid file...")
1022
1023 if ( pgrd%iDataType == datatype_int ) call grid_checkintegergridvalues(pgrd, sfilename)
1024
1025 open ( unit=lu_grid, iostat=istat )
1026 call assert ( istat == 0, "Failed to close grid file" )
1027
1028end subroutine grid_readsurfergrid_sub
1029
1030!------------------------------------------------------------------------------
1031
1032subroutine grid_writegrid(sFilename, pGrd, iOutputFormat)
1033
1034 ! [ ARGUMENTS ]
1035 character (len=*),intent(in) :: sfilename
1036 type (general_grid_t), pointer :: pgrd
1037 integer (c_int) :: ioutputformat
1038
1039 if ( ioutputformat == output_arc ) then
1040
1041 call grid_writearcgrid(trim(output_grid_directory_name)//sfilename, pgrd)
1042
1043 elseif ( ioutputformat == output_surfer ) then
1044
1045 call grid_writesurfergrid(trim(output_grid_directory_name)//sfilename, pgrd)
1046
1047 endif
1048
1049end subroutine grid_writegrid
1050
1051!------------------------------------------------------------------------------
1052
1053subroutine grid_writearcgrid(sFilename, pGrd)
1054
1055 ! [ ARGUMENTS ]
1056 character (len=*),intent(in) :: sfilename
1057 type (general_grid_t), pointer :: pgrd
1058
1059
1060 ! [ LOCALS ]
1061 integer (c_int) :: icol,irow
1062 integer (c_int) :: inumcols, inumrows
1063 integer (c_int) :: istat
1064 character(len=256) :: sbuf
1065
1066 if ( pgrd%iDataType == datatype_int ) then
1067 inumcols = size(pgrd%iData,1)
1068 inumrows = size(pgrd%iData,2)
1069 elseif ( pgrd%iDataType == datatype_real ) then
1070 inumcols = size(pgrd%rData,1)
1071 inumrows = size(pgrd%rData,2)
1072 elseif ( pgrd%iDataType == datatype_double ) then
1073 inumcols = size(pgrd%dpData,1)
1074 inumrows = size(pgrd%dpData,2)
1075 else
1076 call assert(false, "Internal programming error - Unsupported grid type", &
1077 __file__, __line__)
1078 endif
1079
1080 ! dynamically create the Fortran output format
1081 write(sbuf,fmt="(a,a,a)") '(',trim( ascharacter(inumcols)),'(a,1x))'
1082
1083 open ( lu_temp, file=trim(output_grid_directory_name)//sfilename, iostat=istat, status="REPLACE" )
1084 call assert( istat==0, "Could not open output file "//dquote(sfilename), &
1085 __file__,__line__)
1086
1087 write ( unit=lu_temp, fmt="('NCOLS ',i10)", iostat=istat ) inumcols
1088 call assert( istat==0, "Error writing grid file header", __file__, __line__)
1089
1090 write ( unit=lu_temp, fmt="('NROWS ',i10)", iostat=istat ) inumrows
1091 call assert( istat==0, "Error writing grid file header", __file__, __line__)
1092
1093 write ( unit=lu_temp, fmt="('XLLCORNER ',f14.3)", iostat=istat ) pgrd%rX0
1094 call assert( istat==0, "Error writing X limits", __file__, __line__)
1095
1096 write ( unit=lu_temp, fmt="('YLLCORNER ',f14.3)", iostat=istat ) pgrd%rY0
1097 call assert( istat==0, "Error writing Y limits", __file__, __line__)
1098
1099 write ( unit=lu_temp, fmt="('CELLSIZE ',f14.3)", iostat=istat ) pgrd%rGridCellSize
1100 call assert( istat==0, "Error writing cell size", __file__, __line__)
1101
1102 if ( pgrd%iDataType == datatype_int ) then
1103
1104 write ( unit=lu_temp, fmt="('NODATA_VALUE ',i14)", iostat=istat ) pgrd%iNoDataValue
1105 call assert( istat==0, "Error writing NODATA value", __file__, __line__)
1106 do irow=1,inumrows
1107 write( unit=lu_temp, fmt=trim(sbuf), iostat=istat ) &
1108 (trim( ascharacter(pgrd%iData(icol,irow) ) ),icol=1,inumcols)
1109 call assert( istat==0, "Error writing Arc ASCII INTEGER grid data", &
1110 __file__, __line__)
1111 end do
1112
1113 elseif ( pgrd%iDataType == datatype_real ) then
1114
1115 write ( unit=lu_temp, fmt="('NODATA_VALUE ',g0.4)", iostat=istat ) pgrd%rNoDataValue
1116 call assert( istat==0, "Error writing NODATA value", __file__, __line__)
1117 do irow=1,inumrows
1118 write( unit=lu_temp, fmt=trim(sbuf), iostat=istat ) &
1119 (trim(ascharacter( pgrd%rData(icol,irow), fmt_string="g0.5" )),icol=1,inumcols)
1120 call assert( istat==0, "Error writing Arc ASCII REAL grid data", &
1121 __file__, __line__)
1122 end do
1123
1124 elseif ( pgrd%iDataType == datatype_double ) then
1125
1126 write ( unit=lu_temp, fmt="('NODATA_VALUE ',g0.4)", iostat=istat ) pgrd%dpNoDataValue
1127 call assert( istat==0, "Error writing NODATA value", __file__, __line__)
1128 do irow=1,inumrows
1129 write( unit=lu_temp, fmt=trim(sbuf), iostat=istat ) &
1130 (trim(ascharacter( pgrd%dpData(icol,irow), fmt_string="g0.5" )),icol=1,inumcols)
1131 call assert( istat==0, "Error writing Arc ASCII REAL grid data", &
1132 __file__, __line__)
1133 end do
1134
1135 endif
1136
1137 close (unit=lu_temp)
1138
1139end subroutine grid_writearcgrid
1140
1141!------------------------------------------------------------------------------
1142
1143subroutine grid_writesurfergrid(sFilename, pGrd)
1144
1145 ! [ ARGUMENTS ]
1146 character (len=*),intent(in) :: sfilename
1147 type (general_grid_t), pointer :: pgrd
1148
1149 ! [ LOCALS ]
1150 integer (c_int) :: icol,irow
1151 integer (c_int) :: inumcols, inumrows
1152 integer (c_int) :: istat
1153 character(len=256) :: sbuf
1154 real (c_float) :: fhalfcell
1155
1156 if ( pgrd%iDataType == datatype_int ) then
1157 inumcols = size(pgrd%iData,1)
1158 inumrows = size(pgrd%iData,2)
1159 elseif ( pgrd%iDataType == datatype_real ) then
1160 inumcols = size(pgrd%rData,1)
1161 inumrows = size(pgrd%rData,2)
1162 elseif ( pgrd%iDataType == datatype_double ) then
1163 inumcols = size(pgrd%dpData,1)
1164 inumrows = size(pgrd%dpData,2)
1165 else
1166 call assert(false, "Internal programming error - Unsupported grid type", &
1167 __file__, __line__)
1168 endif
1169
1170 fhalfcell = real(pgrd%rGridCellSize * 0.5_c_double, c_float)
1171
1172 ! dynamically create the Fortran output format
1173 write(sbuf,fmt="(a,a,a)") '(',trim( ascharacter(inumcols)),'(a,1x))'
1174
1175 open ( lu_temp, file=trim(output_grid_directory_name)//sfilename, iostat=istat, status="REPLACE" )
1176 call assert( istat==0, "Could not open output file "//dquote(sfilename), &
1177 __file__,__line__)
1178
1179 write ( unit=lu_temp, fmt="('DSAA')", iostat=istat )
1180 call assert( istat==0, "Error writing SURFER header", &
1181 __file__, __line__)
1182
1183 write ( unit=lu_temp, fmt="(2i8)", iostat=istat ) inumcols, inumrows
1184 call assert( istat==0, "Error writing SURFER dimensions", &
1185 __file__, __line__)
1186
1187 write ( unit=lu_temp, fmt="(2f14.3)", iostat=istat ) &
1188 pgrd%rX0 + fhalfcell , pgrd%rX1 - fhalfcell
1189 call assert( istat==0, "Error writing SURFER X limits", &
1190 __file__, __line__)
1191
1192 write ( unit=lu_temp, fmt="(2f14.3)", iostat=istat ) &
1193 pgrd%rY0 + fhalfcell, pgrd%rY1 - fhalfcell
1194 call assert( istat==0, "Error writing SURFER Y limits", &
1195 __file__, __line__)
1196
1197 if ( pgrd%iDataType == datatype_int ) then
1198
1199 write ( unit=lu_temp, fmt="(2i14)", iostat=istat ) minval(pgrd%iData),maxval(pgrd%iData)
1200 call assert( istat==0, "Error writing SURFER Z limits", &
1201 __file__, __line__)
1202
1203
1204 do irow=inumrows,1,-1
1205 write( unit=lu_temp, fmt=trim(sbuf), iostat=istat ) &
1206 (trim(ascharacter(pgrd%iData(icol,irow) ) ),icol=1,inumcols)
1207 call assert( istat==0, "Error writing SURFER grid data" , &
1208 __file__, __line__)
1209 end do
1210
1211 elseif ( pgrd%iDataType == datatype_real ) then
1212
1213 write ( unit=lu_temp, fmt="(2f14.3)", iostat=istat ) minval(pgrd%rData),maxval(pgrd%rData)
1214 call assert( istat==0, "Error writing SURFER Z limits", &
1215 __file__, __line__)
1216
1217 do irow=inumrows,1,-1
1218 write( unit=lu_temp, fmt=trim(sbuf), iostat=istat ) &
1219 (trim( ascharacter(pgrd%rData(icol,irow) ) ),icol=1,inumcols)
1220 call assert( istat==0, "Error writing SURFER grid data" , &
1221 __file__, __line__)
1222 end do
1223
1224 elseif ( pgrd%iDataType == datatype_double ) then
1225
1226 write ( unit=lu_temp, fmt="(2f14.3)", iostat=istat ) minval(pgrd%dpData),maxval(pgrd%dpData)
1227 call assert( istat==0, "Error writing SURFER Z limits", &
1228 __file__, __line__)
1229
1230 do irow=inumrows,1,-1
1231 write( unit=lu_temp, fmt=trim(sbuf), iostat=istat ) &
1232 (trim( ascharacter(pgrd%dpData(icol,irow) ) ),icol=1,inumcols)
1233 call assert( istat==0, "Error writing SURFER grid data" , &
1234 __file__, __line__)
1235 end do
1236
1237 endif
1238
1239 close (unit=lu_temp)
1240
1241end subroutine grid_writesurfergrid
1242!!***
1243!--------------------------------------------------------------------------
1244!!****f* grid/grid_Conform
1245! NAME
1246! grid_WriteArcGrid - Reads an Surfer ASCII grid of a specified type.
1247!
1248! SYNOPSIS
1249! Write an ARC ASCII grid of specified data type
1250!
1251! INPUTS
1252!
1253! pGrd1 - Pointer to first grid
1254! pGrd2 - Pointer to second grid
1255! rTolerance - OPTIONAL - Allowable tolerance between the two grids
1256!
1257! OUTPUTS
1258! lConform - Logical, TRUE if the grids are compatible within the
1259! specified tolerance, FALSE if otherwise
1260!
1261! SOURCE
1262
1263function grid_conform ( pGrd1, pGrd2, rTolerance ) result ( lConform )
1264 !! Returns .true. if the T_GRID objects conform (in terms of cell sizes and extents)
1265 !! The optional argument rTolerance is the precision for checking the (floating-point)
1266 !! extent coordinates (this defaults to rDEFAULT_TOLERANCE, below). The tolerance
1267 !! is set to abs(rTolerance * ( rX1-rX1 ) )
1268 ! ARGUMENTS
1269 type (general_grid_t), pointer :: pgrd1,pgrd2
1270 real (c_float), intent(in), optional :: rtolerance
1271 ! RETURN VALUE
1272 logical (c_bool) :: lconform
1273 ! LOCALS
1274 real (c_float) :: rtol
1275 real (c_float), parameter :: rdefault_tolerance = 0.5_c_float
1276
1277 if ( present ( rtolerance ) ) then
1278 rtol = real(rtolerance * ( pgrd1%rX1 - pgrd1%rX0 ), c_float)
1279 else
1280 rtol = real(rdefault_tolerance * ( pgrd1%rX1 - pgrd1%rX0 ), c_float)
1281 end if
1282
1283 lconform = true
1284
1285 if ( pgrd1%iNX /= pgrd2%iNX ) then
1286 lconform = false
1287 call logs%write("Unequal number of columns between grids", iloglevel=log_all)
1288 endif
1289
1290 if ( pgrd1%iNY /= pgrd2%iNY ) then
1291 lconform = false
1292 call logs%write("Unequal number of rows between grids", iloglevel=log_all)
1293 endif
1294
1295 if( abs( pgrd1%rX0 - pgrd2%rX0 ) > rtol ) then
1296 call logs%write("Lower left-hand side X coordinates do not match:", iloglevel=log_all)
1297 call logs%write("Grid 1 value: "//ascharacter(pgrd1%rX0) &
1298 //"; grid 2 value: "//ascharacter(pgrd2%rX0), iloglevel=log_all)
1299 lconform = false
1300 endif
1301
1302 if( abs( pgrd1%rY0 - pgrd2%rY0 ) > rtol ) then
1303 call logs%write("Lower left-hand side Y coordinates do not match:", iloglevel=log_all)
1304 call logs%write("Grid 1 value: "//ascharacter(pgrd1%rY0) &
1305 //"; grid 2 value: "//ascharacter(pgrd2%rY0), iloglevel=log_all)
1306 lconform = false
1307 endif
1308
1309 if( abs( pgrd1%rX1 - pgrd2%rX1 ) > rtol ) then
1310 call logs%write("Upper right-hand side X coordinates do not match:", iloglevel=log_all)
1311 call logs%write("Grid 1 value: "//ascharacter(pgrd1%rX1) &
1312 //"; grid 2 value: "//ascharacter(pgrd2%rX1), iloglevel=log_all)
1313 lconform = false
1314 endif
1315
1316 if( abs( pgrd1%rY1 - pgrd2%rY1 ) > rtol ) then
1317 call logs%write("Upper right-hand side Y coordinates do not match:", iloglevel=log_all)
1318 call logs%write("Grid 1 value: "//ascharacter(pgrd1%rY1) &
1319 //"; grid 2 value: "//ascharacter(pgrd2%rY1), iloglevel=log_all)
1320 lconform = false
1321 endif
1322
1323end function grid_conform
1324
1325!------------------------------------------------------------------------------
1326
1327function grid_completelycover( pBaseGrd, pOtherGrd, rTolerance ) result ( lCompletelyCover )
1328 !! Returns .true. if the T_GRID objects conform (in terms of cell sizes and extents)
1329 !! The optional argument rTolerance is the precision for checking the (floating-point)
1330 !! extent coordinates (this defaults to rDEFAULT_TOLERANCE, below). The tolerance
1331 !! is set to abs(rTolerance * ( rX1-rX1 ) )
1332 ! ARGUMENTS
1333 type (general_grid_t), pointer :: pbasegrd,pothergrd
1334 real (c_float), intent(in), optional :: rtolerance
1335 ! RETURN VALUE
1336 logical (c_bool) :: lcompletelycover
1337 ! LOCALS
1338 real (c_float) :: rtol
1339 real (c_float) :: rdefault_tolerance
1340
1341 rdefault_tolerance = real(-pbasegrd%rGridCellSize / 100.0_c_double, c_float)
1342 !rDEFAULT_TOLERANCE = rZERO
1343
1344 if ( present ( rtolerance ) ) then
1345 rtol = rtolerance
1346 else
1347 rtol = rdefault_tolerance
1348 end if
1349
1350 if ( (pbasegrd%rX0 - pothergrd%rX0 >= rtol) .and. &
1351 (pbasegrd%rY0 - pothergrd%rY0 >= rtol) .and. &
1352 (pothergrd%rX1 - pbasegrd%rX1 >= rtol) .and. &
1353 (pothergrd%rY1 - pbasegrd%rY1 >= rtol) ) then
1354 lcompletelycover = true
1355 else
1356 lcompletelycover = false
1357
1358 call logs%write("Extents of the data grid file "//dquote(pothergrd%sFilename) &
1359 //" do not cover the base grid extents.")
1360
1361 call logs%write("BASE GRID EXTENTS DATA GRID EXTENTS DIFFERENCE", &
1362 itab=18, ilinesbefore=1)
1363 call logs%write("X (lower-left): "//trim(ascharacter(pbasegrd%rX0)) &
1364 //" "//trim(ascharacter(pothergrd%rX0)) &
1365 //" "//trim(ascharacter(pbasegrd%rX0 - pothergrd%rX0)) )
1366 call logs%write("Y (lower-left): "//trim(ascharacter(pbasegrd%rY0)) &
1367 //" "//trim(ascharacter(pothergrd%rY0)) &
1368 //" "//trim(ascharacter(pbasegrd%rY0 - pothergrd%rY0)) )
1369 call logs%write("X (upper-right): "//trim(ascharacter(pbasegrd%rX1)) &
1370 //" "//trim(ascharacter(pothergrd%rX1)) &
1371 //" "//trim(ascharacter(pothergrd%rX1 - pbasegrd%rX1)) )
1372 call logs%write("Y (upper-right): "//trim(ascharacter(pbasegrd%rY1)) &
1373 //" "//trim(ascharacter(pothergrd%rY1)) &
1374 //" "//trim(ascharacter(pothergrd%rY1 - pbasegrd%rY1)), ilinesafter=1 )
1375
1376 end if
1377
1378end function grid_completelycover
1379
1380!!***
1381!--------------------------------------------------------------------------
1382!!****s* grid/grid_LookupColumn
1383! NAME
1384! grid_LookupColumn - Finds the column position of the value rXval
1385! in the grid
1386!
1387! SYNOPSIS
1388! Finds the column position of the value rXval in the grid and returns:
1389! iBefore = column before the value rYval
1390! iAfter = column after the value rYval
1391! rFrac = fraction of the distance between iBefore and iAfter
1392! If iBefore or iAfter are outside the grid, they're set to -1
1393!
1394! INPUTS
1395! pGrd - Pointer to a data grid
1396! rXval - Real value to be interpolated from the gridded data.
1397! rTolerance - OPTIONAL - Allowable tolerance between the two grids
1398!
1399! OUTPUTS
1400! iBefore - column before the value rXval
1401! iAfter - column after the value rXval
1402! rFrac - fraction of the distance between iBefore and iAfter
1403!
1404! If iBefore or iAfter are outside the grid, they're set to -1
1405!
1406! NOTES
1407! The grid lookup functions were designed specifically to work with the
1408! Thornthwaite-Mather soil moisture retention table. The functions will
1409! work with any data that are similarly arranged.
1410!
1411! If X=0.65, for example, and the following table were queried, the subroutine
1412! will return: iBefore=1, iAfter=2, rFrac=0.3
1413!
1414! X=Max SM Capacity
1415!
1416! 0.50 1.00 1.50
1417! _________________l
1418! 0.0 | 0.50 1.00 1.50
1419! Y=APWL 0.1 | 0.45 0.90 1.40
1420! 0.2 | 0.40 0.80 1.30
1421!
1422! SOURCE
1423
1424subroutine grid_lookupcolumn(pGrd,rXval,iBefore,iAfter,rFrac)
1425
1426 ! [ ARGUMENTS ]
1427 type ( GENERAL_GRID_T ),pointer :: pGrd
1428 real (c_float),intent(in) :: rXval
1429 integer (c_int),intent(out) :: iBefore,iAfter
1430 real (c_float),intent(out) :: rFrac
1431 ! [ LOCALS ]
1432 real (c_float) :: rColPosition
1433
1434 rcolposition = real((pgrd%iNX - 1) * (rxval - pgrd%rX0) / (pgrd%rX1 - pgrd%rX0), c_float)
1435 rfrac = rzero
1436 ibefore = floor(rcolposition) + 1
1437 iafter = ibefore + 1
1438 if ( ibefore > pgrd%iNX .or. ibefore < 1) then
1439 ibefore = -1
1440 iafter = -1
1441 else if ( iafter > pgrd%iNX ) then
1442 iafter = -1
1443 else
1444 rfrac = mod(rcolposition, 1.0_c_float )
1445 end if
1446
1447end subroutine grid_lookupcolumn
1448!!***
1449!--------------------------------------------------------------------------
1450!!****s* grid/grid_LookupRow
1451! NAME
1452! grid_LookupRow - Finds the row position of the value rYval
1453! in the grid
1454!
1455! SYNOPSIS
1456! Finds the column position of the value rYval in the grid and returns:
1457! iBefore = column before the value rYval
1458! iAfter = column after the value rYval
1459! rFrac = fraction of the distance between iBefore and iAfter
1460! If iBefore or iAfter are outside the grid, they're set to -1
1461!
1462! INPUTS
1463! pGrd - Pointer to a data grid
1464! rYval - Real value to be interpolated from the gridded data.
1465! rTolerance - OPTIONAL - Allowable tolerance between the two grids
1466!
1467! OUTPUTS
1468! iBefore - column before the value rXval
1469! iAfter - column after the value rXval
1470! rFrac - fraction of the distance between iBefore and iAfter
1471!
1472! If iBefore or iAfter are outside the grid, they're set to -1
1473!
1474! NOTES
1475! The grid lookup functions were designed specifically to work with the
1476! Thornthwaite-Mather soil moisture retention table. The functions will
1477! work with any data that are similarly arranged.
1478!
1479! If Y=0.15, for example, and the following table were queried, the subroutine
1480! will return: iBefore=2, iAfter=3, rFrac=0.5
1481!
1482! X=Max SM Capacity
1483!
1484! 0.50 1.00 1.50
1485! _________________l
1486! 0.0 | 0.50 1.00 1.50
1487! Y=APWL 0.1 | 0.45 0.90 1.40
1488! 0.2 | 0.40 0.80 1.30
1489!
1490! SOURCE
1491
1492subroutine grid_lookuprow(pGrd,rYval,iBefore,iAfter,rFrac)
1493
1494 ! [ ARGUMENTS ]
1495 type ( GENERAL_GRID_T ),pointer :: pGrd
1496 real (c_float),intent(in) :: rYval
1497 integer (c_int),intent(out) :: iBefore,iAfter
1498 real (c_float),intent(out) :: rFrac
1499 ! [ LOCALS ]
1500 real (c_float) :: rColPosition
1501
1502 rcolposition = real((pgrd%iNY - 1) * (pgrd%rY1 - ryval) / (pgrd%rY1 - pgrd%rY0), c_float)
1503 rfrac = rzero
1504 ibefore = floor(rcolposition) + 1
1505 iafter = ibefore + 1
1506 if ( ibefore > pgrd%iNY ) then
1507 ibefore = -1
1508 iafter = -1
1509 else if ( ibefore < 1) then
1510 ibefore = 1
1511 iafter = 1
1512 else if ( iafter > pgrd%iNY ) then
1513 iafter = ibefore
1514 else
1515 rfrac = mod(rcolposition, 1.0_c_float )
1516 end if
1517
1518end subroutine grid_lookuprow
1519
1520!------------------------------------------------------------------------------
1521
1522!> Call PROJ4 to transform coordinates.
1523!! @details This subroutine calls a Fortran wrapper to the C library
1524!! PROJ4. A set of input coordinates is transformed to the base SWB
1525!! coordinate system.
1526!!
1527!! The idea is to first create a 2D array of the X and Y coordinates in
1528!! the projection system of the input data set. These coordinate values are
1529!! then fed to PROJ4, which modifies the values and returns the coordinate
1530!! values in the projection system of the base grid. Thereafter, mapping
1531!! source to target is a matter of finding the value of the cell closest to
1532!! a SWB grid coordinate pair.
1533!! @param[inout] pGrd
1534!! @param[in] sFromPROJ4
1535!! @param[in] sToPROJ4
1536
1537subroutine grid_transform(pGrd, sFromPROJ4, sToPROJ4, rX, rY )
1538
1539 type ( general_grid_t ),pointer :: pgrd
1540 character (len=*) :: sfromproj4, stoproj4
1541 character (kind=c_char, len=len_trim(sFromPROJ4)) :: csfromproj4
1542 character (kind=c_char, len=len_trim(sToPROJ4)) :: cstoproj4
1543 real (c_double), optional :: rx(:)
1544 real (c_double), optional :: ry(:)
1545
1546 ! [ LOCALS ]
1547 integer (c_int) :: iretval
1548
1549 if (present(rx) .and. present(ry)) then
1550
1551 !> supply X and Y for the untransformed data
1552 call grid_populatexy(pgrd, rx, ry)
1553
1554 else
1555
1556 !> calculate X and Y for the untransformed data
1557 call grid_populatexy(pgrd)
1558
1559 endif
1560
1561 csfromproj4 = trim(sfromproj4)
1562 cstoproj4 = trim(stoproj4)
1563
1564 !> PROJ4 expects unprojected coordinates (i.e. lat lon) to be provided
1565 !> in RADIANS. Therefore, we convert to radians prior to the call...
1566
1567 if ( ( csfromproj4 .containssimilar. "latlon" ) &
1568 .or. ( csfromproj4 .containssimilar. "latlong" ) &
1569 .or. ( csfromproj4 .containssimilar. "lonlat" ) &
1570 .or. ( csfromproj4 .containssimilar. "longlat" ) ) then
1571
1572 pgrd%rX = pgrd%rX * degrees_to_radians
1573 pgrd%rY = pgrd%rY * degrees_to_radians
1574
1575 endif
1576
1577 iretval = pj_init_and_transform(csfromproj4//c_null_char, &
1578 cstoproj4//c_null_char, &
1579 __file__//c_null_char, &
1580 __line__, &
1581 int(pgrd%iNumGridCells, c_long), pgrd%rX, pgrd%rY)
1582
1583 call grid_checkforproj4error(iretval, sfromproj4, stoproj4)
1584
1585
1586 !> If the coordinates have been converted TO latlon, convert back to degrees
1587
1588 if( index(string=cstoproj4, substring="latlon") > 0 &
1589 .or. index(string=cstoproj4, substring="lonlat") > 0 ) then
1590
1591 pgrd%rX = pgrd%rX * radians_to_degrees
1592 pgrd%rY = pgrd%rY * radians_to_degrees
1593
1594 endif
1595
1596 ! now update the grid boundaries based on the transformed coordinate values
1597 pgrd%rGridCellSize = ( maxval(pgrd%rX) - minval(pgrd%rX) ) &
1598 / real(pgrd%iNX - 1, c_double)
1599
1600 pgrd%rX0 = minval(pgrd%rX) - pgrd%rGridCellSize / 2.0_c_double
1601 pgrd%rX1 = maxval(pgrd%rX) + pgrd%rGridCellSize / 2.0_c_double
1602 pgrd%rY0 = minval(pgrd%rY) - pgrd%rGridCellSize / 2.0_c_double
1603 pgrd%rY1 = maxval(pgrd%rY) + pgrd%rGridCellSize / 2.0_c_double
1604
1605 ! finally, change the projection string to reflect the new coordinate system
1606 pgrd%sPROJ4_string = trim(stoproj4)
1607
1608 ! at this point, what we have is the same old values that were read into the
1609 ! input grid, but a new set of coordinate and boundary definitions that reference
1610 ! the SWB grid projection.
1611
1612end subroutine grid_transform
1613
1614!--------------------------------------------------------------------------
1615
1616subroutine grid_checkintegergridvalues(pGrd, sFilename)
1617
1618 type ( GENERAL_GRID_T ),pointer :: pGrd
1619 character (len=*), intent(in) :: sFilename
1620
1621 ! [ LOCALS ]
1622 integer (c_int) :: iRunningSum
1623 integer (c_int) :: iRecord
1624
1625 irunningsum = 0
1626 irecord = 0
1627
1628! call LOGS%set_echo(FALSE)
1629! call LOGS%write("### Summary of integer grid data values for file "//dquote(sFilename)//" ###", &
1630! iLogLevel=LOG_DEBUG, iLinesBefore=1, iLinesAfter=1 )
1631
1632! call LOGS%write("number | count | value ", iLogLevel=LOG_DEBUG )
1633! call LOGS%write("---------- | -------------- | ----------", iLogLevel=LOG_DEBUG )
1634! do iIndex=0,maxval(pGrd%iData)
1635! iCount=COUNT( pGrd%iData==iIndex )
1636! if ( iCount > 0 ) then
1637! iRecord = iRecord + 1
1638! iRunningSum = iRunningSum + iCount
1639! write (sBuf0, fmt="(i10)") iRecord
1640! write (sBuf1, fmt="(i14)") iCount
1641! write (sBuf2, fmt="(i10)") iIndex
1642! write (sBuf3, fmt="(a10,' | ', a14,' | ',a10)") adjustl(sBuf0), adjustl(sBuf1), adjustl(sBuf2)
1643! call LOGS%write( sBuf3, iLogLevel=LOG_DEBUG )
1644! end if
1645! end do
1646
1647! call LOGS%write(" Total number of grid cells with value NODATA: " &
1648! //asCharacter( COUNT(pGrd%iData == pGrd%iNoDataValue ) ), iLinesBefore=1, iLogLevel=LOG_ALL )
1649
1650! call LOGS%write(" Total number of grid cells: "//asCharacter( size(pGrd%iData) ), iLogLevel=LOG_ALL )
1651
1652! call LOGS%write(" Total number of grid cells with value >= 0: "//asCharacter(iRunningSum), iLogLevel=LOG_ALL )
1653
1654
1655! if (size(pGrd%iData) /= iRunningSum) then
1656! call LOGS%write(repeat("*",80), iLogLevel=LOG_ALL)
1657! call LOGS%write("Possible illegal or missing values in integer grid file: "//trim(sFileName), iLogLevel=LOG_ALL)
1658! call LOGS%write(repeat("*",80), iLogLevel=LOG_ALL)
1659! endif
1660
1661end subroutine grid_checkintegergridvalues
1662
1663!--------------------------------------------------------------------------
1664
1665subroutine grid_checkforproj4error(iRetVal, sFromPROJ4, sToPROJ4)
1666
1667 ! [ ARGUMENTS ]
1668 integer (c_int) :: iretval
1669 character (len=*) :: sfromproj4
1670 character (len=*) :: stoproj4
1671
1672 ! [ LOCALS ]
1673 integer (c_int) :: i
1674 logical (c_bool) :: lfound
1675 character (len=256) :: serrormessage
1676
1677 serrormessage = ""
1678
1679 if (iretval /= 0) then
1680
1681 write(serrormessage,fmt="(a,a,a,a,a)") &
1682 "There was an error transforming a grid from this projection:~", &
1683 dquote(sfromproj4), " to:~", &
1684 dquote(stoproj4)
1685
1686 do i=1,49
1687 if(iretval == terrormessage(i)%iErrorNumber) then
1688 serrormessage = serrormessage &
1689 //"~ PROJ4 Error number "// ascharacter(iretval)//" reported.~" &
1690 //" ==> error description: "//trim(terrormessage(i)%cErrorMessageText)
1691 lfound = true
1692 exit
1693 endif
1694 enddo
1695
1696 call assert(false, trim(serrormessage), __file__, __line__)
1697
1698 endif
1699
1700end subroutine grid_checkforproj4error
1701
1702!!***
1703!--------------------------------------------------------------------------
1704!!****f* grid/grid_Interpolate
1705! NAME
1706! grid_Interpolate - Returns an interpolated table value given a row
1707! and column position
1708!
1709
1710! SYNOPSIS
1711! Returns an interpolated table value given a row and column position
1712!
1713! INPUTS
1714! pGrd - Pointer to a data grid
1715! rXval - X value for which we want to interpolate table data.
1716! rYval - Y value for which we want to interpolate table data.
1717!
1718! OUTPUTS
1719! rValue - Value interpolated from input grid.
1720!
1721! NOTES
1722! The grid lookup functions were designed specifically to work with the
1723! Thornthwaite-Mather soil moisture retention table. The functions will
1724! work with any data that are similarly arranged.
1725!
1726! SOURCE
1727
1728function grid_interpolate(pGrd,rXval,rYval) result ( rValue )
1729 !! Interpolates values from the grid 'grd' for the row position 'yval' and
1730 !! the column position 'xval'. Assumes that the row and column spacing
1731 !! are constant. Applicable only to DATATYPE_REAL grids.
1732 ! [ ARGUMENTS ]
1733 type ( general_grid_t ),pointer :: pgrd
1734 real (c_float),intent(in) :: rxval,ryval
1735 ! [ RETURN VALUE ]
1736 real (c_float) :: rvalue
1737 ! [ LOCALS ]
1738 integer (c_int) :: ib,jb,ia,ja
1739 real (c_float) :: ylocal,u,v
1740
1741 call grid_lookupcolumn(pgrd,rxval,ib,ia,u)
1742
1743 call assert (ib>0 .and. ia>0 .and. ib <= ubound(pgrd%rData,1) &
1744 .and. ia <= ubound(pgrd%rData,1), &
1745 "Internal programming error: illegal bounds caught~requested column value " &
1746 //trim( ascharacter(rxval, fmt_string="F0.3")) &
1747 //" out of range", __file__, __line__)
1748
1749 ! In some cases, when things really dry out, the y value
1750 ! goes out of range - enforce bounds.
1751 if ( ryval < pgrd%rY0 ) then
1752 ylocal = real(pgrd%rY0, c_float)
1753 else if ( ryval > pgrd%rY1 ) then
1754 ylocal = real(pgrd%rY1, c_float)
1755 else
1756 ylocal = ryval
1757 end if
1758
1759 call grid_lookuprow(pgrd=pgrd, &
1760 ryval=ylocal, &
1761 ibefore=jb, &
1762 iafter=ja, &
1763 rfrac=v)
1764
1765 call assert (jb>0 .and. ja>0 .and. jb <= ubound(pgrd%rData,2) &
1766 .and. ja <= ubound(pgrd%rData,2), &
1767 "Internal programming error: illegal bounds caught~requested row value " &
1768 //trim( ascharacter(rxval, fmt_string="F0.3")) &
1769 //" out of range", __file__, __line__)
1770
1771 rvalue = ( 1.0_c_float -u) * ( 1.0_c_float -v) * pgrd%rData(ib,jb) + &
1772 u * ( 1.0_c_float -v) * pgrd%rData(ib,ja) + &
1773 ( 1.0_c_float -u) * v * pgrd%rData(ia,jb) + &
1774 u * v * pgrd%rData(ia,ja)
1775
1776end function grid_interpolate
1777!!***
1778!--------------------------------------------------------------------------
1779!!****f* grid/grid_SearchColumn
1780! NAME
1781! grid_SearchColumn - Returns the value of table Y value given the X value
1782! and the table Z value
1783!
1784! SYNOPSIS
1785! Searches in the y-direction, given the x value and the table value
1786! z, and returns the value of y. Assumes that the row and
1787! column spacing are constant. The search begins at the _top_l of the
1788! table (y=pGrd%rY1). Applicable only to DATATYPE_REAL grids.
1789! Parameter 'rNoData' is the missing value code for the grid.
1790
1791! INPUTS
1792! pGrd - Pointer to a data grid
1793! rXval - X value corresponding to one or more rows of the table.
1794! rZval - Table value to scan for.
1795! rNoData - Real value assigned to "No Data" cells.
1796!
1797! OUTPUTS
1798! rYval - Returned value associated with the input x and z values.
1799!
1800! NOTES
1801! The grid lookup functions were designed specifically to work with the
1802! Thornthwaite-Mather soil moisture retention table. The functions will
1803! work with any data that are similarly arranged.
1804!
1805! If X=1.00 and Z=0.90, for example, and the following table were queried,
1806! the subroutine will return: rValue=0.1
1807!
1808! X=Max SM Capacity
1809!
1810! 0.50 1.00 1.50
1811! _________________l
1812! 0.0 | 0.50 1.00 1.50
1813! Y=APWL 0.1 | 0.45 0.90 1.40
1814! 0.2 | 0.40 0.80 1.30
1815!
1816! SOURCE
1817
1818function grid_searchcolumn(pGrd,rXval,rZval,rNoData) result ( rValue )
1819 !! Searches in the y-direction, given the x value 'xval', for the value
1820 !! 'zval', and returns the value of y in 'ry'. Assumes that the row and
1821 !! column spacing are constant. The search begins at the _top_l of the
1822 !! table (y=pGrd%rY1). Applicable only to DATATYPE_REAL grids.
1823 !! Parameter 'rmv' is the missing value code for the grid.
1824 ! [ ARGUMENTS ]
1825 type ( general_grid_t ),pointer :: pgrd
1826 real (c_float),intent(in) :: rxval,rzval,rnodata
1827 ! [ RETURN VALUE ]
1828 real (c_float) :: rvalue
1829 ! [ LOCALS ]
1830 integer (c_int) :: ib,ia,istat,irow
1831 real (c_float) :: u,v,frac,rprev
1832 real (c_float),dimension(:),allocatable :: rcol
1833
1834 ! allocate space for all ROWS of data that come from a single COLUMN
1835 allocate ( rcol(pgrd%iNY), stat=istat )
1836 call assert( LOGICAL(istat==0,c_bool), &
1837 "Couldn't allocate temporary array" )
1838
1839 call grid_lookupcolumn(pgrd=pgrd, &
1840 rxval=rxval, &
1841 ibefore=ib, &
1842 iafter=ia, &
1843 rfrac=u)
1844
1845 call assert (ib>0 .and. ia>0 .and. ib <= ubound(pgrd%rData,1) &
1846 .and. ia <= ubound(pgrd%rData,1), &
1847 "Internal programming error: requested X value " &
1848 //trim( ascharacter(rxval, fmt_string="F0.3")) &
1849 //" out of range", __file__, __line__)
1850
1851 call assert(ubound(rcol,1) == ubound(pgrd%rData,2), &
1852 "Internal programming error: upper bound of rCol /= upper " &
1853 //"bound of first array element of pGrd%rData~" &
1854 //"ubound(rCol)="//trim( ascharacter(ubound(rcol,1))) &
1855 //"~ubound(pGrd%rData)="//trim( ascharacter(ubound(pgrd%rData,1))), &
1856 __file__, __line__)
1857
1858 ! interpolate the column of values based on the columns of values
1859 ! that bracket the real value rXval
1860! rCol = u*pGrd%rData(:,ia) + ( 1.0_c_float -u)*pGrd%rData(:,ib)
1861 rcol = u * pgrd%rData(ia,:) + ( 1.0_c_float -u)*pgrd%rData(ib,:)
1862 ! Fix missing values
1863
1864 do irow=1,pgrd%iNY
1865 if ( pgrd%rData(ia, irow) <= rnodata .and. pgrd%rData(ib, irow) <= rnodata ) then
1866 rcol(irow) = rnodata
1867 else if ( pgrd%rData(ib,irow) <= rnodata .and. u>0.9_c_float ) then
1868 rcol(irow) = pgrd%rData(ia,irow)
1869 else if ( pgrd%rData(ia,irow) <= rnodata .and. u<0.1_c_float ) then
1870 rcol(irow) = pgrd%rData(ib,irow)
1871 else if ( pgrd%rData(ia,irow) <= rnodata .or. pgrd%rData(ib,irow) <= rnodata ) then
1872 rcol(irow) = rnodata
1873 end if
1874 end do
1875
1876 ! Search for the specified value, assuming that the data vary monotonically
1877 ! Skip MVs
1878 rprev = rnodata
1879 rvalue = rnodata
1880 do irow=1,pgrd%iNY
1881 if ( rcol(irow) <= rnodata ) then
1882 if ( rprev <= rnodata ) then
1883 ! keep skipping missing values...
1884 continue
1885 else
1886 ! end of the line -- we didn't find the value, so return the limit
1887 v = real(irow-1, c_float) / real(pgrd%iNY-1, c_float)
1888 rvalue = real(( 1.0_c_float -v)*pgrd%rY0 + v*pgrd%rY1, c_float)
1889 exit
1890 endif
1891 else
1892 if ( rprev <= rnodata ) then
1893 rprev = rcol(irow)
1894 else
1895 ! Detect sign change: a negative product means rCol(iRow) and rprev
1896 ! straddle rZval (one above, one below), indicating the target value
1897 ! lies between these two consecutive table entries.
1898 if ( (rcol(irow) - rzval) * (rprev - rzval) < 0.0_c_float ) then
1899 ! Found it!
1900 frac = (rzval-rcol(irow-1)) / (rcol(irow)-rcol(irow-1))
1901 ! Note: it's 'iRow-2' in the next expr to account for one-based indexing
1902 v = ( real(irow-2, c_float)+frac ) / real(pgrd%iNY-1, c_float)
1903 rvalue = real(v*pgrd%rY0 + ( 1.0_c_float -v)*pgrd%rY1, c_float)
1904 exit
1905 else if ( rzval < rprev .and. rcol(irow) > rprev ) then
1906 ! does not exist, choose the limit
1907 v = real(irow-2, c_float) / real(pgrd%iNY-1, c_float)
1908 rvalue = real(v*pgrd%rY0 + ( 1.0_c_float -v)*pgrd%rY1, c_float)
1909 exit
1910 else if ( rzval > rprev .and. rcol(irow) < rprev ) then
1911 v = real(irow-2, c_float) / real(pgrd%iNY-1, c_float)
1912 rvalue = real(v*pgrd%rY0 + ( 1.0_c_float -v)*pgrd%rY1, c_float)
1913 exit
1914 else if ( irow == pgrd%iNY ) then
1915 v = 1.0_c_float
1916 rvalue = real(v*pgrd%rY0 + ( 1.0_c_float -v)*pgrd%rY1, c_float)
1917 exit
1918 end if
1919 end if
1920 end if
1921 end do
1922
1923 deallocate ( rcol )
1924
1925end function grid_searchcolumn
1926!!***
1927!--------------------------------------------------------------------------
1928!!****f* grid/grid_LookupReal
1929! NAME
1930! grid_LookupReal - Returns the grid cell value nearest to that for a
1931! given a row and column position
1932!
1933! SYNOPSIS
1934! Returns the grid cell value nearest to that for a given a row and
1935! column position. No interpolation is performed.
1936!
1937! INPUTS
1938! pGrd - Pointer to a data grid
1939! rXval - X value for which we want to interpolate table data.
1940! rYval - Y value for which we want to interpolate table data.
1941!
1942! OUTPUTS
1943! rValue - Value of the grid cell nearest the given row, column combination.
1944!
1945! NOTES
1946! The grid lookup functions were designed specifically to work with the
1947! Thornthwaite-Mather soil moisture retention table. The functions will
1948! work with any data that are similarly arranged.
1949!
1950! SOURCE
1951
1952function grid_lookupreal(pGrd,rXval,rYval) result(rValue)
1953 !! Returns the grid value for the cell containing (rXval,rYval)
1954 type ( general_grid_t ),pointer :: pgrd
1955 real (c_float),intent(in) :: rxval,ryval
1956 real (c_float) :: rvalue
1957 integer (c_int) :: icol,irow
1958
1959 irow = int( pgrd%iNY * (pgrd%rY0 - ryval) / (pgrd%rY0 - pgrd%rY1) ) + 1
1960 if ( irow > pgrd%iNY ) irow = pgrd%iNY
1961 icol = int( pgrd%iNX * (rxval - pgrd%rX0) / (pgrd%rX1 - pgrd%rX0) ) + 1
1962 if ( icol > pgrd%iNX ) icol = pgrd%iNX
1963 rvalue = pgrd%rData(icol,irow)
1964
1965end function grid_lookupreal
1966!!***
1967
1968!----------------------------------------------------------------------
1969
1970function grid_coordinatesfallinsidegrid(pGrd, rX, rY) result(coordinates_fall_inside_grid)
1971
1972 type ( general_grid_t ),pointer :: pgrd
1973 real (c_double) :: rx
1974 real (c_double) :: ry
1975 logical (c_bool) :: coordinates_fall_inside_grid
1976
1977 coordinates_fall_inside_grid = ( ( rx >= pgrd%rX0 ) .and. ( rx <= pgrd%rX1 ) &
1978 .and. ( ry >= pgrd%rY0 ) .and. ( ry <= pgrd%rY1 ) )
1979
1980end function grid_coordinatesfallinsidegrid
1981
1982!----------------------------------------------------------------------
1983
1984function grid_rowcolfallsinsidegrid(pGrd, iRow, iCol) result(row_col_falls_inside_grid)
1985
1986 type ( general_grid_t ),pointer :: pgrd
1987 integer (c_int) :: irow
1988 integer (c_int) :: icol
1989 logical (c_bool) :: row_col_falls_inside_grid
1990
1991 row_col_falls_inside_grid = ( ( irow >= 1 ) .and. ( irow <= pgrd%iNY ) &
1992 .and. ( icol >= 1 ) .and. ( icol <= pgrd%iNX ) )
1993
1994end function grid_rowcolfallsinsidegrid
1995
1996!----------------------------------------------------------------------
1997
1998function grid_getgridcolnum(pGrd,rX) result(iColumnNumber)
1999
2000 type ( general_grid_t ),pointer :: pgrd
2001 real (c_double) :: rx
2002 integer (c_int) :: icolumnnumber
2003
2004 !! this only works if the data are in the originally supplied projection. If the coordinates have been
2005 !! transformed, there is no guarantee that the assumption used in this calculation will hold.
2006 icolumnnumber = nint(real(pgrd%iNX, c_double) &
2007 * ( rx - pgrd%rX0 ) / (pgrd%rX1 - pgrd%rX0) + 0.5_c_double, c_int)
2008
2009! * ( rX - pGrd%rX0 ) / (pGrd%rX1 - pGrd%rX0) + 0.5_c_double, c_int)
2010
2011 ! print *, "iColumnNumber = ", iColumnNumber
2012 ! print *, "calc: ", real(pGrd%iNX, c_double) &
2013 ! * ( rX - pGrd%rX0 ) / (pGrd%rX1 - pGrd%rX0) + 0.5_c_double
2014 ! print *, "numerator: ", real(pGrd%iNX, c_double) * ( rX - pGrd%rX0 )
2015 ! print *, "numerator (LHS): ", real(pGrd%iNX, c_double)
2016 ! print *, "numerator (RHS): ", rX, pGrd%rX0, ( rX - pGrd%rX0 )
2017 ! print *, "denominator: ", (pGrd%rX1 - pGrd%rX0)
2018 !
2019! if ( iColumnNumber < 1 .or. iColumnNumber > pGrd%iNX ) then
2020! call grid_DumpGridExtent(pGrd)
2021! write(*, fmt="(a)") "was attempting to find column associated with X: "//trim(asCharacter(rX))
2022! call assert(FALSE, "INTERNAL PROGRAMMING ERROR: Column number out of bounds (value: " &
2023! //trim( asCharacter(iColumnNumber))//")", &
2024! __FILE__, __LINE__)
2025! endif
2026
2027end function grid_getgridcolnum
2028
2029!----------------------------------------------------------------------
2030
2031function grid_getgridrownum(pGrd,rY) result(iRowNumber)
2032
2033 type ( general_grid_t ),pointer :: pgrd
2034 real (c_double) :: ry
2035 integer (c_int) :: irownumber
2036
2037 irownumber = pgrd%iNY - nint(real(pgrd%iNY, c_double) &
2038 * ( ry - pgrd%rY0 ) / (pgrd%rY1 - pgrd%rY0) - 0.5_c_double, c_int)
2039
2040! if ( iRowNumber < 1 .or. iRowNumber > pGrd%iNY ) then
2041! call grid_DumpGridExtent(pGrd)
2042! write(*, fmt="(a)") "was attempting to find row associated with Y: "//trim(asCharacter(rY))
2043! call assert(FALSE, "INTERNAL PROGRAMMING ERROR: Row number out of bounds (value: " &
2044! //trim( asCharacter(iRowNumber))//")", &
2045! __FILE__, __LINE__)
2046! endif
2047
2048end function grid_getgridrownum
2049
2050!----------------------------------------------------------------------
2051
2052function grid_getgridcolrownum(pGrd, rX, rY) result(iColRow)
2053
2054 type ( general_grid_t ),pointer :: pgrd
2055 real (c_double) :: rx, ry
2056 integer (c_int), dimension(2) :: icolrow
2057
2058 ! [ LOCALS ]
2059 real (c_float) :: rdist, rmindistance, rdist2
2060
2061 integer (c_int), save :: ilastcolnum, ilastrownum
2062 integer (c_int) :: istartingcolnum, istartingrownum
2063 integer (c_int) :: irowboundlower, irowboundupper
2064 integer (c_int) :: icolboundlower, icolboundupper
2065 integer (c_int) :: icol, irow
2066 integer (c_int) :: icandidaterow, icandidatecol
2067 logical (c_bool) :: lchanged
2068
2069 ! this will get us close to where we need to be; however, since the
2070 ! transformed grid may contain significant nonlinearity between
2071 ! adjacent grid coordinates, we need to do a more thorough search
2072 ! in the neighborhood of these points
2073
2074 istartingcolnum = grid_getgridcolnum(pgrd,rx)
2075 istartingrownum = grid_getgridrownum(pgrd,ry)
2076
2077 if (istartingcolnum > 0 .and. istartingcolnum <= pgrd%iNX &
2078 .and. istartingrownum > 0 .and. istartingrownum <= pgrd%iNY) then
2079 rdist = real(hypot(pgrd%rX(istartingcolnum,istartingrownum) - rx, &
2080 pgrd%rY(istartingcolnum,istartingrownum) - ry), c_float)
2081
2082 ! need to ensure that the leftover column and row numbers from
2083 ! perhaps an entirely different grid are not used as indices
2084 if (ilastcolnum > 0 .and. ilastcolnum <= pgrd%iNX &
2085 .and. ilastrownum > 0 .and. ilastrownum <= pgrd%iNY) then
2086 rdist2 = real(hypot(pgrd%rX(ilastcolnum,ilastrownum) - rx, &
2087 pgrd%rY(ilastcolnum,ilastrownum) - ry), c_float)
2088 else
2089 rdist2 = rbigval
2090 endif
2091
2092 if (rdist > rdist2) then
2093 icandidaterow = ilastrownum
2094 icandidatecol = ilastcolnum
2095 else
2096 icandidaterow = istartingrownum
2097 icandidatecol = istartingcolnum
2098 endif
2099
2100! print *, 'filename: ', trim(pGrd%sFilename)
2101! print *, 'start of iteration: ', iStartingColNum, iStartingRowNum, rDist, rDist2, iCandidateCol, iCandidateRow
2102
2103 rmindistance = rbigval
2104
2105 do
2106
2107 !> need to ensure that whatever bound is calculated
2108 !> is within the declared array bounds or we get a segfault
2109 !iRowBoundLower = min(max( 1, iCandidateRow - 1), pGrd%iNY)
2110 irowboundlower = clip(value=icandidaterow - 1, minval=1, maxval=pgrd%iNY)
2111 irowboundupper = clip(value=icandidaterow + 1, minval=1, maxval=pgrd%iNY)
2112
2113 !iColBoundLower = min(max( 1, iCandidateCol - 1), pGrd%iNX)
2114 icolboundlower = clip(value=icandidatecol - 1, minval=1, maxval=pgrd%iNX)
2115 icolboundupper = clip(value=icandidatecol + 1, minval=1, maxval=pgrd%iNX)
2116
2117 lchanged = false
2118
2119 do irow=irowboundlower,irowboundupper
2120 do icol=icolboundlower,icolboundupper
2121
2122 rdist = real(hypot(pgrd%rX(icol,irow) - rx, pgrd%rY(icol,irow) - ry), c_float)
2123
2124 if (rdist < rmindistance ) then
2125
2126 rmindistance = rdist
2127 icandidatecol = icol
2128 icandidaterow = irow
2129 lchanged = true
2130
2131 endif
2132
2133 enddo
2134 enddo
2135
2136! print *, 'iterating: rdist: ', rDist, ' cand X: ', &
2137! pGrd%rX(iCandidateCol,iCandidateRow), ' X: ', rX, &
2138! ' cand Y: ', pGrd%rY(iCandidateCol,iCandidateRow), ' Y: ', rY
2139
2140 if (.not. lchanged ) exit
2141
2142 enddo
2143
2144 ilastcolnum = icandidatecol
2145 ilastrownum = icandidaterow
2146
2147 icolrow(column) = icandidatecol
2148 icolrow(row) = icandidaterow
2149
2150 else
2151
2152 icolrow(column) = istartingcolnum
2153 icolrow(row) = istartingrownum
2154 endif
2155
2156end function grid_getgridcolrownum
2157
2158!----------------------------------------------------------------------
2159subroutine grid_set_nodata_value( pGrd, iValue, fValue )
2160 type ( general_grid_t ),pointer :: pgrd
2161 integer (c_int), intent(in), optional :: ivalue
2162 real (c_float), intent(in), optional :: fvalue
2163
2164 if ( present( ivalue ) ) pgrd%iNoDataValue = ivalue
2165 if ( present( fvalue ) ) pgrd%rNoDataValue = fvalue
2166end subroutine grid_set_nodata_value
2167!--------------------------------------------------------------------------------------------------
2168
2169function grid_getgridx(pGrd,iColumnNumber) result(rX)
2170
2171 type ( general_grid_t ),pointer :: pgrd
2172 real (c_double) :: rx
2173 integer (c_int) :: icolumnnumber
2174
2175 rx = pgrd%rX0 + pgrd%rGridCellSize * (real(icolumnnumber, c_double) - 0.5_c_double )
2176
2177end function grid_getgridx
2178
2179!----------------------------------------------------------------------
2180
2181function grid_getgridy(pGrd,iRowNumber) result(rY)
2182
2183 type ( general_grid_t ),pointer :: pgrd
2184 real (c_double) :: ry
2185 integer (c_int) :: irownumber
2186
2187 ry = pgrd%rY1 &
2188 - pgrd%rGridCellSize * (real(irownumber, c_double) - 0.5_c_double )
2189
2190end function grid_getgridy
2191
2192!----------------------------------------------------------------------
2193
2194subroutine grid_populatexy(pGrd, rX, rY)
2195
2196 ! [ ARGUMENTS ]
2197 type ( general_grid_t ),pointer :: pgrd
2198 real (c_double), optional :: rx(:)
2199 real (c_double), optional :: ry(:)
2200
2201 ! [ LOCALS ]
2202 integer (c_int) :: icol,irow
2203 integer (c_int) :: istat
2204 integer (c_int) :: nx, ny
2205
2206 if ( .not. allocated(pgrd%rX) ) then
2207
2208 ALLOCATE (pgrd%rX(pgrd%iNX, pgrd%iNY), stat=istat)
2209 call assert( istat == 0, &
2210 "Could not allocate memory for x-coordinates within grid data structure", &
2211 __file__, __line__)
2212 endif
2213
2214 if ( .not. allocated(pgrd%rY) ) then
2215 ALLOCATE (pgrd%rY(pgrd%iNX, pgrd%iNY), stat=istat)
2216 call assert( istat == 0, &
2217 "Could not allocate memory for y-coordinates within grid data structure", &
2218 __file__, __line__)
2219 endif
2220
2221 if (present(rx) .and. present(ry)) then
2222
2223 nx = size(rx)
2224 ny = size(ry)
2225
2226 call assert( nx == pgrd%iNX, &
2227 "Internal programming error: number of 'X' values does not match number of x values in grid", &
2228 __file__, __line__)
2229
2230 call assert( ny == pgrd%iNY, &
2231 "Internal programming error: number of 'X' values does not match number of x values in grid", &
2232 __file__, __line__)
2233
2234 do irow=1,pgrd%iNY
2235 do icol=1,pgrd%iNX
2236 pgrd%rX(icol, irow) = rx(icol)
2237 pgrd%rY(icol, irow) = ry(irow)
2238 enddo
2239 enddo
2240
2241 else
2242
2243 do irow=1,pgrd%iNY
2244 do icol=1,pgrd%iNX
2245 pgrd%rX(icol, irow) = grid_getgridx(pgrd, icol)
2246 pgrd%rY(icol, irow) = grid_getgridy(pgrd, irow)
2247 enddo
2248 enddo
2249
2250 endif
2251
2252 ! print *, 'grid_PopulateXY: obtained grid coordinates, routine ',__FILE__,' line ', __LINE__
2253 ! print *, '(col=1, row=1): ', pGrd%rX(1,1), pGrd%rY(1,1)
2254 ! print *, '(col=1, row=', pGrd%iNY,'): ', pGrd%rX(1, pGrd%iNY),pGrd%rY(1, pGrd%iNY)
2255 ! print *, '(col=',pGrd%iNX,', row=1): ', pGrd%rX(pGrd%iNX, 1),pGrd%rY(pGrd%iNX, 1)
2256 ! print *, '(col=',pGrd%iNX,', row=', pGrd%iNY,'): ', pGrd%rX(pGrd%iNX, pGrd%iNY),pGrd%rY(pGrd%iNX, pGrd%iNY)
2257
2258
2259end subroutine grid_populatexy
2260
2261!----------------------------------------------------------------------
2262
2263subroutine grid_dumpgridextent(pGrd)
2264
2265 type ( general_grid_t ),pointer :: pgrd
2266
2267 if ( associated( pgrd) ) then
2268
2269 call logs%write("---------------------------------------------------")
2270 call logs%write("GRID DETAILS:")
2271 call logs%write("---------------------------------------------------")
2272 call logs%write("file: "//dquote(pgrd%sFilename) )
2273 call logs%write("nx: "//trim( ascharacter(pgrd%iNX) ) )
2274 call logs%write("ny: "//trim( ascharacter(pgrd%iNY) ) )
2275 call logs%write("cellsize: "//trim(ascharacter(pgrd%rGridCellSize) ) )
2276 call logs%write("X0: "//trim(ascharacter(pgrd%rX0) ) )
2277 call logs%write("Y0: "//trim(ascharacter(pgrd%rY0) ) )
2278 call logs%write("X1: "//trim(ascharacter(pgrd%rX1) ) )
2279 call logs%write("Y1: "//trim(ascharacter(pgrd%rY1) ) )
2280 call logs%write("Type: "//trim(ascharacter(pgrd%iDataType) ) )
2281 call logs%write("PROJ4 string: "//dquote(pgrd%sPROJ4_string) )
2282 call logs%write("---------------------------------------------------")
2283
2284 else
2285
2286 call logs%write("---------------------------------------------------")
2287 call logs%write(" ** grid not allocated ** ")
2288 call logs%write("---------------------------------------------------")
2289
2290 endif
2291
2292
2293end subroutine grid_dumpgridextent
2294
2295!----------------------------------------------------------------------
2296function grid_gridtopoint_int(pGrdFrom, pGrdTo, iCol, iRow) result(iValue)
2297
2298 type ( general_grid_t ), pointer :: pgrdfrom
2299 type ( general_grid_t ), pointer :: pgrdto
2300 integer (c_int), intent(in) :: icol
2301 integer (c_int), intent(in) :: irow
2302 integer (c_int) :: ivalue
2303
2304 ! [ LOCALS ]
2305 integer (c_int), dimension(2) :: icolrow
2306 integer (c_int) :: ispread
2307
2308 ! must ensure that there are coordinates associated with the "to" grid...
2309 ! by default, these are left unpopulated during a "normal" swb run
2310 if( .not. allocated(pgrdto%rX) ) call grid_populatexy(pgrdto)
2311 if( .not. allocated(pgrdfrom%rX) ) call grid_populatexy(pgrdfrom)
2312
2313 icolrow = grid_getgridcolrownum(pgrd=pgrdfrom, &
2314 rx=real(pgrdto%rX(icol, irow), c_double), &
2315 ry=real(pgrdto%rY(icol, irow), c_double))
2316
2317 if ( icolrow(column) < 1 .or. icolrow(column) > pgrdfrom%iNX ) &
2318 call die( "Illegal column number supplied: "//trim(ascharacter(icolrow(column))), &
2319 __file__, __line__)
2320
2321 if ( icolrow(row) < 1 .or. icolrow(row) > pgrdfrom%iNY ) &
2322 call die( "Illegal row number supplied: "//trim(ascharacter(icolrow(row))), &
2323 __file__, __line__)
2324
2325 !> @todo check the logic here...intent is to ensure that the majority filter is
2326 !! searching an area that adequately covers corresponding gridcell areas
2327 ispread = max(1, nint(pgrdto%rGridCellSize / pgrdfrom%rGridCellSize / 2.))
2328
2329 ivalue = grid_majorityfilter_int(pgrdfrom=pgrdfrom, &
2330 itargetcol=icolrow(column), &
2331 itargetrow=icolrow(row), &
2332 inodatavalue=pgrdfrom%iNoDataValue, &
2333 ispread=ispread)
2334
2335
2336end function grid_gridtopoint_int
2337
2338
2339!----------------------------------------------------------------------
2340
2341subroutine grid_gridtogrid_int( pGrdFrom, pGrdTo, lUseMajorityFilter )
2342
2343 type ( general_grid_t ), pointer :: pgrdfrom ! pointer to source grid
2344 type ( general_grid_t ), pointer :: pgrdto ! pointer to destination grid
2345 logical (c_bool), intent(in) :: lusemajorityfilter
2346
2347 ! [ LOCALS ]
2348 integer (c_int) :: icol, irow
2349 integer (c_int), dimension(2) :: icolrow
2350 integer (c_int) :: ispread
2351 real (c_float) :: fgridcellratio
2352
2353 ! must ensure that there are coordinates associated with the "to" grid...
2354 ! by default, these are left unpopulated during a "normal" swb run
2355 if(.not. allocated(pgrdto%rX) ) call grid_populatexy(pgrdto)
2356 if(.not. allocated(pgrdfrom%rX) ) call grid_populatexy(pgrdfrom)
2357
2358 fgridcellratio = real(pgrdto%rGridCellSize / pgrdfrom%rGridCellSize, c_float)
2359
2360 ! Allow USER to trigger whether the majority filter is employed or not
2361 if ( lusemajorityfilter ) then
2362
2363 call logs%write( "** Majority filter in use for data from grid file "//dquote(pgrdfrom%sFilename), lecho=true )
2364
2365 ispread = max( 1, nint( fgridcellratio / 2.0_c_float ) )
2366
2367 do irow=1,pgrdto%iNY
2368 do icol=1,pgrdto%iNX
2369
2370 icolrow = grid_getgridcolrownum(pgrd=pgrdfrom, &
2371 rx=real(pgrdto%rX(icol, irow), c_double), &
2372 ry=real(pgrdto%rY(icol, irow), c_double) )
2373
2374 call assert(icolrow(column) > 0 .and. icolrow(column) <= pgrdfrom%iNX, &
2375 "Illegal column number supplied: "//trim(ascharacter(icolrow(column))), &
2376 __file__, __line__)
2377
2378 call assert(icolrow(row) > 0 .and. icolrow(row) <= pgrdfrom%iNY, &
2379 "Illegal row number supplied: "//trim(ascharacter(icolrow(row))), &
2380 __file__, __line__)
2381
2382 pgrdto%iData(icol,irow) = grid_majorityfilter_int( pgrdfrom=pgrdfrom, &
2383 itargetcol=icolrow(column), &
2384 itargetrow=icolrow(row), &
2385 inodatavalue=pgrdfrom%iNoDataValue, &
2386 ispread=ispread)
2387
2388 enddo
2389 enddo
2390
2391 else ! target grid resolution similar to or more dense than source grid resolution: NEAREST NEIGHBOR [DEFAULT]
2392
2393 do irow=1,pgrdto%iNY
2394 do icol=1,pgrdto%iNX
2395
2396 icolrow = grid_getgridcolrownum(pgrd=pgrdfrom, &
2397 rx=real(pgrdto%rX(icol, irow), c_double), &
2398 ry=real(pgrdto%rY(icol, irow), c_double))
2399
2400 call assert(icolrow(column) > 0 .and. icolrow(column) <= pgrdfrom%iNX, &
2401 "Illegal column number supplied: "//trim(ascharacter(icolrow(column))), &
2402 __file__, __line__)
2403
2404 call assert(icolrow(row) > 0 .and. icolrow(row) <= pgrdfrom%iNY, &
2405 "Illegal row number supplied: "//trim(ascharacter(icolrow(row))), &
2406 __file__, __line__)
2407
2408 pgrdto%iData(icol,irow) = pgrdfrom%iData( icolrow(column), icolrow(row) )
2409
2410 enddo
2411 enddo
2412
2413 endif
2414
2415end subroutine grid_gridtogrid_int
2416
2417!----------------------------------------------------------------------
2418
2419subroutine grid_gridtogrid_sgl(pGrdFrom, pGrdTo )
2420
2421 ! [ ARGUMENTS ]
2422 type ( general_grid_t ),pointer :: pgrdfrom ! pointer to source grid
2423 type ( general_grid_t ),pointer :: pgrdto ! pointer to destination grid
2424
2425 ! [ LOCALS ]
2426 integer (c_int), dimension(2) :: icolrow
2427 integer (c_int) :: icol, irow
2428 real (c_float), dimension(3,3) :: rkernel
2429
2430 rkernel = 1.
2431 rkernel(2,2) = 8.
2432
2433 ! must ensure that there are coordinates associated with the "to" grid...
2434 ! by default, these are left unpopulated during a "normal" swb run
2435 if(.not. allocated(pgrdto%rX) ) call grid_populatexy(pgrdto)
2436 if(.not. allocated(pgrdfrom%rX) ) call grid_populatexy(pgrdfrom)
2437
2438!!! *$OMP PARALLEL DO ORDERED PRIVATE(iRow, iCol, iColRow)
2439
2440 do irow=1,pgrdto%iNY
2441 do icol=1,pgrdto%iNX
2442
2443 icolrow = grid_getgridcolrownum(pgrd=pgrdfrom, &
2444 rx=real(pgrdto%rX(icol, irow), c_double), &
2445 ry=real(pgrdto%rY(icol, irow), c_double))
2446
2447 if ( icolrow(column) < 1 .or. icolrow(column) > pgrdfrom%iNX) cycle
2448 if ( icolrow(row) < 1 .or. icolrow(row) > pgrdfrom%iNY) cycle
2449
2450 pgrdto%rData(icol,irow) = pgrdfrom%rData( icolrow(column), icolrow(row) )
2451
2452
2453 ! printout = (iRow==1 .and. iCol==1) .or. (iRow==1 .and. iCol==pGrdTo%iNX) &
2454 ! .or. (iRow==pGrdTo%iNY .and. iCol==1) .or. (iRow==pGrdTo%iNY .and. iCol==pGrdTo%iNX)
2455
2456 ! if (printout) then
2457 ! print *, trim(__FILE__), ': ', __LINE__
2458 ! print *, 'getting data from the file: ', trim(pGrdFrom%sFilename)
2459 ! print *, ' coord-to : ', pGrdTo%rX(iCol, iRow), pGrdTo%rY(iCol, iRow)
2460 ! print *, ' to: (iCol, iRow): ', iCol, iRow
2461 ! print *, ' from: (col, row): ', iColRow(COLUMN), iColRow(ROW)
2462 ! endif
2463
2464 enddo
2465 enddo
2466
2467!!! *$OMP END PARALLEL DO
2468
2469end subroutine grid_gridtogrid_sgl
2470
2471!----------------------------------------------------------------------
2472
2473function grid_majorityfilter_int(pGrdFrom, iTargetCol, &
2474 iTargetRow, iNoDataValue, iSpread) result(iRetVal)
2475
2476 type ( general_grid_t ),pointer :: pgrdfrom ! pointer to source grid
2477 integer (c_int) :: itargetcol ! column number of target cell
2478 integer (c_int) :: itargetrow ! row number of target cell
2479 integer (c_int) :: inodatavalue
2480 integer (c_int) :: iretval
2481 integer (c_int) :: ispread ! integer representing how many
2482 ! cells should be searched
2483
2484 ! [ LOCALS ]
2485 integer (c_int), dimension(625) :: ivalue
2486 integer (c_int), dimension(625) :: icount
2487 logical (c_bool) :: lmatch
2488 integer (c_int) :: irow
2489 integer (c_int) :: icol
2490 integer (c_int) :: i, ilast, iindex, icellnum
2491
2492 ivalue = 0
2493 icount = 0
2494 ilast = 0
2495 icellnum = 0
2496
2497 ! need to set a reasonable upper bound on the spread to consider
2498 ispread = min(ispread, 12)
2499
2500 do irow=max(1,itargetrow - ispread), min(pgrdfrom%iNY, itargetrow + ispread)
2501 do icol=max(1,itargetcol - ispread), min(pgrdfrom%iNX, itargetcol + ispread)
2502!
2503 icellnum = icellnum + 1
2504!
2505 lmatch = false
2506!
2507 do i=1, ilast
2508 if (ivalue(i) == pgrdfrom%iData(icol, irow) ) then
2509 icount(i) = icount(i) + 1
2510 lmatch = true
2511 exit
2512 endif
2513 enddo
2514
2515 ! this integer value has not been found previously; add the value to the next position
2516 ! and increment counter
2517 if (.not. lmatch) then
2518 ilast = ilast + 1
2519 ivalue(ilast) = pgrdfrom%iData(icol, irow)
2520 icount(ilast) = icount(ilast) + 1
2521 lmatch = true
2522 endif
2523
2524 enddo
2525 enddo
2526
2527 iindex = maxloc(icount, dim=1)
2528 iretval = ivalue(iindex)
2529
2530end function grid_majorityfilter_int
2531
2532function grid_convolve_sgl( pGrdFrom, iTargetCol, &
2533 iTargetRow, rKernel, rNoDataValue) result(rRetVal)
2534
2535 type ( general_grid_t ),pointer :: pgrdfrom ! pointer to source grid
2536 integer (c_int) :: inx
2537 integer (c_int) :: iny
2538 integer (c_int) :: itargetcol ! column number of target cell
2539 integer (c_int) :: itargetrow ! row number of target cell
2540 real (c_float), dimension(:,:) :: rkernel
2541 real (c_float), optional :: rnodatavalue
2542 real (c_float) :: rretval
2543
2544 ! [ LOCALS ]
2545 integer (c_int) :: irowmin, irowmax
2546 integer (c_int) :: icolmin, icolmax
2547 integer (c_int) :: ikernelsize ! i.e. 3, 5, 7, 9, 11, etc.
2548 integer (c_int) :: iincvalue
2549 integer (c_int) :: icol, irow, irownum, icolnum
2550 real (c_float) :: rkernelsum
2551
2552 rkernelsum = rzero
2553 rretval = rzero
2554
2555 ikernelsize = size(rkernel, dim=1)
2556 iincvalue = (ikernelsize - 1) / 2
2557
2558 inx = ubound(pgrdfrom%rData,1)
2559 iny = ubound(pgrdfrom%rData,2)
2560
2561 irowmin = max(1,itargetrow - iincvalue)
2562 irowmax = min(iny, itargetrow + iincvalue)
2563
2564 icolmin = max(1,itargetcol - iincvalue)
2565 icolmax = min(iny, itargetcol + iincvalue)
2566
2567 if( (irowmax - irowmin + 1) /= ikernelsize &
2568 .or. (icolmax - icolmin + 1) /= ikernelsize ) then
2569
2570 ! This is a simple, but less than desirable way to treat cells that
2571 ! fall near the edge of the source grid. Ideally, the source grid will
2572 ! be greater than the target grid by a wide enough margin that this
2573 ! code will rarely be used.
2574 rretval = pgrdfrom%rData(itargetcol, itargetrow)
2575
2576 else
2577
2578 do icol=0,ikernelsize-1
2579 do irow=0,ikernelsize-1
2580
2581 irownum = irow + irowmin
2582 icolnum = icol + icolmin
2583
2584 if (irownum > iny .or. icolnum > inx) then
2585 cycle ! our calculated row or column number is outside of the
2586 ! bounds of the rValues array
2587 else
2588 rretval = rretval + pgrdfrom%rData(icolnum,irownum) * rkernel(icol+1, irow+1)
2589 rkernelsum = rkernelsum + rkernel(icol+1, irow+1)
2590 endif
2591 enddo
2592 enddo
2593
2594 if (rkernelsum > 0.) rretval = rretval / rkernelsum
2595
2596 endif
2597
2598end function grid_convolve_sgl
2599
2600end module grid
constants_and_conversions
This module contains physical constants and convenience functions aimed at performing unit conversion...
Definition constants_and_conversions.F90:8
constants_and_conversions::true
logical(c_bool), parameter, public true
Definition constants_and_conversions.F90:38
constants_and_conversions::rzero
real(c_float), parameter, public rzero
Definition constants_and_conversions.F90:54
constants_and_conversions::datatype_real
integer(c_int), parameter datatype_real
Definition constants_and_conversions.F90:17
constants_and_conversions::false
logical(c_bool), parameter, public false
Definition constants_and_conversions.F90:39
constants_and_conversions::degrees_to_radians
real(c_double), parameter, public degrees_to_radians
Definition constants_and_conversions.F90:30
constants_and_conversions::rbigval
real(c_float), parameter, public rbigval
Definition constants_and_conversions.F90:42
constants_and_conversions::datatype_int
integer(c_int), parameter datatype_int
Definition constants_and_conversions.F90:15
constants_and_conversions::radians_to_degrees
real(c_double), parameter, public radians_to_degrees
Definition constants_and_conversions.F90:31
constants_and_conversions::datatype_double
integer(c_int), parameter datatype_double
Definition constants_and_conversions.F90:18
exceptions::die
subroutine, public die(smessage, smodule, iline, shints, scalledby, icalledbyline)
Definition exceptions.F90:140
grid
Provides support for input and output of gridded ASCII data, as well as for creation and destruction ...
Definition grid.F90:8
grid::grid_createcomplete
type(general_grid_t) function, pointer, public grid_createcomplete(inx, iny, rx0, ry0, rx1, ry1, idatatype)
Creates a grid of a specified type.
Definition grid.F90:211
grid::grid_getgridcolnum
integer(c_int) function, public grid_getgridcolnum(pgrd, rx)
Definition grid.F90:1999
grid::grid_coordinatesfallinsidegrid
logical(c_bool) function, public grid_coordinatesfallinsidegrid(pgrd, rx, ry)
Definition grid.F90:1971
grid::terrormessage
type(error_message_t), dimension(49) terrormessage
Definition grid.F90:99
grid::grid_populatexy
subroutine, public grid_populatexy(pgrd, rx, ry)
Definition grid.F90:2195
grid::grid_nodata_real
real(c_float), parameter, public grid_nodata_real
Definition grid.F90:21
grid::grid_searchcolumn
real(c_float) function grid_searchcolumn(pgrd, rxval, rzval, rnodata)
Definition grid.F90:1819
grid::output_grid_directory_name
character(len=256) output_grid_directory_name
Definition grid.F90:176
grid::grid_completelycover
logical(c_bool) function, public grid_completelycover(pbasegrd, pothergrd, rtolerance)
Definition grid.F90:1328
grid::eight_cells
integer(c_int), parameter eight_cells
Definition grid.F90:168
grid::grid_active_cell
integer(c_int), parameter, public grid_active_cell
Definition grid.F90:32
grid::grid_read
type(general_grid_t) function, pointer, public grid_read(sfilename, sfiletype, idatatype)
Definition grid.F90:447
grid::grid_lookupreal
real(c_float) function grid_lookupreal(pgrd, rxval, ryval)
Definition grid.F90:1953
grid::four_cells
integer(c_int), parameter four_cells
Definition grid.F90:167
grid::grid_convolve_sgl
real(c_float) function grid_convolve_sgl(pgrdfrom, itargetcol, itargetrow, rkernel, rnodatavalue)
Definition grid.F90:2534
grid::grid_writegrid
subroutine, public grid_writegrid(sfilename, pgrd, ioutputformat)
Definition grid.F90:1033
grid::nc_fill_float
real(c_float), parameter nc_fill_float
Definition grid.F90:34
grid::grid_dumpgridextent
subroutine, public grid_dumpgridextent(pgrd)
Definition grid.F90:2264
grid::grid_getgridrownum
integer(c_int) function, public grid_getgridrownum(pgrd, ry)
Definition grid.F90:2032
grid::lu_grid
integer(c_int) lu_grid
Definition grid.F90:174
grid::grid_readarcgrid_fn
type(general_grid_t) function, pointer grid_readarcgrid_fn(sfilename, idatatype)
Definition grid.F90:516
grid::grid_createsimple
type(general_grid_t) function, pointer, public grid_createsimple(inx, iny, rx0, ry0, rgridcellsize, idatatype)
Definition grid.F90:281
grid::grid_datatype_int
integer(c_int), parameter, public grid_datatype_int
Definition grid.F90:25
grid::grid_destroy
subroutine, public grid_destroy(pgrd)
Definition grid.F90:366
grid::grid_writearcgrid
subroutine, public grid_writearcgrid(sfilename, pgrd)
Definition grid.F90:1054
grid::row
integer(c_int), parameter, private row
Definition grid.F90:171
grid::grid_lookuprow
subroutine grid_lookuprow(pgrd, ryval, ibefore, iafter, rfrac)
Definition grid.F90:1493
grid::lu_temp
integer(c_int) lu_temp
Definition grid.F90:174
grid::grid_nodata_double
real(c_double), parameter, public grid_nodata_double
Definition grid.F90:22
grid::grid_set_nodata_value
subroutine, public grid_set_nodata_value(pgrd, ivalue, fvalue)
Definition grid.F90:2160
grid::grid_lookupcolumn
subroutine grid_lookupcolumn(pgrd, rxval, ibefore, iafter, rfrac)
Definition grid.F90:1425
grid::grid_gridtopoint_int
integer(c_int) function grid_gridtopoint_int(pgrdfrom, pgrdto, icol, irow)
Definition grid.F90:2297
grid::grid_gridtogrid_int
subroutine, public grid_gridtogrid_int(pgrdfrom, pgrdto, lusemajorityfilter)
Definition grid.F90:2342
grid::grid_getgridcolrownum
integer(c_int) function, dimension(2) grid_getgridcolrownum(pgrd, rx, ry)
Definition grid.F90:2053
grid::grid_rowcolfallsinsidegrid
logical(c_bool) function, public grid_rowcolfallsinsidegrid(pgrd, irow, icol)
Definition grid.F90:1985
grid::grid_nodata_int
integer(c_int), parameter, public grid_nodata_int
Definition grid.F90:20
grid::grid_datatype_real
integer(c_int), parameter, public grid_datatype_real
Definition grid.F90:26
grid::grid_readexisting
subroutine, public grid_readexisting(sfilename, sfiletype, pgrd)
Definition grid.F90:472
grid::output_arc
integer(c_int), parameter, public output_arc
Definition grid.F90:30
grid::grid_interpolate
real(c_float) function, public grid_interpolate(pgrd, rxval, ryval)
Definition grid.F90:1729
grid::column
integer(c_int), parameter, private column
Definition grid.F90:170
grid::output_surfer
integer(c_int), parameter, public output_surfer
Definition grid.F90:29
grid::grid_readsurfergrid_sub
subroutine grid_readsurfergrid_sub(sfilename, pgrd)
Definition grid.F90:955
grid::nc_fill_int
integer(c_int), parameter nc_fill_int
Definition grid.F90:33
grid::nc_fill_double
real(c_double), parameter nc_fill_double
Definition grid.F90:35
grid::grid_gridtogrid_sgl
subroutine, public grid_gridtogrid_sgl(pgrdfrom, pgrdto)
Definition grid.F90:2420
grid::grid_getgridy
real(c_double) function, public grid_getgridy(pgrd, irownumber)
Definition grid.F90:2182
grid::grid_writesurfergrid
subroutine, public grid_writesurfergrid(sfilename, pgrd)
Definition grid.F90:1144
grid::grid_checkintegergridvalues
subroutine grid_checkintegergridvalues(pgrd, sfilename)
Definition grid.F90:1617
grid::grid_conform
logical(c_bool) function, public grid_conform(pgrd1, pgrd2, rtolerance)
Definition grid.F90:1264
grid::grid_checkforproj4error
subroutine, public grid_checkforproj4error(iretval, sfromproj4, stoproj4)
Definition grid.F90:1666
grid::grid_set_output_directory_name
subroutine, public grid_set_output_directory_name(sdirname)
Definition grid.F90:182
grid::grid_transform
subroutine, public grid_transform(pgrd, sfromproj4, stoproj4, rx, ry)
Call PROJ4 to transform coordinates.
Definition grid.F90:1538
grid::grid_getgridx
real(c_double) function, public grid_getgridx(pgrd, icolumnnumber)
Definition grid.F90:2170
grid::grid_readsurfergrid_fn
type(general_grid_t) function, pointer grid_readsurfergrid_fn(sfilename, idatatype)
Definition grid.F90:873
grid::grid_datatype_double
integer(c_int), parameter, public grid_datatype_double
Definition grid.F90:27
grid::grid_readarcgrid_sub
subroutine grid_readarcgrid_sub(sfilename, pgrd)
Definition grid.F90:685
grid::grid_majorityfilter_int
integer(c_int) function grid_majorityfilter_int(pgrdfrom, itargetcol, itargetrow, inodatavalue, ispread)
Definition grid.F90:2475
logfiles::logs
type(logfile_t), public logs
Definition logfiles.F90:62
grid::pj_init_and_transform
interface to C code that provides a simplified entry point to PROJ4 capabilities: it has been modifie...
Definition grid.F90:42