Surfer grd文件读写Fortran代码(面向对象编程示例)

地球屋里老师 2021-08-21 14:49:22

    module GlodenGrd 用于读写Golden Surfer采用的网格文件(*.grd)，其中定义了一个父类(超类) grd及三个子类 grd_DSAA，grd_DSBB， grd_DSRB (三个子类对应于三种不同的数据存储格式，详见Grd文件格式说明)，并为每个类提供了两个绑定的过程(read: 用于读取文件；convert: 文件类型转换)。为便于管理代码，子类的过程存放于对应的子模块中：submodule(GlodenGrd) GlodenGrd_DSAA，submodule(GlodenGrd) GlodenGrd_DSBB，submodule(GlodenGrd) GlodenGrd_DSRB。

    该代码采用面向对象(OOP)设计思想，实现了读取三种格式的通用算法，并可在任意两种格式之间进行转换。学习该代码可初探Fortran的OOP特性。用户亦可根据自身需求添加功能。

module GlodenGrd
!****************************************************************************
! 处理Grd文件
! author: LI xingwang, Chang''an university, China.
! email : lixingwang@chd.edu.cn
!****************************************************************************
implicit none
private
public read_grd
real(4), parameter:: noData = 1.71041e38
type tag
  character(4):: sectionName = ''
  integer(4):: sectionSize = 0
end type
!超类，抽象类
type, abstract, public:: grd
contains
procedure(sub_read), private, deferred:: read
procedure(sub_convert), public, deferred:: convert
end type
abstract interface
subroutine sub_read(this, filename)
import
implicit none
class(grd), intent(out):: this
character(*), intent(in):: filename
end subroutine
subroutine sub_convert(this, filename, fileFormat)
import
implicit none
class(grd), target, intent(in):: this
character(*), intent(in):: filename
character(4), optional, intent(in):: fileFormat
end subroutine
end interface
! Surfer 6 Text
type, public, extends(grd):: grd_DSAA
  character(4):: id = 'DSAA'
  integer nx, ny
  real xlo, xhi, ylo, yhi, zlo, zhi
  real, allocatable:: z(:,:)
contains
procedure, pass:: read => read_grd_DSAA
procedure, pass:: convert => convert_grd_DSAA
end type
! Surfer 6 Binary
type, public, extends(grd):: grd_DSBB
  character(4):: id = 'DSBB'
  integer(2) nx, ny
  real(8) xlo, xhi, ylo, yhi, zlo, zhi
  real(4), allocatable:: z(:,:)
contains
procedure, pass:: read => read_grd_DSBB
procedure, pass:: convert => convert_grd_DSBB
end type
! Surfer 7
type, public, extends(grd):: grd_DSRB
  type(tag):: tHeader=tag('DSRB',4), tGrid=tag('GRID',72), tData=tag('DATA',0), tFlti=tag('FLTI',0)
  character(4):: id = 'DSRB'
  integer(4):: version=2
  integer(4) ny, nx
  real(8):: xlo, ylo, xSize, ySize, zlo, zhi, rotation=0.0, blankValue=noData
  real(8), allocatable:: z(:,:)
  integer(4):: nTraces=0, nVertices=0
  integer(4), allocatable:: traceIndex(:,:)  ! 2 × nTraces
  real(8), allocatable:: verticeLoc(:,:)  ! 2 × nVertices
contains
procedure, pass:: read => read_grd_DSRB
procedure, pass:: convert => convert_grd_DSRB
end type

