Reading an ESRI ascii raster file

Hello

I am trying to set a code to read the grid data from an ESRI ascii raster file and output the data as X,Y,Z. The file structure is simple, but not found an easy solution to extract the header variable names and associated values on read.

This snippet works extracting the header but rather convoluted, used real*8 but real(kind=2) is also valid:

      PROGRAM ReadHeader

      IMPLICIT NONE

      ! Define the structure for the grid header
      type ascii_raster
         integer  :: ncols
         integer  :: nrows
         REAL*8  :: xllcorner
         REAL*8  :: yllcorner
         REAL*8  :: cellsize
         REAL*8  :: NODATA_value
      END type ascii_raster

      TYPE(ascii_raster) :: grid
      CHARACTER :: label 

      ! Open the file for reading
      OPEN(UNIT=10, FILE='header2.asc', STATUS='OLD')

      ! Read data from the file into the structure fields
      READ(10,*) label, grid%ncols
      READ(10,*) label, grid%nrows
      READ(10,*) label, grid%xllcorner
      READ(10,*) label, grid%yllcorner
      READ(10,*) label, grid%cellsize
      READ(10,*) label, grid%NODATA_value

      ! Close the file after reading
      CLOSE(10)

      ! Print the values read to verify
      PRINT *, 'ncols = ', grid%ncols
      PRINT *, 'nrows = ', grid%nrows
      PRINT *, 'xllcorner = ', grid%xllcorner
      PRINT *, 'yllcorner = ', grid%yllcorner
      PRINT *, 'cellsize = ', grid%cellsize
      PRINT *, 'NODATA_value = ', grid%NODATA_value

      END PROGRAM ReadHeader

The input file is small as a tester:

ncols 10
nrows 10
xllcorner 563435.00
yllcorner 5258305.00
cellsize 10.00
NODATA_value -9999.00
103 102 102 101 100 99 97 95 94 92
100 100 100 99 97 96 94 92 90 88
97 97 96 95 93 92 91 88 85 82
94 93 92 90 89 87 84 82 79 76
90 89 87 86 84 81 78 75 74 72
86 84 83 80 77 75 72 70 69 69
81 79 76 74 72 69 68 67 68 71
76 73 71 69 68 66 67 68 71 73
71 69 68 66 65 66 68 71 74 76
68 66 65 65 66 68 71 74 77 81

Is there a simpler methodology to extract the header variable names with their values and the zdata(ncols,nrows)?

Given the grid dimensions (ncols, nrows) is it an easy task to get the data into X,Y, Z ?

Ideally I want to build this into a subroutine where the only input is the data file, eg asc_2_xyz(filename, X, Y, Z) .

Thanks for any pointers

Lester

This works to get the header data and z-grid, but need to know the positions in the file where numerical data start in the header

program testreadesriasciiraster
  implicit none

  ! variables
  character(len=256) :: filename
  integer :: ncols, nrows
  real(kind=2) :: xllcorner, yllcorner, cellsize, nodata_value
  real(kind=2), allocatable :: zdata(:,:)
  integer :: i, j

  ! set the filename to the test input file
  filename = 'header2.asc'

  ! call the subroutine to read the file
  call readesriasciiraster(filename, ncols, nrows, xllcorner, yllcorner, cellsize, nodata_value, zdata)

  ! print the header values
  print *, 'ncols = ', ncols
  print *, 'nrows = ', nrows
  print *, 'xllcorner = ', xllcorner
  print *, 'yllcorner = ', yllcorner
  print *, 'cellsize = ', cellsize
  print *, 'nodata_value = ', nodata_value

  xmax = xllcorner + cellsize * ncols
  ymax = yllcorner + cellsize * nrows 

  ! print the z-data
  print *, 'z-data:'
  do i = 1, nrows
    do j = 1, ncols
      write(*, '(f8.2)', advance='no') zdata(i, j)
      if (j < ncols) write(*, '(a)', advance='no') ' '
    end do
    print *
  end do

  ! deallocate the array
  deallocate(zdata)

print *, xmax, ymax

contains

