从地形数据中提取一个矩形区域

梦纯

各位大虾，先谢谢了！
我现在想从一个280*200的中国及其周边地形格点数据topo.txt中提取出一个矩形区域，所以限定了i和j的范围，见以下程序：
program read_topo
implicit none

integer(kind=4),parameter::nx=280,ny=200
integer(kind=4) :: i,j
real topo(nx,ny),xgrid(ny),ygrid(nx)
!
!***open file
!
open(unit=12,file='topo.txt',form='formatted')
  do i=1,nx
     read(12,*) (topo(i,j),j=1,ny)
  end do
close(12)
!
!***select topo
!
  do i=1,nx
    do j=1,ny
        if(((i.ge.89).and.(i.le.198)) .and. ((j.ge.5).and.(j.le.98)))then
        topo(i,j)=-99
        end if
    end do
  end do

open(unit=44,file='after.txt',form='formatted')
  do i=1,nx
   write(44,'(200F10.4)') (topo(i,j),j=1,ny)
  end do
close(44)
!
!***write files
!

open(unit=13,file="topo_tp_new.dat",form='unformatted',access='direct',&
      status='replace',recl=4*nx*ny)
     write(13,rec=1)((topo(i,j),j=1,ny),i=1,nx)
     print*,topo(1,1)
close(13)
print*,topo(1,1),topo(200,200),topo(161,3)
end program read_topo
然后，我又写了ctl文件，如下：
dset topo_tp_new.dat
undef -99
options little_endian
xdef 280 linear 75.3785  0.25
ydef 200 linear 5.41628   0.25
zdef 1   levels 1
tdef 1   linear jan2008 1mon
vars 1
topo 0 99 topography
endvars
但画出来的图很诡异，见附图，不知道怎么回事，请大家帮忙看看，多谢了！

zhang710758774

你想要你取得矩形里的内容不显示？

mofangbao

http://bbs.06climate.com/forum.php?mod=viewthread&tid=18693
楼主可以参考下这个帖子，提问前多搜索~

梦纯

zhang710758774 发表于 2014-8-9 10:59
你想要你取得矩形里的内容不显示？

我就是想在这个topo.txt文件中限定一个矩形区域，为了直观，直接把这个区域内的值全部设成了缺省值。
就不知道这个是怎么回事，画出来的图这个么诡异。

梦纯

mofangbao 发表于 2014-8-9 18:43
http://bbs.06climate.com/forum.php?mod=viewthread&tid=18693
楼主可以参考下这个帖子，提问前多搜索~

谢了！不过我就是想请教一下我哪里出错了，嘿嘿。anyway，多谢哈！

mofangbao

梦纯 发表于 2014-8-13 07:48
谢了！不过我就是想请教一下我哪里出错了，嘿嘿。anyway，多谢哈！

那还得看看你的gs

zhang710758774

梦纯 发表于 2014-8-13 07:45
我就是想在这个topo.txt文件中限定一个矩形区域，为了直观，直接把这个区域内的值全部设成了缺省值。
就 ...

确实诡异。。
我照着你的样子做了一遍（可想而知这无聊的暑假），如图：


。。。。。。。
怎么没有图

梦纯

zhang710758774 发表于 2014-8-13 18:36
确实诡异。。
我照着你的样子做了一遍（可想而知这无聊的暑假），如图：

我试了好多遍，也出现过和你一样的情形，我觉得fortran程序里面的那个if语句可能有问题。但又找不出问题，真捉急！

梦纯

mofangbao 发表于 2014-8-13 08:55
那还得看看你的gs

我米有写gs文件，直接open *.ctl文件画了个图。清风大侠，是不是fortran程序里面的if语句有问题？多谢了

zhang710758774

梦纯 发表于 2014-8-14 06:24
我试了好多遍，也出现过和你一样的情形，我觉得fortran程序里面的那个if语句可能有问题。但又找不出问题 ...

你要是想要赋值的话，而且还是一块一块的，，直接topo(89:198,5:98)=-99
不用写if的

你就没想过是GrADS或者是ctl之类的有问题。。。