!接口, submodule GlodenGrd_DSAA
interface
module subroutine read_grd_DSAA(this, filename)
implicit none
class(grd_DSAA), intent(out):: this
character(*), intent(in):: filename
end subroutine
module subroutine convert_grd_DSAA(this, filename, fileFormat)
implicit none
class(grd_DSAA), target, intent(in):: this
character(*), intent(in):: filename
character(4), optional, intent(in):: fileFormat
end subroutine
end interface
!接口, submodule GlodenGrd_DSBB
interface
module subroutine read_grd_DSBB(this, filename)
implicit none
class(grd_DSBB), intent(out):: this
character(*), intent(in):: filename
end subroutine
module subroutine convert_grd_DSBB(this, filename, fileFormat)
implicit none
class(grd_DSBB), target, intent(in):: this
character(*), intent(in):: filename
character(4), optional, intent(in):: fileFormat
end subroutine
end interface
!接口, submodule GlodenGrd_DSRB
interface
module subroutine read_grd_DSRB(this, filename)
implicit none
class(grd_DSRB), intent(out):: this
character(*), intent(in):: filename
end subroutine
module subroutine convert_grd_DSRB(this, filename, fileFormat)
implicit none
class(grd_DSRB), target, intent(in):: this
character(*), intent(in):: filename
character(4), optional, intent(in):: fileFormat
end subroutine
end interface
contains
!****************************************************************************
! 读取 grd 格式文件
! 输入:
! filename     ... 文件名
! 输出：
! this         ... 数据文件
!****************************************************************************
subroutine read_grd(this, filename)
implicit none
class(grd), pointer:: this
character(*), intent(in):: filename
character(4) fileFormat
type(grd_DSAA), pointer:: temp1
type(grd_DSBB), pointer:: temp2
type(grd_DSRB), pointer:: temp3
integer unt
open(newunit=unt,file=trim(filename),status='old',access='stream')
read(unt) fileFormat
close(unt)
select case(fileFormat)
case('DSAA')
  allocate(temp1); this=>temp1
case('DSBB')
  allocate(temp2); this=>temp2
case('DSRB')
  allocate(temp3); this=>temp3
end select
call this%read(filename)
end subroutine

!****************************************************************************
! 将 grd 转为指定格式
! 输入:
! this         ... 数据
! filename     ... 输出文件名
! fileFormat   ... 输出文件的格式, 忽略该参数按原格式输出
!****************************************************************************
subroutine convert_grd(this, filename, fileFormat)
implicit none
class(grd), intent(in):: this
character(*), intent(in):: filename
character(4), optional,intent(in):: fileFormat
character(4) fmt
integer(4) nx, ny, version, nTraces, nVertices
real(8) xmin, xmax, ymin, ymax, zmin, zmax, rotation, blankValue
real(8), allocatable:: z(:,:)
integer(4), allocatable:: traceIndex(:,:)  ! 2 × nTraces
real(8), allocatable:: verticeLoc(:,:)  ! 2 × nVertices
integer i, unt

version = 2; rotation = 0.0; blankValue = noData
nTraces=0; nVertices=0
!转为大写字母
if(present(fileFormat)) then
  do concurrent (i=1:len(fileFormat))
    if(fileFormat(i:i)>='a'.and.fileFormat(i:i)<='z') then
      fmt(i:i) = achar( iachar(fileFormat(i:i))-32 )
    else
      fmt(i:i) = fileFormat(i:i)
    end if
  end do
end if
!提取数据
select type(this)
type is(grd_DSAA)
  if(.not.present(fileFormat)) fmt = 'DSAA'
  nx = this%nx;    ny = this%ny
  xmin = this%xlo; xmax = this%xhi
  ymin = this%ylo; ymax = this%yhi
  allocate(z(nx,ny)); z = this%z
  zmin = minval(z); zmax = maxval(z)
type is(grd_DSBB)
  if(.not.present(fileFormat)) fmt = 'DSBB'
  nx = this%nx;    ny = this%ny
  xmin = this%xlo; xmax = this%xhi
  ymin = this%ylo; ymax = this%yhi
  allocate(z(nx,ny)); z = this%z
  zmin = minval(z); zmax = maxval(z)
