 |
forums.silverfrost.com
Welcome to the Silverfrost forums
|
| View previous topic :: View next topic |
| Author |
Message |
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
 Posted: Tue Jan 21, 2025 4:45 pm Post subject: 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:
| Code: | 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:
| Code: | 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 |
|
| Back to top |
|
 |
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Tue Jan 21, 2025 6:43 pm Post subject: |
|
|
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
| Code: | 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. |
|
| Back to top |
|
|
Kenneth_Smith
Joined: 18 May 2012
Posts: 801
Location: Hamilton, Lanarkshire, Scotland.
|
| Posted: Tue Jan 21, 2025 8:55 pm Post subject: |
|
|
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.
| Code: | 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 |
|
|
| Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006
Posts: 2615
Location: Sydney
|
| Posted: Wed Jan 22, 2025 1:06 am Post subject: |
|
|
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 ?
| Code: | ! 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. |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 146
Location: United Kingdom
|
| Posted: Thu Jan 23, 2025 12:12 pm Post subject: |
|
|
| Thanks for the pointers |
|
| Back to top |
|
|
|
|
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|
FAQ
Search
Memberlist
Usergroups
Register
Profile
Log in to check your private messages
Log in
