滑动平均法
PROGRAM HDPJ CHARACTER *(100) CMDFILE,INFILE,REGIONFILE,LOCALFILE,FILETYPE REAL EIGVAL INTEGER ITERATION WRITE(*,*)'INPUT CMDFILE-NAME: ' READ(*,*) CMDFILE CALL READCMD(CMDFILE,FILETYPE,INFILE,REGIONFILE,LOCALFILE,MLINE, NLINE,RADX,R ADY ,ITERA TION,EPS,EIGVAL) IF((FILETYPE.EQ.'GRD').OR.(FILETYPE.EQ.'GRD'))THEN CALL PROCESSGRD(INFILE,REGIONFILE,LOCALFILE,MLINE,NLINE,RADX, * * RADY ,ITERATION,EPS,EIGVAL) ELSE IF((FILETYPE.EQ.'BLN').OR.(FILETYPE.EQ.'BLN'))THEN
CALL PROCESSBLN(INFILE,REGIONFILE,LOCALFILE,MLINE,RADX, * ITERATION,EPS,EIGVAL)
! 参数输入子程序
SUBROUTINE READCMD(CMDFILE,FILETYPE,INFILE,REGIONFILE,LOCALFILE, MLINE,NLINE,RADX,RAD Y ,ITERA TION,EPS,EIGVAL) CHARACTER *(*)CMDFILE,FILETYPE,INFILE,REGIONFILE,LOCALFILE OPEN(10,FILE=CMDFILE,STATUS='OLD', ACCESS='SEQUENTIAL',FORM='FORMATTED') READ(10,*) FILETYPE READ(10,*) INFILE READ(10,*) REGIONFILE READ(10,*) LOCALFILE READ(10,*) RADX,RADY READ(10,*) ITERATION READ(10,*) EPS READ(10,*) EIGVAL CLOSE(10)
* ELSE IF((FILETYPE.EQ.'XYZ').OR.(FILETYPE.EQ.'XYZ'))THEN CALL PROCESSXYZ(INFILE,REGIONFILE,LOCALFILE,MLINE,RADX, ELSE WRITE(*,*)'无此类型' END IF END PROGRAM * RADY ,ITERATION,EPS,EIGVAL) *
IF((FILETYPE.EQ.'GRD').OR.(FILETYPE.EQ.'GRD'))THEN CALL MNGRD(INFILE,MLINE,NLINE) ELSE IF((FILETYPE.EQ.'BLN').OR.(FILETYPE.EQ.'BLN'))THEN ELSE IF((FILETYPE.EQ.'XYZ').OR.(FILETYPE.EQ.'XYZ'))THEN CALL MNXYZ(INFILE,MLINE) ELSE END IF END SUBROUTINE CALL MNBLN(INFILE,MLINE)
!*******************************************GRD************************************************** ! 个数输入子程序
!GRD 类型数据处理子程序
SUBROUTINE PROCESSGRD(INFILE,REGIONFILE,LOCALFILE,MLINE,NLINE, CHARACTER *(*)INFILE,REGIONFILE,LOCALFILE REAL,ALLOCATABLE::FIRST(:,:),REGION(:,:),LOCAL(:,:), ALLOCATE(FIRST(MLINE,NLINE),REGION(MLINE,NLINE), CALL READGRD(INFILE,MLINE,NLINE,XMIN,XMAX,YMIN, YMAX,FIRST) DX=(XMAX-XMIN)/(MLINE-1) DY=(YMAX-YMIN)/(NLINE-1) MR=INT(RADX/DX+0.5) NR=INT(RADY/DY+0.5) CALL ITERATION_GRD(MLINE,NLINE,MR,NR,ITERATION,EPS,EIGVAL, * RADX,RADY ,ITERA TION,EPS,EIGVAL) SUBROUTINE MNGRD(INFILE,MLINE,NLINE) CHARACTER *(*)INFILE OPEN(10,FILE=INFILE,STATUS='OLD',ACCESS='SEQUENTIAL', FORM='FORMATTED') READ(10,*) READ(10,*) MLINE,NLINE CLOSE(10) END SUBROUTINE * * REFIRST(:,:) * LOCAL(MLINE,NLINE),REFIRST(MLINE,NLINE)) * * FIRST,REFIRST,REGION,LOCAL)
CALL WRITEGRD(REGIONFILE,MLINE,NLINE,XMIN,XMAX,
* YMIN,YMAX,EIGVAL,REGION)
CALL WRITEGRD(LOCALFILE,MLINE,NLINE,XMIN,XMAX,
* YMIN,YMAX,EIGVAL,LOCAL)
DEALLOCA TE(FIRST,REGION,LOCAL,REFIRST)
END SUBROUTINE
! 读入GRD 文件
SUBROUTINE READGRD(INFILE,MLINE,NLINE,XMIN,XMAX,YMIN,
* YMAX,FIRST)
CHARACTER *(*)INFILE
REAL FIRST(MLINE,NLINE)
OPEN(10,FILE=INFILE,STATUS='OLD',ACCESS='SEQUENTIAL',
* FORM='FORMATTED')
READ(10,*)
READ(10,*)
READ(10,*) XMIN,XMAX
READ(10,*) YMIN,YMAX
READ(10,*)
READ(10,*) (((FIRST(I,J)),I=1,MLINE),J=1,NLINE)
CLOSE(10)
END SUBROUTINE
! 输出GRD 文件
SUBROUTINE WRITEGRD(OUTFILE,MLINE,NLINE,XMIN,XMAX,
* YMIN,YMAX,EIGVAL,OUTGRD)
CHARACTER *(*)OUTFILE
REAL OUTGRD(MLINE,NLINE)
ZMIN=HUGE(ZMIN)
ZMAX=-HUGE(ZMAX)
DO J=1,NLINE,1
DO I=1,MLINE,1
IF(OUTGRD(I,J).LT.0.5*EIGVAL)THEN
ZMIN=MIN(ZMIN,OUTGRD(I,J))
ZMAX=MAX(ZMAX,OUTGRD(I,J))
END IF
END DO
END DO OPEN(11,FILE=OUTFILE,STATUS='UNKNOWN', ACCESS='SEQUENTIAL',FORM='FORMATTED') WRITE(11,'(A)')'DSAA' WRITE(11,*) MLINE,NLINE WRITE(11,*) XMIN,XMAX WRITE(11,*) YMIN,YMAX WRITE(11,*) ZMIN,ZMAX DO J=1,NLINE WRITE(11,*) (OUTGRD(I,J),I=1,MLINE) END DO CLOSE(11) END SUBROUTINE *
! 滑动平均法（GRD ）
SUBROUTINE ITERATION_GRD(MLINE,NLINE,MR,NR,ITERATION,EPS,EIGVAL, REAL FIRST(MLINE,NLINE),REGION(MLINE,NLINE), EPSG=0.;ITE=1;REFIRST=FIRST DO WHILE(EPSG>EPS.OR.ITE
EPSG=0. ITE=ITE+1 DO I=1,MLINE,1 MRMIN=MAX(I-MR,1) MRMAX=MIN(I+MR,MLINE) DO J=1,NLINE,1 IF(REFIRST(I,J)
ENDIF ENDDO ENDDO REGION(I,J)=SUM/NUM EPSG=MAX(EPSG,ABS(REFIRST(I,J)-REGION(I,J))) REGION(I,J)=REFIRST(I,J) ELSE ENDIF ENDDO ENDDO REFIRST=REGION ENDDO DO J=1,NLINE,1 DO I=1,MLINE,1 IF(FIRST(I,J)
!*******************************************BLN************************************************** !BLN 个数输入程序
SUBROUTINE MNBLN(INFILE,MLINE) CHARACTER *(*)INFILE OPEN(10,FILE=INFILE,STATUS='OLD',ACCESS='SEQUENTIAL', FORM='FORMATTED') READ(10,*) MLINE CLOSE(10) END SUBROUTINE !BLN 文件处理子程序 SUBROUTINE PROCESSBLN(INFILE,REGIONFILE,LOCALFILE,MLINE,RADX, CHARACTER *(*)INFILE,REGIONFILE,LOCALFILE
* * ITERATION,EPS,EIGVAL)
REAL,ALLOCATABLE::FIRST(:),REGION(:),LOCAL(:),
* REFIRST(:),X(:)
ALLOCATE(FIRST(MLINE),REGION(MLINE),
* LOCAL(MLINE),REFIRST(MLINE),X(MLINE))
CALL READBLN(INFILE,MLINE,X,FIRST)
CALL ORDER_BLN(X,FIRST,MLINE)
CALL ITERATION_BLN(MLINE,RADX,ITERATION,EPS,EIGVAL,X,
* FIRST,REFIRST,REGION,LOCAL)
CALL WRITEBLN(INFILE,MLINE,X,FIRST)
CALL WRITEBLN(REGIONFILE,MLINE,X,REGION)
CALL WRITEBLN(LOCALFILE,MLINE,X,LOCAL)
DEALLOCA TE(X,FIRST,REGION,LOCAL,REFIRST)
END SUBROUTINE
! 读取BLN 文件
SUBROUTINE READBLN(INFILE,MLINE,X,FIRST)
CHARACTER *(*)INFILE
REAL X(MLINE),FIRST(MLINE)
OPEN(10,FILE=INFILE,STATUS='OLD',ACCESS='SEQUENTIAL',
* FORM='FORMATTED')
READ(10,*)
READ(10,*) ((X(I),FIRST(I)),I=1,MLINE)
CLOSE(10)
END SUBROUTINE
! 对BLN 文件排序
SUBROUTINE ORDER_BLN(X,FIRST,MLINE)
REAL X(MLINE),FIRST(MLINE)
INTEGER FLAG
FLAG=1
DO WHILE(FLAG==1)
FLAG=0
DO I=1,MLINE-1,1
IF(X(I+1)
FLAG=1
TEMP=X(I)
X(I)=X(I+1)
X(I+1)=TEMP
TEMP=FIRST(I)
FIRST(I)=FIRST(I+1)
FIRST(I+1)=TEMP
ENDIF
ENDDO
ENDDO
END SUBROUTINE
! 滑动平均法（BLN ）
SUBROUTINE ITERATION_BLN(MLINE,RADX,ITERATION,EPS,EIGVAL,X, * FIRST,REFIRST,REGION,LOCAL)
REAL X(MLINE),FIRST(MLINE),REGION(MLINE),
* LOCAL(MLINE),REFIRST(MLINE)
EPSG=0.;ITE=1;REFIRST=FIRST
DO WHILE(EPSG>EPS.OR.ITE
EPSG=0.
ITE=ITE+1
DO I=1,MLINE,1
IF(REFIRST(I)
J=I;SUM=0;NUM=0
DO WHILE( (X(J)+RADX)>=X(I) .AND. J>=1)
IF(REFIRST(J)
SUM=SUM+REFIRST(J)
NUM=NUM+1
ENDIF
J=J-1
ENDDO
J=I+1
DO WHILE( (X(I)+RADX)>=X(J) .AND. J
IF(REFIRST(J)
SUM=SUM+REFIRST(J)
NUM=NUM+1
ENDIF
J=J+1
ENDDO
REGION(I)=SUM/NUM
EPSG=MAX(EPSG,ABS(REFIRST(I)-REGION(I))) REGION(I)=REFIRST(I) ELSE ENDIF ENDDO REFIRST=REGION ENDDO DO I=1,MLINE,1 IF(FIRST(I)
! 输出BLN 文件
!*******************************************XYZ************************************************** !XYZ 个数检测程序
SUBROUTINE MNXYZ(FILENAME,NUMBER) CHARACTER*(*)FILENAME LOGICAL UNIT_OPEN,FILE_EXIST INTEGER NUNIT_IN NUNIT_IN=11 UNIT_OPEN=.TRUE.
SUBROUTINE WRITEBLN(OUTFILE,MLINE,X,OUTBLN) CHARACTER *(*)OUTFILE REAL OUTBLN(MLINE),X(MLINE) OPEN(11,FILE=OUTFILE,STATUS='UNKNOWN', ACCESS='SEQUENTIAL',FORM='FORMATTED') WRITE(11,*) MLINE DO I=1,MLINE,1 WRITE(11,*) X(I),OUTBLN(I) ENDDO CLOSE(11) END SUBROUTINE *
DO WHILE (UNIT_OPEN)
NUNIT_IN=NUNIT_IN+1
INQUIRE(UNIT=NUNIT_IN,OPENED=UNIT_OPEN)
END DO
INQUIRE(FILE=FILENAME,EXIST=FILE_EXIST)
IF (FILE_EXIST) THEN
OPEN(UNIT= NUNIT_IN,FILE=FILENAME,STATUS='OLD') NUMBER=0
DO WHILE (.NOT.EOF(NUNIT_IN))
READ(NUNIT_IN,*,END=200,ERR=200)X,F
NUMBER=NUMBER+1
200 END DO
CLOSE(UNIT=NUNIT_IN)
ELSE
WRITE(*,'(A,A,A)') '文件:',TRIM(FILENAME),' 不存在!'
END IF
END SUBROUTINE
!XYZ 文件处理子程序
SUBROUTINE PROCESSXYZ(INFILE,REGIONFILE,LOCALFILE,MLINE, * RADX,RADY ,ITERA TION,EPS,EIGVAL) CHARACTER *(*)INFILE,REGIONFILE,LOCALFILE
REAL,ALLOCATABLE::FIRST(:),REGION(:),LOCAL(:),
* REFIRST(:),XY(:,:)
ALLOCATE(FIRST(MLINE),REGION(MLINE),
* LOCAL(MLINE),REFIRST(MLINE),XY(2,MLINE))
CALL READXYZ(INFILE,MLINE,XY,FIRST)
CALL ORDER_XYZ(XY,FIRST,MLINE)
CALL ITERATION_XYZ(MLINE,RADX,RADY,ITERATION,EPS,EIGVAL,XY , * FIRST,REFIRST,REGION,LOCAL)
CALL WRITEXYZ(REGIONFILE,MLINE,XY,REGION)
CALL WRITEXYZ(LOCALFILE,MLINE,XY,LOCAL)
DEALLOCA TE(XY,FIRST,REGION,LOCAL,REFIRST)
END SUBROUTINE
! 读取XYZ 文件
SUBROUTINE READXYZ(INFILE,MLINE,XY,FIRST)
CHARACTER *(*)INFILE
REAL FIRST(MLINE),XY(2,MLINE)
OPEN(10,FILE=INFILE,STATUS='OLD',ACCESS='SEQUENTIAL', * FORM='FORMATTED')
DO I=1,MLINE,1
READ(10,*) XY(1,I),XY(2,I),FIRST(I)
ENDDO
CLOSE(10)
ENDSUBROUTINE
! 对XYZ 文件排序
SUBROUTINE ORDER_XYZ(XY,FIRST,MLINE)
REAL XY(2,MLINE),FIRST(MLINE)
INTEGER FLAG
FLAG=1
DO WHILE(FLAG==1)
FLAG=0
DO I=1,MLINE-1,1
IF(XY(1,I+1)
FLAG=1
TEMP=XY(1,I)
XY(1,I)=XY(1,I+1)
XY(1,I+1)=TEMP
TEMP=XY(2,I)
XY(2,I)=XY(2,I+1)
XY(2,I+1)=TEMP
TEMP=FIRST(I)
FIRST(I)=FIRST(I+1)
FIRST(I+1)=TEMP
ELSEIF(XY(1,I+1)==XY(1,I) .AND. XY(2,I+1)
TEMP=XY(2,I)
XY(2,I)=XY(2,I+1)
XY(2,I+1)=TEMP
TEMP=FIRST(I)
FIRST(I)=FIRST(I+1)
FIRST(I+1)=TEMP
ENDIF ENDDO END SUBROUTINE ENDDO
! 滑动平均法（XYZ ）
SUBROUTINE ITERATION_XYZ(MLINE,RADX,RADY,ITERATION,EPS,EIGVAL,XY , * FIRST,REFIRST,REGION,LOCAL)
REAL XY(2,MLINE),FIRST(MLINE),REGION(MLINE),
* LOCAL(MLINE),REFIRST(MLINE)
EPSG=0.;ITE=1;REFIRST=FIRST
DO WHILE(EPSG>EPS.OR.ITE
EPSG=0.
ITE=ITE+1
DO I=1,MLINE,1
IF(REFIRST(I)
J=I;SUM=0;NUM=0
DO WHILE( (XY(1,J)+RADX)>=XY(1,I) .AND. J>=1)
IF(ABS(XY(2,J)-XY(2,I))
* REFIRST(J)
SUM=SUM+REFIRST(J)
NUM=NUM+1
ENDIF
J=J-1
ENDDO
J=I+1
DO WHILE( (XY(1,I)+RADX)>=XY(1,J) .AND. J
* REFIRST(J)
SUM=SUM+REFIRST(J)
NUM=NUM+1
ENDIF
J=J+1
ENDDO
REGION(I)=SUM/NUM
EPSG=MAX(EPSG,ABS(REFIRST(I)-REGION(I)))
ELSE
11
REGION(I)=REFIRST(I)
ENDIF
ENDDO
REFIRST=REGION
ENDDO
DO I=1,MLINE,1
IF(FIRST(I)
LOCAL(I)=FIRST(I)-REGION(I)
ELSE
LOCAL(I)=FIRST(I)
ENDIF
ENDDO
ENDSUBROUTINE
! 输出XYZ 文件
SUBROUTINE WRITEXYZ(OUTFILE,MLINE,XY,OUTZ) CHARACTER *(*)OUTFILE
REAL OUTZ(MLINE),XY(2,MLINE)
OPEN(11,FILE=OUTFILE,STATUS='UNKNOWN',
* ACCESS='SEQUENTIAL',FORM='FORMATTED') DO I=1,MLINE,1
WRITE(11,*) XY(1,I),XY(2,I),OUTZ(I)
ENDDO
CLOSE(11)
END SUBROUTINE
12