type is(grd_DSRB)
  if(.not.present(fileFormat)) fmt = 'DSRB'
  nx = this%nx;    ny = this%ny
  xmin = this%xlo; xmax = this%xlo+(this%nx-1)*this%xSize
  ymin = this%ylo; ymax = this%ylo+(this%ny-1)*this%ySize
  allocate(z(nx,ny)); z = this%z
  zmin = minval(z); zmax = maxval(z)
  version = this%version
  rotation = this%rotation
  blankValue = this%blankValue
  nTraces = this%nTraces
  nVertices = this%nVertices
  if(nTraces>0.and.nVertices>0) then
    allocate(traceIndex,source=this%traceIndex)
    allocate(verticeLoc,source=this%verticeLoc)
  else
    nTraces = 0; nVertices = 0
  end if
end select
!输出文件
select case(fmt)
case('DSAA')
  open(newunit=unt,file=filename)
  write(unt,"('DSAA')")
  write(unt,"(g0,2x,g0)") nx, ny
  write(unt,"(g0,2x,g0)") xmin, xmax
  write(unt,"(g0,2x,g0)") ymin, ymax
  write(unt,"(g0,2x,g0)") zmin, zmax
  do i = 1, ny
    write(unt,"(*(g0,2x))") z(:,i)
  end do
  close(unt)
case('DSBB')
  open(newunit=unt,file=filename,access='stream')
  write(unt)'DSBB'
  write(unt) int(nx,2), int(ny,2)
  write(unt) xmin, xmax
  write(unt) ymin, ymax
  write(unt) zmin, zmax
  do i = 1, ny
    write(unt) real(z(:,i),4)
  end do
  close(unt)
case('DSRB')
  open(newunit=unt,file=filename,access='stream')
  write(unt) 'DSRB', 4_4, version  ! Header
  write(unt) 'GRID', 72_4, ny, nx  ! GRID
  write(unt) xmin, ymin
  write(unt) (xmax-xmin)/(nx-1), (ymax-ymin)/(ny-1)
  write(unt) zmin, zmax
  write(unt) rotation, blankValue
  write(unt) 'DATA', size(z)*8, z  ! DATA
  if(nTraces>0.and.nVertices>0) then
    write(unt) 'FLTI', 8+4*size(traceIndex)+8*size(verticeLoc) !FLTI
    write(unt) nTraces, nVertices, traceIndex, verticeLoc
  end if
  close(unt)
end select
end subroutine
end module GlodenGrd

submodule(GlodenGrd) GlodenGrd_DSAA
implicit none
contains
!****************************************************************************
! 读取 DSAA 格式文件
! 输入:
! filename     ... 文件名
! 输出：
! this         ... 数据
!****************************************************************************
module subroutine read_grd_DSAA(this, filename)
implicit none
class(grd_DSAA), intent(out):: this
character(*), intent(in):: filename
integer unt
open(newunit=unt,file=trim(filename),status='old')
read(unt,*) this%id
read(unt,*) this%nx, this%ny
read(unt,*) this%xlo, this%xhi
read(unt,*) this%ylo, this%yhi
read(unt,*) this%zlo, this%zhi
if(allocated(this%z)) deallocate(this%z)
allocate(this%z(this%nx,this%ny))
read(unt,*) this%z
close(unt)
this%zlo = minval(this%z)
this%zhi = maxval(this%z)
end subroutine

!****************************************************************************
! 将 DSAA 转为指定格式
! 输入:
! this         ... 数据
! filename     ... 输出文件名
! fileFormat   ... 输出文件的格式
!****************************************************************************
module subroutine convert_grd_DSAA(this, filename, fileFormat)
implicit none
class(grd_DSAA), target, intent(in):: this
character(*), intent(in):: filename
character(4), optional, intent(in):: fileFormat
class(grd), pointer:: pGrd
pGrd => this
if(present(fileFormat)) then
  call convert_grd(pGrd, filename, fileFormat)
else
  call convert_grd(pGrd, filename)
end if
end subroutine
end submodule