subroutine readesriasciiraster(filename, ncols, nrows, xllcorner, yllcorner, cellsize, nodata_value, zdata)
  implicit none

  ! input arguments
  character(len=*), intent(in) :: filename

  ! output arguments
  integer, intent(out) :: ncols, nrows
  real*8, intent(out) :: xllcorner, yllcorner, cellsize, nodata_value
  real*8, allocatable, intent(out) :: zdata(:,:)

  ! local variables
  character(len=256) :: line
  integer :: i, j, ios
  integer :: unit_number

  ! open the file
  unit_number = 10
  open(unit_number, file=filename, status='old', action='read', iostat=ios)
  if (ios /= 0) then
    print *, 'error opening file: ', filename
    stop
  end if

  ! read and parse the header values
  read(unit_number, '(a)', iostat=ios) line
  if (ios /= 0) stop 'error reading ncols line'
  read(line(7:), *) ncols

  read(unit_number, '(a)', iostat=ios) line
  if (ios /= 0) stop 'error reading nrows line'
  read(line(7:), *) nrows

  read(unit_number, '(a)', iostat=ios) line
  if (ios /= 0) stop 'error reading xllcorner line'
  read(line(11:), *) xllcorner

  read(unit_number, '(a)', iostat=ios) line
  if (ios /= 0) stop 'error reading yllcorner line'
  read(line(11:), *) yllcorner

  read(unit_number, '(a)', iostat=ios) line
  if (ios /= 0) stop 'error reading cellsize line'
  read(line(10:), *) cellsize

  read(unit_number, '(a)', iostat=ios) line
  if (ios /= 0) stop 'error reading nodata_value line'
  read(line(14:), *) nodata_value



  ! allocate the z-data array
  allocate(zdata(nrows, ncols))

  ! read the z-data values
  do i = 1, nrows
      read(unit_number, *, iostat=ios) (zdata(i, j), j=1, ncols)
      if (ios /= 0) then
        print *, 'error reading z-data at (', i, ',', j, ')'
        stop
      end if
  end do

  ! close the file
  close(unit_number)

end subroutine readesriasciiraster

end program testreadesriasciiraster

xllcorner and yllcorner correspond to the lower-left grid point ending top-right.

If I have understood the format of the provided sample file correctly this should give you the information you need in the file output.txt.

    program esri_to_xyz
    implicit none
    integer :: ncols, nrows,  i, j
    real    :: xllcorner, yllcorner, cellsize, nodata_value
    real, allocatable :: dataval(:,:)
    character(len=20) :: label 
    character(len=30), parameter :: fmt1 = '(F12.2, 1X, F12.2, 1X, F12.2)'

    open(10, file='input.txt', status='old', action='read')

    read(10, *) label, ncols
    read(10, *) label, nrows
    read(10, *) label, xllcorner
    read(10, *) label, yllcorner
    read(10, *) label, cellsize
    read(10, *) label, nodata_value

    allocate(dataval(nrows, ncols))

    read(10, *) ((dataval(i, j), j = 1, ncols), i = 1, nrows)
    close(10,status='keep')

    open(10, file='output.txt', status='replace', action='write')
    
    do i = 1, nrows, 1
      do j = 1, ncols, 1
        if (dataval(i, j) .ne. nodata_value) then
          write(10, fmt1) xllcorner + (j-1) * cellsize, yllcorner + (nrows-i) * cellsize, dataval(i, j)
        end if
      end do
    end do

    close(10, status='keep')

    deallocate(dataval)

end program esri_to_xyz

The following is an incomplete template of reading the (6) headers with some flexibility and validation of the labels. It can be help to provide some basic checking of names, although they can vary with read data ?

 !  should be enhanced with 'iostat=' !!
  module grid_data
     integer, parameter :: max_header = 10

     integer :: ncols        = 1
     integer :: nrows        = 1
     real*8  :: xllcorner    = 0.
     real*8  :: yllcorner    = 0.
     real*8  :: cellsize     = 1.
     real*8  :: NODATA_value = -999.

     real*8, allocatable :: grid_values(:,:)

   end module grid_data

     open ( unit=10, file='grid_info.txt' )
     call read_header
     call read_grid
   end

   subroutine read_header
     use grid_data
     integer :: i
     real*8  :: value
     character :: line*80, label*12

     do i = 1, max_header
       read (10, fmt='(a)' ) line
       read ( line,* ) label, value

       select case ( label )
         case ( 'ncols' )
           ncols = nint ( value )
         case ( 'nrows' )
           nrows = nint ( value )
         case ( 'xllcorner' )
           xllcorner = value
         case ( 'yllcorner' )
           yllcorner = value
         case ( 'cellsize' )
           cellsize = value
         case ( 'NODATA_value' )
           NODATA_value = value
           exit                         ! assume this will be the last

         case default
           write (*,*) i, ' label :', label,' is not recognised'
           stop

       end select

     end do

   ! Print the values read to verify
      PRINT *, 'ncols = ', ncols
      PRINT *, 'nrows = ', nrows
      PRINT *, 'xllcorner = ', xllcorner
      PRINT *, 'yllcorner = ', yllcorner
      PRINT *, 'cellsize = ', cellsize
      PRINT *, 'NODATA_value = ', NODATA_value

   end subroutine read_header
     
   subroutine read_grid
     use grid_data
     integer :: i
     real*8, allocatable :: values(:)

     allocate ( values(ncols) )
     allocate ( grid_values(ncols,nrows) )
     
     do i = 1, nrows
       read (10,*) values
       grid_values(:,i) = values
       write (*,fmt='(99(1x,f0.1))') values
     end do
   end subroutine read_grid

I find that it is good to have some validation of the header info, especially for when the header info changes and you did not expect the change. You could also include some basic check of value range. There will probably be a data file that does not fit the mold.

It would be better if there was a header to identify the start of row values in the data structure.

Thanks for the pointers