Fortran Coder

查看: 10473|回复: 2
打印 上一主题 下一主题

[求助] 从DEM中提取数据时一个点距无法正常输出

[复制链接]

22

帖子

4

主题

0

精华

熟手

F 币
319 元
贡献
154 点
跳转到指定楼层
楼主
发表于 2015-4-13 15:02:20 | 只看该作者 |只看大图 回帖奖励 |倒序浏览 |阅读模式
需要从DEM提取相应栅格的属性值,代码如下
[Fortran] 纯文本查看 复制代码
PROGRAM MAIN
    IMPLICIT NONE
    INTEGER::I,J,K
    INTEGER::DEMCOLS,DEMROWS    !DEM数据列数、行数
    INTEGER::DEMNODATA          !DEM无效值
    INTEGER,ALLOCATABLE,DIMENSION(:,:)::PROVALUE
    REAL::CORNERLON,CORNERLAT,CELLSIZE
    REAL,ALLOCATABLE,DIMENSION(:,:)::DEMLON,DEMLAT
    CHARACTER(LEN=30)::DEMTEMP
    OPEN(UNIT=1201503,FILE="DATA.txt")
    READ(1201503,*)DEMTEMP,DEMCOLS
    READ(1201503,*)DEMTEMP,DEMROWS
    READ(1201503,*)DEMTEMP,CORNERLON
    READ(1201503,*)DEMTEMP,CORNERLAT
    READ(1201503,*)DEMTEMP,CELLSIZE
    READ(1201503,*)DEMTEMP,DEMNODATA
    ALLOCATE(PROVALUE(DEMROWS,DEMCOLS),DEMLON(DEMROWS,DEMCOLS),DEMLAT(DEMROWS,DEMCOLS))
    DO I=1,DEMROWS
        DEMLAT(I,:)=CORNERLAT+(DEMROWS-I)*CELLSIZE
    END DO
    DO I=1,DEMCOLS
        DEMLON(:,I)=CORNERLON+(I-1)*CELLSIZE
    END DO
    DO I=1,DEMROWS
        READ(1201503,*)(PROVALUE(I,J),J=1,DEMCOLS)
    END DO
    OPEN(UNIT=2201501,FILE="ResultsCAREPOINTS.txt")
    DO J=1,DEMROWS
        DO K=1,DEMCOLS
            IF((DEMLON(J,K)==101.95.AND.DEMLAT(J,K)==27.05).OR.(DEMLON(J,K)==101.45.AND.DEMLAT(J,K)==26.95).OR.(DEMLON(J,K)==101.75.AND.DEMLAT(J,K)==26.65))THEN
                WRITE(2201501,'(2F10.3,I10)')DEMLON(J,K),DEMLAT(J,K),PROVALUE(J,K)
            END IF
        END DO
    END DO
END PROGRAM

数据如下: DATA.zip (31.26 KB, 下载次数: 1)
我需要提取经纬度坐标为(101.95,27.05),(101.45,26.95),(101.75,26.65)三个点的属性值
但是,结果只能提取出其中的两个点,如图

请各位大侠,帮我看一下
谢谢






分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

954

帖子

0

主题

0

精华

大师

F 币
184 元
贡献
75 点

规矩勋章元老勋章新人勋章水王勋章热心勋章

QQ
沙发
发表于 2015-4-13 16:27:27 | 只看该作者
编程时,时刻记住,浮点数是有误差的。

所以不做浮点数相等判断。

而改用差的绝对值小于某个很小的数。
[Fortran] 纯文本查看 复制代码
Real , parameter :: eps = 1.0e-5
If ( (abs(demlon(j,k)-101.95)<eps.and.abs(demlat(j,k)-27.05)<eps) .or. &
   (abs(demlon(j,k)-101.45)<eps.and.abs(demlat(j,k)-26.95)<eps) .or. &
   (abs(demlon(j,k)-101.75)<eps.and.abs(demlat(j,k)-26.65)<eps) )then
    WRITE(2201501,'(2F10.3,I10)')DEMLON(J,K),DEMLAT(J,K),PROVALUE(J,K)
END IF

22

帖子

4

主题

0

精华

熟手

F 币
319 元
贡献
154 点
板凳
 楼主| 发表于 2015-4-13 16:36:46 | 只看该作者
脑子糊涂了,咋忘了这一条,还忙活了半天,判断其他地方呢

谢谢啦
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-12-24 08:17

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表