submodule(GlodenGrd) GlodenGrd_DSBB
implicit none
contains
!****************************************************************************
! 读取 DSBB 格式文件
! 输入:
! filename     ... 文件名
! 输出：
! this         ... 数据
!****************************************************************************
module subroutine read_grd_DSBB(this, filename)
implicit none
class(grd_DSBB), intent(out):: this
character(*), intent(in):: filename
integer unt
open(newunit=unt,file=trim(filename),status='old',access='stream')
read(unt) this%id
read(unt) this%nx, this%ny
read(unt) this%xlo, this%xhi
read(unt) this%ylo, this%yhi
read(unt) this%zlo, this%zhi
if(allocated(this%z)) deallocate(this%z)
allocate(this%z(this%nx,this%ny))
read(unt) this%z
close(unt)
this%zlo = minval(this%z)
this%zhi = maxval(this%z)
end subroutine

!****************************************************************************
! 将 DSBB 转为指定格式
! 输入:
! this         ... 数据
! filename     ... 输出文件名
! fileFormat   ... 输出文件的格式
!****************************************************************************
module subroutine convert_grd_DSBB(this, filename, fileFormat)
implicit none
class(grd_DSBB), target, intent(in):: this
character(*), intent(in):: filename
character(4), optional, intent(in):: fileFormat
class(grd), pointer:: pGrd
pGrd => this
if(present(fileFormat)) then
  call convert_grd(pGrd, filename, fileFormat)
else
  call convert_grd(pGrd, filename)
end if
end subroutine
end submodule

submodule(GlodenGrd) GlodenGrd_DSRB
implicit none
contains
!****************************************************************************
! 读取 DSRB 格式文件
! 输入:
! filename     ... 文件名
! 输出：
! this         ... 数据
!****************************************************************************
module subroutine read_grd_DSRB(this, filename)
implicit none
class(grd_DSRB), intent(out):: this
character(*), intent(in):: filename
integer unt
integer(8) p, i
open(newunit=unt,file=trim(filename),status='old',access='stream')
!读取 Header section
read(unt) this%tHeader%sectionName, this%tHeader%sectionSize
read(unt) this%version
!读取 Grid section
do
  read(unt) this%tGrid%sectionName, this%tGrid%sectionSize
  if(this%tGrid%sectionName=='GRID') exit
  inquire(unt,pos=p)
  read(unt,pos=p+this%tGrid%sectionSize)
end do
read(unt) this%ny, this%nx
read(unt) this%xlo, this%ylo
read(unt) this%xSize, this%ySize
read(unt) this%zlo, this%zhi
read(unt) this%rotation, this%blankValue
!读取 Data section
read(unt) this%tData%sectionName, this%tData%sectionSize
if(allocated(this%z)) deallocate(this%z)
allocate(this%z(this%nx,this%ny))
read(unt) this%z
this%zlo = minval(this%z)
this%zhi = maxval(this%z)
!读取 Fault Info Section, 不一定存在
if(allocated(this%traceIndex)) deallocate(this%traceIndex)
if(allocated(this%verticeLoc)) deallocate(this%verticeLoc)
read(unt,pos=1) !回到文件头
p = 1
do
  read(unt,pos=p,iostat=i) this%tFlti%sectionName, this%tFlti%sectionSize
  if(i/=0) then
    this%tFlti%sectionName = 'FLTI'
    this%tFlti%sectionSize = 0
    exit
  end if
  if(this%tFlti%sectionName=='FLTI') exit
  inquire(unt,pos=p)
  p = p + this%tFlti%sectionSize
end do
if(this%tFlti%sectionSize/=0) then
  read(unt) this%nTraces, this%nVertices
  allocate(this%traceIndex(2,this%nTraces)); allocate(this%verticeLoc(2,this%nVertices))
  read(unt) this%traceIndex
  read(unt) this%verticeLoc
end if
close(unt)
!白化值
select case(this%version)
case(1)
  where(this%z>=this%blankValue) this%z = noData
case(2)
  where(this%z==this%blankValue) this%z = noData
end select
end subroutine

!****************************************************************************
! 将 DSRB 转为指定格式
! 输入:
! this         ... 数据
! filename     ... 输出文件名
! fileFormat   ... 输出文件的格式
!****************************************************************************
module subroutine convert_grd_DSRB(this, filename, fileFormat)
implicit none
class(grd_DSRB), target, intent(in):: this
character(*), intent(in):: filename
character(4), optional, intent(in):: fileFormat
class(grd), pointer:: pGrd
pGrd => this
if(present(fileFormat)) then
  call convert_grd(pGrd, filename, fileFormat)
else
  call convert_grd(pGrd, filename)
end if
end subroutine
end submodule
program Test
use GlodenGrd
implicit none
class(grd), pointer:: a
call read_grd(a, '3.grd')  !读取grd文件
call a%convert('a2.grd','DSAA') !转存为 DSAA 格式文件

select type(a) !输出grd信息
type is(grd_DSRB)
  write(*,"(a)") a%id
  write(*,"(g0,2x,g0)") a%nx,  a%ny
  write(*,"(g0,2x,g0)") a%xlo, a%xlo + (a%nx-1)*a%xSize
  write(*,"(g0,2x,g0)") a%ylo, a%ylo + (a%ny-1)*a%ySize
  write(*,"(g0,2x,g0)") a%zlo, a%zhi
type is(grd_DSBB)
  write(*,"(a)") a%id
  write(*,"(g0,2x,g0)") a%nx,  a%ny
  write(*,"(g0,2x,g0)") a%xlo, a%xhi
  write(*,"(g0,2x,g0)") a%ylo, a%yhi
  write(*,"(g0,2x,g0)") a%zlo, a%zhi
type is(grd_DSAA)
  write(*,"(a)") a%id
  write(*,"(g0,2x,g0)") a%nx,  a%ny
  write(*,"(g0,2x,g0)") a%xlo, a%xhi
  write(*,"(g0,2x,g0)") a%ylo, a%yhi
  write(*,"(g0,2x,g0)") a%zlo, a%zhi
end select
deallocate(a) !注销内存
end program Test

 

...全文
1430 回复 打赏 收藏 转发到动态 举报
写回复
用AI写文章
回复
切换为时间正序
请发表友善的回复…
发表回复
surfer grd文件生成
用Fortran程序生成用于surfer绘图的grd文件
readsurferGRD_surfergrd文件_matlab、surfer_grdfile_读取grd_grd_
读取surfer 中的grd文件，可以将surfer中的GRD文件读取，或者将矩阵数据通过matlab写入到grd中
Surfer Grid 导入/导出：从 Surfer 和 Grap 等 Golden 软件应用程序导入和导出 GS ASCII（GRD 文件）-matlab开发
Surfer Grid，导入/导出。 Matlab <-> Golden Software Surfer 和 Grapher 该包包含两个简单的例程：grd_write.m和grd_read.m。 他们通过 GRD 文件格式（ASCII 版本）将 Matlab 与 Golden Software Surfer 通信。 grd_write（矩阵，xmin，xmax，ymin，ymax，namefile） 输入： 矩阵 = 要导出的矩阵xmin,xmax,ymin,ymax = 网格限制namefile = 要写入的文件的名称（包括“.grd”扩展名） 输出： 当前目录下的grd文件 [矩阵 xmin xmax ymin ymax]=grd_read（文件名） 输入： nomarch = 要读取的文件的名称，包括“.grd”扩展名输出： 矩阵 = 读取数据的矩阵xmin xmax ym
Surfer软件Grd数字化
Surfer软件目前比较流行，是制作平面图的重要软件。Grd数字化是图形编辑的需要，上传这个软件希望对大家有用
Surfer网格文件转Arcgis网格文件
Surfer网格文件转Arcgis网格文件 采用python语言实现 论文文档
Fortran

187

社区成员

27

社区内容

发帖
与我相关
我的任务
社区描述
社区旨在便于Fortran编程语言交流学习，同时讨论数值计算、并行计算相关问题，共同发展fortran生态。
社区管理员
  • 地球屋里老师
加入社区
  • 近7日
  • 近30日
  • 至今

加载中

查看更多榜单
社区公告

自愿，自律，自强

试试用AI创作助手写篇文章吧

+ 用AI写文章