PROGRAM SHP_toFrom_xyDIG ! Utility program to transform shapefiles containing polylines ! (shapefile format is explained at https://en.wikipedia.org/wiki/Shapefile ! and also at https://www.esri.com/library/whitepapers/pdfs/shapefile.pdf ) ! to produce simple ASCII text files with (x, y) coordinates in meters, ! in Peter Bird's .DIG format, which is explained at: ! http://peterbird.name/guide/dig_format.htm . ! (Can also transform such files in the other direction: _xy.DIG --> .SHP/.DBF .) ! Typically, each (x, y)-type .DIG file will eventually be converted to a ! (longitude, latitude)-type .DIG file by my other utility program Projector, ! which supports a wide selection of common map-projections. ! At present SHP_toFrom_xyDIG primarily converts fault-trace & state-line ! polylines to be used by my other programs: ! Shells/FiniteMap and/or NeoKinema/NeoKineMap and/or Restore/RetroMap; ! up to 4 header records are copied for each polyline. ! The practical significance of this utility program is that it allows ! commercial digitizing programs like Golden Software's Didger5 ! which can Export / Save As Type: ESRI Shapefile (*.shp) ! to be used in digitizing fault-traces and state-lines that will later be used in ! my suite of .DIG-based programs. Also, conversion in the other direction ! allows existing (x, y)-type .DIG files to appear as a basemap layer in Didger5. ! By Peter Bird, UCLA, 2018.10.12 & 2019.03.08 IMPLICIT NONE CHARACTER*80 :: shapefile, token, dBase_file, dig_file CHARACTER*1, DIMENSION(:), ALLOCATABLE :: shp_bytes, dbf_bytes ! N.B. Indexing should start with 0 (0:...). CHARACTER*9 :: c_max_Ma, c_min_Ma CHARACTER*11, DIMENSION(:), ALLOCATABLE :: dbf_field_names CHARACTER*32 :: box1, box2, box3, box4, box5, boxN CHARACTER*50 :: period CHARACTER*64 :: line ! N.B. This is rather short, but that is because any polyline title needs to fit into 2 32-byte fields. CHARACTER*80 :: DIG_pathfile, text1, text2, text3, text4 CHARACTER*80, DIMENSION(4) :: text_store CHARACTER*300 :: age_string, name_string INTEGER :: content_length, & & dbf_bytes_per_record, dbf_fields, dbf_file_last_byte, dbf_file_length_KB, dbf_header_bytes, dbf_record_count, & & i, i_file, i1, i2, iConversionMode, iDotByte, ios, j, k, last_dbf_byte_index, last_shp_byte_index, last_dbf_byte_used, last_shp_byte_used, & & nPoints, nPolylines, numParts, numPoints, & & points_read, pre_record_LBU, & & record_length_bytes, record_number, shapetype, shp_file_last_byte, shp_file_length_KB, & & texts_read INTEGER, DIMENSION(:), ALLOCATABLE :: dbf_field_bytes, Parts LOGICAL :: another_input_file = .FALSE. ! (unless changed by user) LOGICAL :: using_dBase LOGICAL :: gap REAL*8 :: t1, t2, t3, t4, x, xMax, Xmax_m, xMin, Xmin_m, y, yMax, Ymax_m, yMin, Ymin_m REAL*8, DIMENSION(:, :), ALLOCATABLE :: Points ! (1, ... = X_m; 2, ... = Y_m) !Critical type-switching EQUIVALENCES: CHARACTER*1 :: group_C1 CHARACTER*2 :: group_C2 CHARACTER*4 :: group_C4 CHARACTER*8 :: group_C8 INTEGER*1 :: group_I1 INTEGER*2 :: group_I2 INTEGER*4 :: group_I4 REAL*8 :: group_R8 EQUIVALENCE ( group_C1, group_I1 ) EQUIVALENCE ( group_C2, group_I2 ) EQUIVALENCE ( group_C4, group_I4 ) EQUIVALENCE ( group_C8, group_R8 ) WRITE (*, "(' PROGRAM SHP_toFrom_xyDIG')") WRITE (*, "(' ')") WRITE (*, "(' Utility program to transform shapefiles containing polylines')") WRITE (*, "(' (shapefile format is explained at https://en.wikipedia.org/wiki/Shapefile')") WRITE (*, "(' and also at https://www.esri.com/library/whitepapers/pdfs/shapefile.pdf )')") WRITE (*, "(' to produce simple ASCII text files with (x, y) coordinates in meters,')") WRITE (*, "(' in Peter Bird''s .DIG format, which is explained at:')") WRITE (*, "(' http://peterbird.name/guide/dig_format.htm .')") WRITE (*, "(' (Can also transform such files in the other direction: _xy.DIG --> .SHP/.DBF .)')") WRITE (*, "(' Typically, each (x, y)-type .DIG file will eventually be converted to a ')") WRITE (*, "(' (longitude, latitude)-type .DIG file by my other utility program Projector,')") WRITE (*, "(' which supports a wide selection of common map-projections.')") WRITE (*, "(' At present SHP_toFrom_xyDIG primarily converts fault-trace & state-line')") WRITE (*, "(' polylines to be used by my other programs:')") WRITE (*, "(' Shells/FiniteMap and/or NeoKinema/NeoKineMap and/or Restore/RetroMap;')") WRITE (*, "(' up to 4 header records are copied for each polyline.')") WRITE (*, "(' The practical significance of this utility program is that it allows')") WRITE (*, "(' commercial digitizing programs like Golden Software''s Didger5')") WRITE (*, "(' which can Export / Save As Type: ESRI Shapefile (*.shp)')") WRITE (*, "(' to be used in digitizing fault-traces and state-lines that will later be used')") WRITE (*, "(' in my suite of .DIG-based programs. Also, conversion in the other direction')") WRITE (*, "(' allows existing (x, y)-type .DIG files to appear as a basemap layer in Didger5.')") WRITE (*, "(' By Peter Bird, UCLA, 2018.10.12 & 2019.03.08')") WRITE (*, "(' ')") CALL Pause() !-------------------------------------------------------------------------------------------------------- 10 WRITE (*, *) WRITE (*, "(' Which operation would you like to perform?')") WRITE (*, "(' ==========================================================================')") WRITE (*, "(' 1: Convert {token}.SHP and {token}.DBF ---> {token}_xy.DIG')") WRITE (*, "(' 2: Convert {token}_xy.DIG ---> {token}.SHP and {token}.DBF')") WRITE (*, "(' ==========================================================================')") iConversionMode = 1 CALL DPrompt_for_Integer ('Enter integer # from list above:', iConversionMode, iConversionMode) IF ((iConversionMode < 1).OR.(iConversionMode > 2)) THEN WRITE (*, "(' ERROR: Please select 1 or 2; no other legal choices.')") CALL Pause GO TO 10 END IF !============================================================================================================== IF (iConversionMode == 1) THEN ! Convert {token}.SHP and {token}.DBF ---> {token}_xy.DIG !============================================================================================================== shapefile = "YOUR.shp" 100 CONTINUE ! <=== TOP of user-controlled indefinite loop on multiple(?) input .SHP/.DBF file-pairs WRITE (*, *) CALL DPrompt_for_String(' Please name the existing shapefile you wish to convert:', shapefile, shapefile) iDotByte = INDEX (shapefile, '.', .TRUE.) IF (iDotByte == 0) THEN ! user probably left of the ".shp" filename-extension. shapefile = TRIM(shapefile) // ".shp" ! fix the error iDotByte = INDEX (shapefile, '.', .TRUE.) ! re-test END IF token = shapefile(1:(iDotByte-1)) dBase_file = TRIM(token) // ".dbf" !Try opening dBase_file (to see if it exists?): OPEN (UNIT = 2, FILE = dBase_file, STATUS = "OLD", IOSTAT = ios, FORM = "BINARY", RECL = 1) ! 1 byte at a time; ! requires compiler switch assume:byterecl, which can be set from menus ! under Project Properties / Fortran / Data / Use shp_bytes as RECL: YES. IF (ios /= 0) THEN WRITE (*, "(' NOTE: File ', A, ' not found (in current folder).')") TRIM(dBase_file) WRITE (*, "(' Proceeding without any dBase file; therefore')") WRITE (*, "(' any header records that identify your points/polylines will be LOST!')") CALL Pause() using_dBase = .FALSE. ELSE using_dBase = .TRUE. CLOSE (2) END IF !Set name for output .DIG file: dig_file = TRIM(token) // "_xy.dig" !-------------------------------------------------------------------------------------------------------- OPEN (UNIT = 1, FILE = shapefile, STATUS = "OLD", IOSTAT = ios, FORM = "BINARY", RECL = 1) ! 1 byte at a time; ! requires compiler switch assume:byterecl, which can be set from menus ! under Project Properties / Fortran / Data / Use shp_bytes as RECL: YES. IF (ios /= 0) THEN WRITE (*, "(' ERROR: File ', A, ' not found (in current folder).')") TRIM(shapefile) CALL Pause() STOP END IF WRITE (*, *) WRITE (*, "(' ----------------------------------------------------------------')") WRITE (*, *) WRITE (*, "(' ', A, ':')") TRIM(shapefile) IF (ALLOCATED(shp_bytes)) DEALLOCATE(shp_bytes) ! necessary when looping to process multiple shapefiles ALLOCATE ( shp_bytes(0:99) ) ! Just enough for reading the header, which contains file-length. READ (1, IOSTAT = ios) shp_bytes(0:99) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Attempt to READ shp_bytes 0:99 gave ios = ', I8)") CALL Pause() STOP END IF !Bytes 0:3 should contain Bigendian version of 0x0000270a: !Convert Bigendian to Littleendian: group_C4(1:1) = shp_bytes(3) group_C4(2:2) = shp_bytes(2) group_C4(3:3) = shp_bytes(1) group_C4(4:4) = shp_bytes(0) !EQUIVALENCE allows group_C4 to be accessed as INTEGER*4 :: group_I4: WRITE (*, "(' File code = ', Z10, ' (Compare to expected 0x0000270a)')") group_I4 !Bytes 24:47 should contain Bigendian version of file-length (integer, in shp_bytes): !Convert Bigendian to Littleendian: group_C4(1:1) = shp_bytes(27) group_C4(2:2) = shp_bytes(26) group_C4(3:3) = shp_bytes(25) group_C4(4:4) = shp_bytes(24) !EQUIVALENCE allows group_C4 to be accessed as INTEGER*4 :: group_I4: WRITE (*, "(' File length = ', I10, ' 16-bit words')") group_I4 group_I4 = group_I4 * 2 WRITE (*, "(' File length = ', I10, ' shp_bytes')") group_I4 shp_file_last_byte = group_I4 last_shp_byte_index = shp_file_last_byte - 1 ! because of (0:last_shp_byte_index) !Bytes 28:31 should contain Littleendian "Version" integer group_C4(1:1) = shp_bytes(28) group_C4(2:2) = shp_bytes(29) group_C4(3:3) = shp_bytes(30) group_C4(4:4) = shp_bytes(31) !EQUIVALENCE allows group_C4 to be accessed as INTEGER*4 :: group_I4: WRITE (*, "(' Version = ', I10)") group_I4 !Bytes 32:35 should contain Littleendian "shapetype" integer group_C4(1:1) = shp_bytes(32) group_C4(2:2) = shp_bytes(33) group_C4(3:3) = shp_bytes(34) group_C4(4:4) = shp_bytes(35) !EQUIVALENCE allows group_C4 to be accessed as INTEGER*4 :: group_I4: shapetype = group_I4 IF (shapetype == 3) THEN WRITE (*, "(' Shapetype = ', I10, ' (Polyline)')") shapetype ELSE IF (shapetype == 5) THEN WRITE (*, "(' Shapetype = ', I10, ' (Polygon)')") shapetype ELSE WRITE (*, "(' Shapetype = ', I10, ' (UNKNOWN)')") shapetype END IF last_shp_byte_used = 35 !Next 32 shp_bytes expected to contain Bounding Box as 4 consecutive Littleendian doubles: DO j = 1, 4 DO k = 1, 8 last_shp_byte_used = last_shp_byte_used + 1 group_C8(k:k) = shp_bytes(last_shp_byte_used) END DO ! k = 1, 8 !EQUIVALENCE allows group_C8 to be accessed as REAL*8 :: group_R8: IF (j == 1) WRITE (*, "(' Minimum Bounding Rectangle component = ', F15.3, ' (min X)')") group_R8 IF (j == 2) WRITE (*, "(' Minimum Bounding Rectangle component = ', F15.3, ' (min Y)')") group_R8 IF (j == 3) WRITE (*, "(' Minimum Bounding Rectangle component = ', F15.3, ' (max X)')") group_R8 IF (j == 4) WRITE (*, "(' Minimum Bounding Rectangle component = ', F15.3, ' (max Y)')") group_R8 END DO ! j = 1, 4 !skipping some pointless stuff about "Z" and "M", CLOSE the .SHP file: ! (Below, it will be re-OPENed, and read to the end, and memorized.) CLOSE(1) DEALLOCATE ( shp_bytes ) !-------------------------------------------------------------------------------------------------------- IF (using_dBase) THEN OPEN (UNIT = 2, FILE = dBase_file, STATUS = "OLD", IOSTAT = ios, FORM = "BINARY", RECL = 1) ! 1 byte at a time; ! requires compiler switch assume:byterecl, which can be set from menus ! under Project Properties / Fortran / Data / Use shp_bytes as RECL: YES. IF (ios /= 0) THEN WRITE (*, "(' ERROR: File ', A, ' not found (in current folder).')") TRIM(dBase_file) CALL Pause() STOP END IF WRITE (*, *) WRITE (*, "(' HEADER INFORMATION FROM dBase FILE ', A)") TRIM(dBase_file) IF (ALLOCATED(dbf_bytes)) DEALLOCATE(dbf_bytes) ! necessary when looping to process multiple shapefiles ALLOCATE ( dbf_bytes(0:354) ) ! Probably enough to read the header... READ (2, IOSTAT = ios) dbf_bytes(0:354) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Attempt to READ dbf_bytes 0:354 gave ios = ', I8)") ios CALL Pause() STOP END IF !Bytes 4-7 should contain # of records in this file. Assuming Littleendian(?): group_C4(1:1) = dbf_bytes(4) group_C4(2:2) = dbf_bytes(5) group_C4(3:3) = dbf_bytes(6) group_C4(4:4) = dbf_bytes(7) !EQUIVALENCE allows group_C4 to be accessed as INTEGER*4 :: group_I4: dbf_record_count = group_I4 WRITE (*, "(' dbf_record_count = ', I10)") dbf_record_count !Bytes 8-9 should tell how many bytes (INTEGER*2) are in the header. Assuming Littleendian(?): group_C2(1:1) = dbf_bytes(8) group_C2(2:2) = dbf_bytes(9) !EQUIVALENCE allows group_C2 to be accessed as INTEGER*2 :: group_I2: dbf_header_bytes = group_I2 WRITE (*, "(' dbf_header_bytes = ', I10)") dbf_header_bytes ! But, remember that the first of these is numbered "0" in dbf_bytes(). !Bytes 10-11 should tell how many bytes (INTEGER*2) per record. Assuming Littleendian(?): group_C2(1:1) = dbf_bytes(10) group_C2(2:2) = dbf_bytes(11) !EQUIVALENCE allows group_C2 to be accessed as INTEGER*2 :: group_I2: dbf_bytes_per_record = group_I2 WRITE (*, "(' dbf_bytes_per_record = ', I10)") dbf_bytes_per_record !Infer list of dbf fields (spreadsheet columns): dbf_fields = 0 last_dbf_byte_used = 31 ! As first field descriptor should begin in byte 32. fielding: DO IF (dbf_bytes(last_dbf_byte_used + 1) == CHAR(13)) EXIT fielding ! special end-of-list byte: 0x0D dbf_fields = dbf_fields + 1 last_dbf_byte_used = last_dbf_byte_used + 32 END DO fielding WRITE (*, "(' dbf_fields = ', I10)") dbf_fields ALLOCATE ( dbf_field_names(dbf_fields) ) ALLOCATE ( dbf_field_bytes(dbf_fields) ) dbf_fields = 0 last_dbf_byte_used = 31 ! As first field descriptor should begin in byte 32. refielding: DO IF (dbf_bytes(last_dbf_byte_used + 1) == CHAR(13)) EXIT refielding ! special end-of-list byte: 0x0D !First 11 bytes are ACII string for name of field dbf_fields = dbf_fields + 1 DO j = 1, 11 dbf_field_names(dbf_fields)(j:j) = dbf_bytes(last_dbf_byte_used + j) END DO !Byte "16" (from base of 0) contains the field length: group_C1 = dbf_bytes(last_dbf_byte_used + 17) !EQUIVALENCE allows this to be accessed as INTEGER*1 :: group_I1 dbf_field_bytes(dbf_fields) = group_I1 WRITE (*, "(' ', A11, ' (', I3,' bytes)')") dbf_field_names(dbf_fields), dbf_field_bytes(dbf_fields) last_dbf_byte_used = last_dbf_byte_used + 32 END DO refielding dbf_file_last_byte = dbf_header_bytes + (dbf_bytes_per_record * dbf_record_count) last_dbf_byte_index = dbf_file_last_byte - 1 CLOSE(2) ! to be re-OPENed below... DEALLOCATE ( dbf_bytes ) ! to be re-ALLOCATED, much longer, below ... END IF ! using_dBase; must decode the header and field structures. !-------------------------------------------------------------------------------------------------------- CALL Pause() ! Allowing user to read the header information from this/these file(s) ... !-------------------------------------------------------------------------------------------------------- !Re-READ the .SHP file, and memorize it: ALLOCATE ( shp_bytes(0:shp_file_last_byte) ) OPEN (UNIT = 1, FILE = shapefile, STATUS = "OLD", IOSTAT = ios, FORM = "BINARY", RECL = 1) ! 1 byte at a time; ! requires compiler switch assume:byterecl, which can be set from menus ! under Project Properties / Fortran / Data / Use shp_bytes as RECL: YES. READ (1, IOSTAT = ios) shp_bytes(0:last_shp_byte_index) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Attempt to read shp_bytes 0:', I10, ' gave ios = ', I8)") last_shp_byte_index, ios CALL Pause() STOP END IF CLOSE(1) !-------------------------------------------------------------------------------------------------------- IF (using_dBase) THEN !Re-READ the .DBF file, and memorize it: ALLOCATE ( dbf_bytes(0:dbf_file_last_byte) ) OPEN (UNIT = 2, FILE = dBase_file, STATUS = "OLD", IOSTAT = ios, FORM = "BINARY", RECL = 1) ! 1 byte at a time; ! requires compiler switch assume:byterecl, which can be set from menus ! under Project Properties / Fortran / Data / Use shp_bytes as RECL: YES. READ (2, IOSTAT = ios) dbf_bytes(0:last_dbf_byte_index) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Attempt to READ dbf_bytes 0:', I10, ' gave ios = ', I8)") last_dbf_byte_index, ios CALL Pause() STOP END IF CLOSE(2) END IF ! using_dBase !--------------------------------------------------------------------------------------------------------------- !OPEN the output "_xy.dig" file: OPEN (UNIT = 3, FILE = TRIM(dig_file)) ! unconditional OPEN will overwrite any existing file with this name !--------------------------------------------------------------------------------------------------------------- !Start indefinite loop on shapes: last_shp_byte_used = 99 ! at end of the .SHP header get_records: DO !Each record header has 2 parts; first is its "Record number (1-based)" in Bigendian i32: !Convert Bigendian to Littleendian: group_C4(1:1) = shp_bytes(last_shp_byte_used + 4) group_C4(2:2) = shp_bytes(last_shp_byte_used + 3) group_C4(3:4) = shp_bytes(last_shp_byte_used + 2) group_C4(4:4) = shp_bytes(last_shp_byte_used + 1) !EQUIVALENCE allows group_C4 to be accessed as INTEGER*4 :: group_I4: record_number = group_I4 WRITE (*, *) WRITE (*, "(' Record number (1-based) = ', I10)") record_number last_shp_byte_used = last_shp_byte_used + 4 !Second part of record header is its "Record length (in 16-bit words)" in Bigendian i32: !Convert Bigendian to Littleendian: group_C4(1:1) = shp_bytes(last_shp_byte_used + 4) group_C4(2:2) = shp_bytes(last_shp_byte_used + 3) group_C4(3:4) = shp_bytes(last_shp_byte_used + 2) group_C4(4:4) = shp_bytes(last_shp_byte_used + 1) !EQUIVALENCE allows group_C4 to be accessed as INTEGER*4 :: group_I4: WRITE (*, "(' Record length (in 16-bit words) = ', I10)") group_I4 record_length_bytes = 2 * group_I4 last_shp_byte_used = last_shp_byte_used + 4 pre_record_LBU = last_shp_byte_used !- - - - - - Beginning of one record (not including its header, above) - - - - - - - - - - - - - - !Next 4 shp_bytes should contain (redundant) Littleendian "shapetype" integer: group_C4(1:1) = shp_bytes(last_shp_byte_used + 1) group_C4(2:2) = shp_bytes(last_shp_byte_used + 2) group_C4(3:3) = shp_bytes(last_shp_byte_used + 3) group_C4(4:4) = shp_bytes(last_shp_byte_used + 4) !EQUIVALENCE allows group_C4 to be accessed as INTEGER*4 :: group_I4: shapetype = group_I4 last_shp_byte_used = last_shp_byte_used + 4 IF (shapetype == 3) THEN ! Polyline ----------------------------------------------------------- !Next 32 shp_bytes expected to contain Bounding Box as 4 consecutive Littleendian doubles: DO j = 1, 4 DO k = 1, 8 last_shp_byte_used = last_shp_byte_used + 1 group_C8(k:k) = shp_bytes(last_shp_byte_used) END DO ! k = 1, 8 !EQUIVALENCE allows group_C8 to be accessed as REAL*8 :: group_R8: IF (j == 1) Xmin_m = group_R8 IF (j == 2) Ymin_m = group_R8 IF (j == 3) Xmax_m = group_R8 IF (j == 4) Xmax_m = group_R8 END DO ! j = 1, 4 !Next 4 shp_bytes contain Littleendian "NumParts" integer group_C4(1:1) = shp_bytes(last_shp_byte_used + 1) group_C4(2:2) = shp_bytes(last_shp_byte_used + 2) group_C4(3:3) = shp_bytes(last_shp_byte_used + 3) group_C4(4:4) = shp_bytes(last_shp_byte_used + 4) !EQUIVALENCE allows group_C4 to be accessed as INTEGER*4 :: group_I4: NumParts = group_I4 last_shp_byte_used = last_shp_byte_used + 4 !Next 4 shp_bytes contain Littleendian "NumPoints" integer group_C4(1:1) = shp_bytes(last_shp_byte_used + 1) group_C4(2:2) = shp_bytes(last_shp_byte_used + 2) group_C4(3:3) = shp_bytes(last_shp_byte_used + 3) group_C4(4:4) = shp_bytes(last_shp_byte_used + 4) !EQUIVALENCE allows group_C4 to be accessed as INTEGER*4 :: group_I4: NumPoints = group_I4 last_shp_byte_used = last_shp_byte_used + 4 ALLOCATE ( Parts(NumParts) ) ! Offset-index of its first point, in the Points array. First value should be 0. ALLOCATE ( Points(2, NumPoints) ) ! REAL*8; (1, == X; 2, == Y) DO i = 1, NumParts !Next 4 shp_bytes contain 1 of the Littleendian "Parts" integers group_C4(1:1) = shp_bytes(last_shp_byte_used + 1) group_C4(2:2) = shp_bytes(last_shp_byte_used + 2) group_C4(3:3) = shp_bytes(last_shp_byte_used + 3) group_C4(4:4) = shp_bytes(last_shp_byte_used + 4) !EQUIVALENCE allows group_C4 to be accessed as INTEGER*4 :: group_I4: Parts(i) = group_I4 last_shp_byte_used = last_shp_byte_used + 4 END DO ! i = 1, NumParts DO i = 1, NumPoints !Next 16 shp_bytes expected to contain (X_m, Y_m) as 2 consecutive Littleendian doubles: DO j = 1, 2 DO k = 1, 8 last_shp_byte_used = last_shp_byte_used + 1 group_C8(k:k) = shp_bytes(last_shp_byte_used) END DO ! k = 1, 8 !EQUIVALENCE allows group_C8 to be accessed as REAL*8 :: group_R8: Points(j, i) = group_R8 END DO ! j = 1, 2 END DO ! i = 1, Points !--------------Get corresponding (desired) field(s) for this polyline from the .dbf file? ---- IF (using_dBase) THEN !The following block of code is specific to .dbf files produced by Golden Software's Didger5 (using default settings). !Other .dbf files with different fields will require small changes. IF (dbf_fields >= 1) THEN i1 = dbf_header_bytes + (dbf_bytes_per_record * (record_number - 1)) + 1 i2 = i1 + dbf_field_bytes(1) - 1 text1 = ' ' DO j = i1, i2 k = (j - i1) + 1 text1(k:k) = dbf_bytes(j) END DO ELSE text1 = ' ' END IF IF (dbf_fields >= 2) THEN i1 = dbf_header_bytes + (dbf_bytes_per_record * (record_number - 1)) + 1 + dbf_field_bytes(1) i2 = i1 + dbf_field_bytes(2) - 1 text2 = ' ' DO j = i1, i2 k = (j - i1) + 1 text2(k:k) = dbf_bytes(j) END DO ELSE text2 = ' ' END IF IF (dbf_fields >= 3) THEN i1 = i2 + 1 i2 = i1 + dbf_field_bytes(3) - 1 text3 = ' ' DO j = i1, i2 k = (j - i1) + 1 text3(k:k) = dbf_bytes(j) END DO ELSE text3 = ' ' END IF IF (dbf_fields >= 4) THEN i1 = i2 + 1 i2 = i1 + dbf_field_bytes(4) - 1 text4 = ' ' DO j = i1, i2 k = (j - i1) + 1 text4(k:k) = dbf_bytes(j) END DO ELSE text4 = ' ' END IF ELSE ! using_dBase = .FALSE. text1 = ' ' text2 = ' ' ! N.B. Any text-header that remains blank will be skipped, not written-out. text3 = ' ' text4 = ' ' END IF ! using_dBase, or not !----------------------------------------------------------------------- !---------------output section (for Polyline, assuming UNIT = 3 is already OPENed )--- IF (NumParts == 1) THEN IF (LEN_TRIM(text1) > 0) WRITE (3, "(A)") TRIM(text1) IF (LEN_TRIM(text2) > 0) WRITE (3, "(A)") TRIM(text2) IF (LEN_TRIM(text3) > 0) WRITE (3, "(A)") TRIM(text3) IF (LEN_TRIM(text4) > 0) WRITE (3, "(A)") TRIM(text4) DO i = 1, NumPoints WRITE (3, "(1X, SP, ES12.5, ',', ES12.5)") Points(1, i), Points(2, i) END DO ! i = 1, NumPoints WRITE (3, "('*** end of polyline ***')") ELSE ! multiple Parts DO j = 1, NumParts i1 = Parts(j) + 1 ! first value will be 0 + 1 = 1; if first fault has 10 points, then second value will be 10 + 1 = 11. IF (j < NumParts) THEN i2 = Parts(j + 1) ! -1 + 1; == Parts(j + 1) ELSE ! j == NumParts i2 = NumPoints END IF WRITE (3, "(A, ' part ', I5)") TRIM(text1), j IF (LEN_TRIM(text2) > 0) WRITE (3, "(A)") TRIM(text2) IF (LEN_TRIM(text3) > 0) WRITE (3, "(A)") TRIM(text3) IF (LEN_TRIM(text4) > 0) WRITE (3, "(A)") TRIM(text4) DO i = i1, i2 WRITE (3, "(1X, SP, ES12.5, ',', ES12.5)") Points(1, i), Points(2, i) END DO ! i = i1, i2 WRITE (3, "('*** end of polyline ***')") END DO ! j = 1, NumParts END IF ! NumParts == 1, or greater !----------------------------------------------------------------------- DEALLOCATE ( Points ) DEALLOCATE ( Parts ) ELSE IF (shapetype == 5) THEN ! Polygon-------------------------------------------------------- !NOTE that shapetype 5 (Polygon) has SAME file structure as shapetype 3 (Polyline) above, ! so the code in this block is copied from above and only slightly modified. ! The creator of the shapefile is supposed to have (already) ensured that every ! polygon closes with last (x, y) == first (x, y). ! I am choosing to reverse the shapefile convention about point order, by ! reversing outcrop perimeters from CW to CCW (as in all my other graphical and F-E codes). ! It remains to be seen how I should handle "holes" in multi-part polygons, ! but for now, these will be reversed as well, from CCW to CW. !Next 32 shp_bytes expected to contain Bounding Box as 4 consecutive Littleendian doubles: DO j = 1, 4 DO k = 1, 8 last_shp_byte_used = last_shp_byte_used + 1 group_C8(k:k) = shp_bytes(last_shp_byte_used) END DO ! k = 1, 8 !EQUIVALENCE allows group_C8 to be accessed as REAL*8 :: group_R8: IF (j == 1) Xmin_m = group_R8 IF (j == 2) Ymin_m = group_R8 IF (j == 3) Xmax_m = group_R8 IF (j == 4) Xmax_m = group_R8 END DO ! j = 1, 4 !Next 4 shp_bytes contain Littleendian "NumParts" integer group_C4(1:1) = shp_bytes(last_shp_byte_used + 1) group_C4(2:2) = shp_bytes(last_shp_byte_used + 2) group_C4(3:3) = shp_bytes(last_shp_byte_used + 3) group_C4(4:4) = shp_bytes(last_shp_byte_used + 4) !EQUIVALENCE allows group_C4 to be accessed as INTEGER*4 :: group_I4: NumParts = group_I4 last_shp_byte_used = last_shp_byte_used + 4 !Next 4 shp_bytes contain Littleendian "NumPoints" integer group_C4(1:1) = shp_bytes(last_shp_byte_used + 1) group_C4(2:2) = shp_bytes(last_shp_byte_used + 2) group_C4(3:3) = shp_bytes(last_shp_byte_used + 3) group_C4(4:4) = shp_bytes(last_shp_byte_used + 4) !EQUIVALENCE allows group_C4 to be accessed as INTEGER*4 :: group_I4: NumPoints = group_I4 last_shp_byte_used = last_shp_byte_used + 4 ALLOCATE ( Parts(NumParts) ) ! Offset-index of its first point, in the Points array. First value should be 0. ALLOCATE ( Points(2, NumPoints) ) ! REAL*8; (1, == X; 2, == Y) DO i = 1, NumParts !Next 4 shp_bytes contain 1 of the Littleendian "Parts" integers group_C4(1:1) = shp_bytes(last_shp_byte_used + 1) group_C4(2:2) = shp_bytes(last_shp_byte_used + 2) group_C4(3:3) = shp_bytes(last_shp_byte_used + 3) group_C4(4:4) = shp_bytes(last_shp_byte_used + 4) !EQUIVALENCE allows group_C4 to be accessed as INTEGER*4 :: group_I4: Parts(i) = group_I4 last_shp_byte_used = last_shp_byte_used + 4 END DO ! i = 1, NumParts DO i = 1, NumPoints !Next 16 shp_bytes expected to contain (X_m, Y_m) as 2 consecutive Littleendian doubles: DO j = 1, 2 DO k = 1, 8 last_shp_byte_used = last_shp_byte_used + 1 group_C8(k:k) = shp_bytes(last_shp_byte_used) END DO ! k = 1, 8 !EQUIVALENCE allows group_C8 to be accessed as REAL*8 :: group_R8: Points(j, i) = group_R8 END DO ! j = 1, 2 END DO ! i = 1, Points !--------------Get corresponding (desired) field(s) from .dbf file? ---- IF (using_dBase) THEN !The following block of code is specific to .dbf files produced by Golden Software's Didger5 (using default settings). !Other .dbf files with different fields will require small changes. IF (dbf_fields >= 1) THEN i1 = dbf_header_bytes + (dbf_bytes_per_record * (record_number - 1)) + 1 i2 = i1 + dbf_field_bytes(1) - 1 text1 = ' ' DO j = i1, i2 k = (j - i1) + 1 text1(k:k) = dbf_bytes(j) END DO ELSE text1 = ' ' END IF IF (dbf_fields >= 2) THEN i1 = dbf_header_bytes + (dbf_bytes_per_record * (record_number - 1)) + 1 + dbf_field_bytes(1) i2 = i1 + dbf_field_bytes(2) - 1 text2 = ' ' DO j = i1, i2 k = (j - i1) + 1 text2(k:k) = dbf_bytes(j) END DO ELSE text2 = ' ' END IF IF (dbf_fields >= 3) THEN i1 = i2 + 1 i2 = i1 + dbf_field_bytes(3) - 1 text3 = ' ' DO j = i1, i2 k = (j - i1) + 1 text3(k:k) = dbf_bytes(j) END DO ELSE text3 = ' ' END IF IF (dbf_fields >= 4) THEN i1 = i2 + 1 i2 = i1 + dbf_field_bytes(4) - 1 text4 = ' ' DO j = i1, i2 k = (j - i1) + 1 text4(k:k) = dbf_bytes(j) END DO ELSE text4 = ' ' END IF ELSE ! using_dBase = .FALSE. text1 = ' ' text2 = ' ' ! N.B. Any text-header that remains blank will be skipped, not written-out. text3 = ' ' text4 = ' ' END IF ! using_dBase, or not !----------------------------------------------------------------------- !---------------output section (for Polygon) --------------------------- IF (NumParts == 1) THEN IF (LEN_TRIM(text1) > 0) WRITE (3, "(A)") TRIM(text1) IF (LEN_TRIM(text2) > 0) WRITE (3, "(A)") TRIM(text2) IF (LEN_TRIM(text3) > 0) WRITE (3, "(A)") TRIM(text3) IF (LEN_TRIM(text4) > 0) WRITE (3, "(A)") TRIM(text4) DO i = 1, NumPoints WRITE (3, "(1X, SP, ES12.5, ',', ES12.5)") Points(1, i), Points(2, i) END DO ! i = 1, NumPoints WRITE (3, "('*** end of polygon ***')") ELSE ! multiple Parts DO j = 1, NumParts i1 = Parts(j) + 1 ! first value will be 0 + 1 = 1; if first fault has 10 points, then second value will be 10 + 1 = 11. IF (j < NumParts) THEN i2 = Parts(j + 1) ! -1 + 1; == Parts(j + 1) ELSE ! j == NumParts i2 = NumPoints END IF WRITE (3, "(A, ' part ', I5)") TRIM(text1), j IF (LEN_TRIM(text2) > 0) WRITE (3, "(A)") TRIM(text2) IF (LEN_TRIM(text3) > 0) WRITE (3, "(A)") TRIM(text3) IF (LEN_TRIM(text4) > 0) WRITE (3, "(A)") TRIM(text4) DO i = i1, i2 WRITE (3, "(1X, SP, ES12.5, ',', ES12.5)") Points(1, i), Points(2, i) END DO ! i = i1, i2 WRITE (3, "('*** end of polygon ***')") END DO ! j = 1, NumParts END IF ! NumParts == 1, or greater !----------------------------------------------------------------------- DEALLOCATE ( Points ) DEALLOCATE ( Parts ) ELSE ! shapetype is another integer (not 3 Polyline or 5 Polygon) that I have not coded yet! WRITE (*, *) WRITE (*, "(' ERROR: No Fortran currently in main-loop to handle shapetype = ', I6)") shapetype CALL Pause() STOP END IF ! shapetype is something I have coded, or something remaining to be done! !- - - - - - End of one record - - - - - - - - - - - - - - last_shp_byte_used = pre_record_LBU + record_length_bytes ! (should be redundant; just for safety) !WRITE (*, *) !CALL Pause() IF (last_shp_byte_used >= last_shp_byte_index) EXIT get_records END DO get_records IF (using_dBase) DEALLOCATE ( dbf_field_bytes ) IF (using_dBase) DEALLOCATE ( dbf_field_names ) WRITE (*, *) CALL Pause() CLOSE (3) ! output file !================================================================================================================ WRITE (*, *) CALL DPrompt_for_Logical('Do you wish to process another shapefile?', another_input_file, another_input_file) IF (another_input_file) GO TO 100 ! <=== BOTTOM of user-controlled indefinate loop. !============================================================================================================== ELSE IF (iConversionMode == 2) THEN ! Convert {token}_xy.DIG ---> {token}.SHP and {token}.DBF !============================================================================================================== dig_file = "YOUR_xy.dig" 2000 CONTINUE ! <=== TOP of user-controlled indefinite loop on multiple(?) input _xy.DIG files WRITE (*, *) CALL DPrompt_for_String(' Please name the existing .DIG file you wish to convert:', dig_file, dig_file) iDotByte = INDEX (dig_file, '.', .TRUE.) IF (iDotByte == 0) THEN ! user probably left off the ".dig" filename-extension. dig_file = TRIM(dig_file) // ".dig" ! fix the error iDotByte = INDEX (shapefile, '.', .TRUE.) ! re-test END IF token = dig_file(1:(iDotByte-1)) !But, if input file name already includes "_xy", then omit this from the token: IF ((INDEX(token, "_xy", .FALSE.) > 0).OR.(INDEX(token, "_XY", .FALSE.) > 0)) THEN token = dig_file(1:(iDotByte-4)) END IF OPEN (UNIT = 1, FILE = TRIM(dig_file), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Input file ', A, ' not found (in current folder).')") TRIM(dig_file) CALL Pause GO TO 2000 END IF !Just count vertices (digitized points) and polylines in .DIG input file, !and simultaneously determine its bounding-box coordinates (xMin ~ xMax, yMin ~ yMax), in meters. nPoints = 0 nPolylines = 0 xMin = +9.99D99 xMax = -9.99D99 yMin = +9.99D99 yMax = -9.99D99 scanning: DO READ (1, "(A)", IOSTAT = ios) line IF (ios /= 0) EXIT scanning ! probably hit the EOF IF (LEN_TRIM(line) == 0) THEN !Take no action for blank line. (These are sometimes accidentally included at the end.) ELSE IF (line(1:3) == "***") THEN nPolylines = nPolylines + 1 ELSE IF (line(1:1) == ' ') THEN ! probably a vertex, but let's be sure: READ (line, *, IOSTAT = ios) x, y IF (ios == 0) THEN nPoints = nPoints + 1 xMin = MIN(xMin, x) xMax = MAX(xMax, x) yMin = MIN(yMin, y) yMax = MAX(yMax, y) END IF ELSE ! line is probably some kind of title or data !No action required; each polyline will get exactly 4 title records in the .DBF END IF END DO scanning CLOSE (1) ! input _xy.DIG file !Report on the count that was just concluded: WRITE (*, *) WRITE (*, "(' Input file contains ', I10, ' polylines, including ', I10, ' total points.')") nPolylines, nPoints !ALLOCATE memory in which to build the shapefile: shp_file_last_byte = 99 + & ! <--- N.B. I have already subtracted 1 because indexing starts at 0; header is really 100 bytes (0:99). & ( 8 * nPolylines) + & ! 8 bytes per record-header (independent of object-class) & (48 * nPolylines) + & ! 48 bytes per polyline-header & (16 * nPoints) ! 8 bytes per REAL*8, so 16 bytes per vertex shp_file_length_KB = NINT((1 + shp_file_last_byte) / 1024.0D0) shp_file_length_KB = MAX(1, shp_file_length_KB) WRITE (*, "(' Allocating shapefile-sized memory block of ', I12, ' bytes = ', I8, ' KB.')") (1 + shp_file_last_byte), shp_file_length_KB !ALLOCATE memory in which to build the dBase file: dbf_file_last_byte = 193 + & ! header bytes, copied from an example output by Didger5 (includes generic field names). & nPolylines * (1 + 5*32) & ! leading byte for each record, plus 5 fields of 32 bytes each & +1 & ! final byte 0x1A to close the file & -1 ! to account for initial index of 0, not 1 dbf_file_length_KB = NINT((1 + dbf_file_last_byte) / 1024.0D0) dbf_file_length_KB = MAX(1, dbf_file_length_KB) WRITE (*, "(' Allocating dBasefile-sized memory block of ', I12, ' bytes = ', I8, ' KB.')") (1 + dbf_file_last_byte), dbf_file_length_KB CALL Pause ALLOCATE ( shp_bytes(0:shp_file_last_byte) ) ! Note that formulas above already allowed for starting index of 0. ALLOCATE ( dbf_bytes(0:dbf_file_last_byte) ) ! Note that formulas above already allowed for starting index of 0. !---------------------------------------------------------------------------------------------------------- !Actually BUILD the .SHP and .DBF files in memory: !- - - - - - Build the header for the shapefile: - - - - - - - - - - - - - - - - - - - - - - - - !Bytes 0:3 should contain Bigendian version of 0x0000270a: group_I4 = 10 + (7 * 16**2) + (2 * 16**3) !EQUIVALENCE allows INTEGER*4 :: group_I4 to be accessed as group_C4: !Convert Lettleendian to Bigendian: shp_bytes(3) = group_C4(1:1) shp_bytes(2) = group_C4(2:2) shp_bytes(1) = group_C4(3:3) shp_bytes(0) = group_C4(4:4) !Bytes 4:23 are unused; should contain 5 Bigendian copies of INTEGER*4 zero: group_I4 = 0 !EQUIVALENCE allows INTEGER*4 :: group_I4 to be accessed as group_C4: !Convert Lettleendian to Bigendian: shp_bytes(7) = group_C4(1:1) shp_bytes(6) = group_C4(2:2) shp_bytes(5) = group_C4(3:3) shp_bytes(4) = group_C4(4:4) shp_bytes(11) = group_C4(1:1) shp_bytes(10) = group_C4(2:2) shp_bytes( 9) = group_C4(3:3) shp_bytes( 8) = group_C4(4:4) shp_bytes(15) = group_C4(1:1) shp_bytes(14) = group_C4(2:2) shp_bytes(13) = group_C4(3:3) shp_bytes(12) = group_C4(4:4) shp_bytes(19) = group_C4(1:1) shp_bytes(18) = group_C4(2:2) shp_bytes(17) = group_C4(3:3) shp_bytes(16) = group_C4(4:4) shp_bytes(23) = group_C4(1:1) shp_bytes(22) = group_C4(2:2) shp_bytes(21) = group_C4(3:3) shp_bytes(20) = group_C4(4:4) !Bytes 24:27 should contain Bigendian version of file-length: group_I4 = shp_file_last_byte / 2 ! because NOW we are supposed to measure in 16-bit words, not 8-bit bytes. !EQUIVALENCE allows INTEGER*4 :: group_I4 to be accessed as group_C4: !Convert Littleendian to Bigendian: shp_bytes(27) = group_C4(1:1) shp_bytes(26) = group_C4(2:2) shp_bytes(25) = group_C4(3:3) shp_bytes(24) = group_C4(4:4) !Bytes 28:31 should contain Littleendian "Version" integer = 1000 group_I4 = 1000 !EQUIVALENCE allows INTEGER*4 :: group_I4 to be accessed as group_C4: shp_bytes(28) = group_C4(1:1) shp_bytes(29) = group_C4(2:2) shp_bytes(30) = group_C4(3:3) shp_bytes(31) = group_C4(4:4) !Bytes 32:35 should contain Littleendian "shapetype" integer shapetype = 3 ! Polyline group_I4 = shapetype !EQUIVALENCE allows INTEGER*4 :: group_I4 to be accessed as group_C4: shp_bytes(32) = group_C4(1:1) shp_bytes(33) = group_C4(2:2) shp_bytes(34) = group_C4(3:3) shp_bytes(35) = group_C4(4:4) last_shp_byte_used = 35 !Next 32 shp_bytes expected to contain Bounding Box as 4 consecutive Littleendian doubles: DO j = 1, 4 IF (j == 1) THEN group_R8 = xMin ELSE IF (j == 2) THEN group_R8 = yMin ELSE IF (j == 3) THEN group_R8 = xMax ELSE IF (j == 4) THEN group_R8 = yMax END IF !EQUIVALENCE allows REAL*8 :: group_R8 to be accessed as group_C8: DO k = 1, 8 last_shp_byte_used = last_shp_byte_used + 1 shp_bytes(last_shp_byte_used) = group_C8(k:k) END DO ! k = 1, 8 END DO ! j = 1, 4 (parameters of bounding-box) !Bytes 68:99 should have 4 copies of REAL*8 0.0D0: DO j = 1, 4 group_R8 = 0.0D0 ! zMin, zMax, mMin, mMax !EQUIVALENCE allows REAL*8 :: group_R8 to be accessed as group_C8: DO k = 1, 8 last_shp_byte_used = last_shp_byte_used + 1 shp_bytes(last_shp_byte_used) = group_C8(k:k) END DO ! k = 1, 8 END DO ! j = 1, 4 (Z and M dimensions of bounding-box; not in use) !- - - - - - Build the header for the dBase file: - - - - - - - - - - - - - - - - - - - - - - - - !First 4 bytes (0:3) copied from an example output by Didger5: dbf_bytes(0) = CHAR( 3) ! 03 in hexadecimal dbf_bytes(1) = CHAR(99) ! 63 in hexadecimal dbf_bytes(2) = CHAR( 1) ! 01 in hexadecimal dbf_bytes(3) = CHAR( 1) ! 01 in hexadecimal !Next 4 bytes (4:7) are the record count as a Littleendian(?) INTEGER*4 group_I4 = nPolylines !EQUIVALENCE allows INTEGER*4 group_I4 to be accessed as group_C4: dbf_bytes(4) = group_C4(1:1) dbf_bytes(5) = group_C4(2:2) dbf_bytes(6) = group_C4(3:3) dbf_bytes(7) = group_C4(4:4) !Next 2 bytes (8:9) are the the number of bytes in the .DBF header, as a Littleendian(?) INTEGER*2: group_I2 = 193 ! copied from an example output by Didger5 !EQUIVALENCE allows INTEGER*2 group_I2 to be accessed as group_C2: dbf_bytes(8) = group_C2(1:1) dbf_bytes(9) = group_C2(2:2) !Next 2 bytes (10:11) are the the number of bytes in each .DBF record, as a Littleendian(?) INTEGER*2: group_I2 = 161 ! copied from an example output by Didger5; == 1 + 5 * 32 !EQUIVALENCE allows INTEGER*2 group_I2 to be accessed as group_C2: dbf_bytes(10) = group_C2(1:1) dbf_bytes(11) = group_C2(2:2) !Bytes (12:31) are arcane; just copying from an example output by Didger5, where they are (0C:1F) in hex: DO j = 12, 28 dbf_bytes(j) = CHAR(0) ! 00 in hex END DO dbf_bytes(29) = CHAR(87) ! 57 in hex dbf_bytes(30) = CHAR( 0) ! 00 in hex dbf_bytes(31) = CHAR( 0) ! 00 in hex last_dbf_byte_used = 31 !Give names of 5 fields (32 bytes each): DO j = 1, 5 DO k = 1, 32 box1(k:k) = CHAR(0) ! all 32 bytes, as initialization END DO IF (j == 1) THEN box1(1:7) = "Primary" box1(12:12) = 'C' ELSE IF (j == 2) THEN box1(1:9) = "Secondary" box1(12:12) = 'C' ELSE IF (j == 3) THEN box1(1:8) = "Tertiary" box1(12:12) = 'C' ELSE IF (j == 4) THEN box1(1:10) = "Quaternary" box1(12:12) = 'C' ELSE IF (j == 5) THEN box1(1:6) = "ZLevel" box1(12:12) = 'N' END IF box1(17:17) = CHAR(32) ! hex 20; length of field as an INTEGER*1 DO k = 1, 32 last_dbf_byte_used = last_dbf_byte_used + 1 dbf_bytes(last_dbf_byte_used) = box1(k:k) END DO END DO !Terminate the .DBF header: last_dbf_byte_used = last_dbf_byte_used + 1 dbf_bytes(last_dbf_byte_used) = CHAR(13) ! hex 0x0D !--------------------------- Process the bulk of the .dig file: ------------------------------------------- OPEN (UNIT = 1, FILE = TRIM(dig_file), STATUS = "OLD", IOSTAT = ios) ! Read the .DIG file again, for all its data: ALLOCATE ( Points(2, nPoints) ) ! big enough, even if the .DIG file has only one polyline record_number = 0 ! Record numbers should begin at 1, but this will be bumped-up before first WRITE. polylining: DO text_store(1) = ' ' ! initialize as blank, before looking at file contents text_store(2) = ' ' text_store(3) = ' ' text_store(4) = ' ' texts_read = 0 points_read = 0 Points = 0.0D0 ! whole array (1:2, 1:nPoints) xMin = +9.99D99 ! to recompute a new bounding-box for this single polyline xMax = -9.99D99 yMin = xMin yMax = xMax !- - - - - - Read and memorize a single polyline from .DIG: - - - - - - - - - - - - - - - - - - - - - - - - internal: DO READ (1, "(A)", IOSTAT = ios) line IF (ios /= 0) EXIT polylining ! This is the exit taken at EOF, when no more polylines can be found in the .DIG file. IF (LEN_TRIM(line) == 0) THEN ! empty line !take no action ELSE IF (line(1:3) == "***") THEN EXIT internal ! This is the most common way to end the inner loop, at the end of each polyline. ELSE IF (line(1:1) == ' ') THEN ! most likely a vertex READ (line, *, IOSTAT = ios) x, y IF (ios == 0) THEN ! normal case points_read = points_read + 1 Points(1, points_read) = x Points(2, points_read) = y xMin = MIN(xMin, x) ! update the bounding-box for this single polyline xMax = MAX(xMax, x) yMin = MIN(yMin, y) yMax = MAX(yMax, y) ELSE ! treat as a (defective) header/parameter line texts_read = texts_read + 1 IF (texts_read <= 4) text_store(texts_read) = TRIM(line) END IF ELSE ! probably a header/parameter line: texts_read = texts_read + 1 IF (texts_read <= 4) text_store(texts_read) = TRIM(line) END IF ! different kinds of line in .DIG file END DO internal !- - - - - - Add this content to the .SHP and .DBF files: - - - - - - - - - - - - - - - - - - - - - - - - !Record header for the .SHP file (8 bytes): record_number = record_number + 1 ! So, 1 or higher, up to nPolylines. content_length = 24 + (8 * points_read) ! (In 16-bit (2-byte) words. Not counting the 8-byte record header we are creating now.) group_I4 = record_number !EQUIVALENCE allows INTEGER*4 :: group_I4 to be accessed as group_C4: !Convert Littleendian to Bigendian: shp_bytes(last_shp_byte_used + 1) = group_C4(4:4) shp_bytes(last_shp_byte_used + 2) = group_C4(3:3) shp_bytes(last_shp_byte_used + 3) = group_C4(2:2) shp_bytes(last_shp_byte_used + 4) = group_C4(1:1) last_shp_byte_used = last_shp_byte_used + 4 group_I4 = content_length !EQUIVALENCE allows INTEGER*4 :: group_I4 to be accessed as group_C4: !Convert Littleendian to Bigendian: shp_bytes(last_shp_byte_used + 1) = group_C4(4:4) shp_bytes(last_shp_byte_used + 2) = group_C4(3:3) shp_bytes(last_shp_byte_used + 3) = group_C4(2:2) shp_bytes(last_shp_byte_used + 4) = group_C4(1:1) last_shp_byte_used = last_shp_byte_used + 4 !Polyline record for the .SHP file: shapetype = 3 ! Polyline group_I4 = shapetype !EQUIVALENCE allows INTEGER*4 :: group_I4 to be accessed as group_C4: shp_bytes(last_shp_byte_used + 1) = group_C4(1:1) shp_bytes(last_shp_byte_used + 2) = group_C4(2:2) shp_bytes(last_shp_byte_used + 3) = group_C4(3:3) shp_bytes(last_shp_byte_used + 4) = group_C4(4:4) last_shp_byte_used = last_shp_byte_used + 4 group_R8 = xMin !EQUIVALENCE allows REAL*8 group_R8 to be accessed as group_C8: DO j = 1, 8 last_shp_byte_used = last_shp_byte_used + 1 shp_bytes(last_shp_byte_used) = group_C8(j:j) END DO group_R8 = yMin !EQUIVALENCE allows REAL*8 group_R8 to be accessed as group_C8: DO j = 1, 8 last_shp_byte_used = last_shp_byte_used + 1 shp_bytes(last_shp_byte_used) = group_C8(j:j) END DO group_R8 = xMax !EQUIVALENCE allows REAL*8 group_R8 to be accessed as group_C8: DO j = 1, 8 last_shp_byte_used = last_shp_byte_used + 1 shp_bytes(last_shp_byte_used) = group_C8(j:j) END DO group_R8 = yMax !EQUIVALENCE allows REAL*8 group_R8 to be accessed as group_C8: DO j = 1, 8 last_shp_byte_used = last_shp_byte_used + 1 shp_bytes(last_shp_byte_used) = group_C8(j:j) END DO numParts = 1 group_I4 = numParts !EQUIVALENCE allows INTEGER*4 :: group_I4 to be accessed as group_C4: shp_bytes(last_shp_byte_used + 1) = group_C4(1:1) shp_bytes(last_shp_byte_used + 2) = group_C4(2:2) shp_bytes(last_shp_byte_used + 3) = group_C4(3:3) shp_bytes(last_shp_byte_used + 4) = group_C4(4:4) last_shp_byte_used = last_shp_byte_used + 4 numPoints = points_read group_I4 = numPoints !EQUIVALENCE allows INTEGER*4 :: group_I4 to be accessed as group_C4: shp_bytes(last_shp_byte_used + 1) = group_C4(1:1) shp_bytes(last_shp_byte_used + 2) = group_C4(2:2) shp_bytes(last_shp_byte_used + 3) = group_C4(3:3) shp_bytes(last_shp_byte_used + 4) = group_C4(4:4) last_shp_byte_used = last_shp_byte_used + 4 group_I4 = 0 ! = relative offset index of first "part" within total list of points for this polyline !EQUIVALENCE allows INTEGER*4 :: group_I4 to be accessed as group_C4: shp_bytes(last_shp_byte_used + 1) = group_C4(1:1) shp_bytes(last_shp_byte_used + 2) = group_C4(2:2) shp_bytes(last_shp_byte_used + 3) = group_C4(3:3) shp_bytes(last_shp_byte_used + 4) = group_C4(4:4) last_shp_byte_used = last_shp_byte_used + 4 DO j = 1, points_read group_R8 = Points(1, j) !EQUIVALENCE allows REAL*8 group_R8 to be accessed as group_C8: DO k = 1, 8 last_shp_byte_used = last_shp_byte_used + 1 shp_bytes(last_shp_byte_used) = group_C8(k:k) END DO group_R8 = Points(2, j) !EQUIVALENCE allows REAL*8 group_R8 to be accessed as group_C8: DO k = 1, 8 last_shp_byte_used = last_shp_byte_used + 1 shp_bytes(last_shp_byte_used) = group_C8(k:k) END DO END DO !- - - - - - end of adding this polyline to the memory-image of the shapefile - - - - - - - !- - - - Append new data to growing memory-image of dBase-file: - - - - - - - - - - - - - - - box1 = ' ' box2 = ' ' box3 = ' ' box4 = ' ' box5 = ' ' ! contents of (unused) ZLevel field IF (LEN_TRIM(text_store(1)) > 32) THEN ! split it into 2 boxes: box1 = text_store(1)( 1:32) box2 = text_store(1)(33:64) ! (and, possibly truncate it?) IF (texts_read >= 2) box3 = text_store(2)( 1:32) IF (texts_read >= 3) box4 = text_store(3)( 1:32) ELSE ! short first-title line, such as "F3217N" box1 = text_store(1) IF (texts_read >= 2) box2 = text_store(2) IF (texts_read >= 3) box3 = text_store(3) IF (texts_read >= 4) box4 = text_store(4) END IF !Writing to the DBF memory-image always starts with a 0x20 byte to mark start-of-record: last_dbf_byte_used = last_dbf_byte_used + 1 dbf_bytes(last_dbf_byte_used) = CHAR(32) ! hex 20 !Now, write in the contents of the 4 32-byte boxes (and an empty 5th box for unused ZLevel): DO j = 1, 5 IF (j == 1) THEN boxN = box1 ELSE IF (j == 2) THEN boxN = box2 ELSE IF (j == 3) THEN boxN = box3 ELSE IF (j == 4) THEN boxN = box4 ELSE IF (j == 5) THEN boxN = box5 END IF DO k = 1, 32 last_dbf_byte_used = last_dbf_byte_used + 1 dbf_bytes(last_dbf_byte_used) = boxN(k:k) END DO END DO END DO polylining !Finish the .DBF memory-image with a final 0x1A byte: last_dbf_byte_used = last_dbf_byte_used + 1 dbf_bytes(last_dbf_byte_used) = CHAR(26) ! hex 1A CLOSE (1) ! input .DIG file DEALLOCATE ( Points ) !Output the shapefile that has been built in memory: shapefile = TRIM(token) // ".shp" OPEN (UNIT = 2, FILE = shapefile, IOSTAT = ios, FORM = "BINARY", RECL = 1) ! 1 byte at a time; ! requires compiler switch assume:byterecl, which can be set from menus ! under Project Properties / Fortran / Data / Use shp_bytes as RECL: YES. WRITE (2) shp_bytes(0:shp_file_last_byte) CLOSE (2) DEALLOCATE ( shp_bytes ) !Output the dBase file that has been built in memory: dBase_file = TRIM(token) // ".dbf" OPEN (UNIT = 3, FILE = dBase_file, IOSTAT = ios, FORM = "BINARY", RECL = 1) ! 1 byte at a time; ! requires compiler switch assume:byterecl, which can be set from menus ! under Project Properties / Fortran / Data / Use shp_bytes as RECL: YES. WRITE (3) dbf_bytes(0:dbf_file_last_byte) CLOSE (3) DEALLOCATE ( dbf_bytes ) WRITE (*, *) WRITE (*, "(' Shapefile ', A, ' and dBase-file ', A, ' have been written.')") TRIM(shapefile), TRIM(dBase_file) !================================================================================================================ WRITE (*, *) CALL DPrompt_for_Logical('Do you wish to process another {token}_xy.DIG file?', another_input_file, another_input_file) IF (another_input_file) GO TO 2000 ! <=== BOTTOM of user-controlled indefinate loop. !============================================================================================================== END IF ! iConversionMode == 1, or == 2 !============================================================================================================== WRITE (*, "(' ')") WRITE (*, "(' Job completed.')") CALL Pause() CONTAINS INTEGER FUNCTION DInt_Below (x) ! Returns integer equal to, or less than, REAL*8 x. ! (Note: INT() is different; always truncates toward zero.) IMPLICIT NONE REAL*8, INTENT(IN) :: x INTEGER :: i REAL*8 :: y i = DINT(x) IF (x >= 0.D0) THEN DInt_Below = i ELSE ! x < 0. y = 1.0D0 * i IF (y <= x) THEN DInt_Below = i ELSE ! most commonly DInt_Below = i - 1 END IF END IF END FUNCTION DInt_Below CHARACTER*10 FUNCTION DASCII10(x) ! Returns a right-justified 10-byte (or shorter) version of a ! REAL*8 number, in "human" format, with no more ! than 4 significant digits. IMPLICIT NONE REAL*8, INTENT(IN) :: x CHARACTER*10 :: temp10 CHARACTER*20 :: temp20 INTEGER :: j, k1, k10, zeros LOGICAL :: punt REAL*8 :: x_log DOUBLE PRECISION :: y IF (x == 0.0D0) THEN DASCII10=' 0' RETURN ELSE IF (x > 0.0D0) THEN punt = (x >= 999999999.5D0).OR.(x < 0.0000100D0) ELSE ! x < 0.0 punt = (x <= -99999999.5D0).OR.(x > -0.000100D0) END IF IF (punt) THEN ! need exponential notation; use Fortran utility WRITE (temp10,'(1P,E10.3)') x !consider possible improvements, from left to right: IF (temp10(3:6) == '.000') THEN ! right-shift 4 spaces over it temp20(7:10) = temp10(7:10) temp20(5:6) = temp10(1:2) temp20(1:4) = ' ' temp10 = temp20(1:10) ELSE IF (temp10(5:6) == '00') THEN ! right-shift 2 spaces over it temp20(7:10) = temp10(7:10) temp20(3:6) = temp10(1:4) temp20(1:2) = ' ' temp10 = temp20(1:10) ELSE IF (temp10(6:6) == '0') THEN ! right-shift 1 space over it temp20(7:10) = temp10(7:10) temp20(2:6) = temp10(1:5) temp20(1:1) = ' ' temp10 = temp20(1:10) END IF IF (temp10(8:8) == '+') THEN ! right-shift over + sign in exponent temp20(9:10) = temp10(9:10) temp20(2:8) = temp10(1:7) temp20(1:1) = ' ' temp10 = temp20(1:10) END IF IF (temp10(9:9) == '0') THEN ! right-shift over leading 0 in exponent temp20(10:10) = temp10(10:10) temp20(2:9) = temp10(1:8) temp20(1:1) = ' ' temp10 = temp20(1:10) END IF DASCII10 = temp10 ELSE ! can represent without exponential notation x_log = DLOG10(DABS(x)) zeros = DInt_Below(x_log) - 3 y = (10.D0**zeros) * NINT(DABS(x) / (10.D0**zeros)) IF (x < 0.0D0) y = -y WRITE (temp20,"(F20.9)") y ! byte 11 is the '.' !Avoid results like "0.7400001" due to rounding error! IF (temp20(19:20) == '01') temp20(19:20) = '00' !Find first important byte from right; change 0 -> ' ' k10 = 10 ! (if no non-0 found to right of .) right_to_left: DO j = 20, 12, -1 IF (temp20(j:j) == '0') THEN temp20(j:j) = ' ' ELSE k10 = j EXIT right_to_left END IF END DO right_to_left !put leading (-)0 before . of fractions, if it fits IF (x > 0.0D0) THEN IF (temp20(10:11) == ' .') temp20(10:11) = '0.' ELSE ! x < 0.0 IF (k10 <= 18) THEN IF (temp20(9:11) == ' -.') temp20(9:11) = '-0.' END IF END IF k1 = k10 - 9 DASCII10 = temp20(k1:k10) END IF ! punt, or not END FUNCTION DASCII10 SUBROUTINE DPrompt_for_Integer (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with an integer value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! This also occurs IF (mt_flashby), without waiting for user. ! Note that prompt_text should usually end with '?'. ! It can be more than 52 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text INTEGER, INTENT(IN) :: default INTEGER, INTENT(OUT) :: answer CHARACTER*11 :: instring, suggested INTEGER :: blank_at, bytes, ios, trial, written LOGICAL :: finished WRITE (suggested,"(I11)") default suggested = ADJUSTL(suggested) bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) written = 0 DO WHILE ((bytes - written) > 52) blank_at = written + INDEX(prompt_text((written+1):(written+52)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 52 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(suggested) finished = .TRUE. ! unless changed below READ (*,"(A)") instring IF (LEN_TRIM(instring) == 0) THEN answer = default ELSE !The following lead to occoasional abends !under Digital Visual Fortran 5.0D !(memory violations caught by WinNT): !READ (instring, *, IOSTAT = ios) trial !The following fix leads to a compiler error: !BACKSPACE (*) !READ (*, *, IOSTAT = ios) trial !and the following fix lead to an immediate abend: !BACKSPACE (5) !READ (*, *, IOSTAT = ios) trial !So, I am creating and then reading a dummy file: OPEN (UNIT = 72, FILE = 'trash') WRITE (72, "(A)") instring CLOSE (72) OPEN (UNIT = 72, FILE = 'trash') READ (72, *, IOSTAT = ios) trial CLOSE (72, STATUS = 'DELETE') IF (ios /= 0) THEN ! bad string WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") TRIM(instring) WRITE (*,"(' Enter an integer of 9 digits or less.')") WRITE (*,"(' Please try again:')") finished = .FALSE. ELSE answer = trial END IF ! problem with string, or not? END IF ! some bytes were entered END DO ! until finished END SUBROUTINE DPrompt_for_Integer SUBROUTINE DPrompt_for_Logical (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a logical value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! The same happens IF (mt_flashby), without waiting for the user. ! Note that prompt_text should usually end with '?'. ! It can be more than 70 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text LOGICAL, INTENT(IN) :: default LOGICAL, INTENT(OUT) :: answer CHARACTER*1 :: inbyte CHARACTER*3 :: yesno INTEGER :: blank_at, bytes, ios, written LOGICAL :: finished bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) IF (default) THEN yesno = 'Yes' ELSE yesno = 'No' END IF written = 0 DO WHILE ((bytes - written) > 70) blank_at = written + INDEX(prompt_text((written+1):(written+70)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 70 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(yesno) finished = .TRUE. ! unless changed below READ (*,"(A)") inbyte IF (LEN_TRIM(inbyte) == 0) THEN answer = default ELSE SELECT CASE (inbyte) CASE ('Y') answer = .TRUE. CASE ('y') answer = .TRUE. CASE ('T') answer = .TRUE. CASE ('t') answer = .TRUE. CASE ('R') answer = .TRUE. CASE ('r') answer = .TRUE. CASE ('O') answer = .TRUE. CASE ('o') answer = .TRUE. CASE ('N') answer = .FALSE. CASE ('n') answer = .FALSE. CASE ('F') answer = .FALSE. CASE ('f') answer = .FALSE. CASE ('W') answer = .FALSE. CASE ('w') answer = .FALSE. CASE DEFAULT WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") inbyte WRITE (*,"(' (Only the first letter of your answer is used.)')") WRITE (*,"(' To agree, enter Y, y, T, t, O, o, R, or r.')") WRITE (*,"(' To disagree, enter N, n, F, f, W, or w.')") WRITE (*,"(' Please try again:')") finished = .FALSE. END SELECT END IF ! a byte was entered END DO ! until finished END SUBROUTINE DPrompt_for_Logical SUBROUTINE DPrompt_for_Real (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a real value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! The same happens IF (mt_flashby), without waiting for the user. ! Note that prompt_text should usually end with '?'. ! It can be more than 52 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text REAL*8, INTENT(IN) :: default REAL*8, INTENT(OUT) :: answer CHARACTER*11 :: instring, suggested INTEGER :: blank_at, bytes, ios, point, written LOGICAL :: finished REAL*8 :: trial !------------------------------------------------------------------------------------ !This code worked (provided 4 significant digits), but left unecessary trailing zeros ("20.00"; "6.000E+07") !IF (((ABS(default) >= 0.1).AND.(ABS(default) < 1000.)).OR.(default == 0.0)) THEN ! ! Provide 4 significant digits by using Gxx.4 (the suffix shows significant digits, NOT digits after the decimal point!) ! WRITE (suggested,"(G11.4)") default !ELSE ! ! Use 1P,E because it avoids wasted and irritating leading 0 ("0.123E+4"). ! WRITE (suggested,"(1P,E11.3)") default !END IF !------------------------------------------------------------------------------------ !So I replaced it with the following: !(1) Use ASCII10 to get 4 significant digits (but no unecessary trailing zeroes): suggested = DASCII10(default) !(2) Be sure that the number contains some sign that it is floating-point, not integer: IF (INDEX(suggested, '.') == 0) THEN IF ((INDEX(suggested, 'E') == 0).AND.(INDEX(suggested, 'e') == 0).AND. & & (INDEX(suggested, 'D') == 0).AND.(INDEX(suggested, 'd') == 0)) THEN suggested = ADJUSTL(suggested) point = LEN_TRIM(suggested) + 1 suggested(point:point) = '.' END IF END IF !------------------------------------------------------------------------------------ suggested = ADJUSTL(suggested) bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) written = 0 DO WHILE ((bytes - written) > 52) blank_at = written + INDEX(prompt_text((written+1):(written+52)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 52 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(suggested) finished = .TRUE. ! unless changed below READ (*,"(A)") instring IF (LEN_TRIM(instring) == 0) THEN answer = default ELSE !The following lead to occoasional abends !under Digital Visual Fortran 5.0D !(memory violations caught by WinNT): !READ (instring, *, IOSTAT = ios) trial !The following fix leads to a compiler error: !BACKSPACE (*) !READ (*, *, IOSTAT = ios) trial !and the following fix lead to an immediate abend: !BACKSPACE (5) !READ (*, *, IOSTAT = ios) trial !So, I am creating and then reading a dummy file: OPEN (UNIT = 72, FILE = 'trash') WRITE (72, "(A)") instring CLOSE (72) OPEN (UNIT = 72, FILE = 'trash') READ (72, *, IOSTAT = ios) trial CLOSE (72, STATUS = 'DELETE') IF (ios /= 0) THEN ! bad string WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") TRIM(instring) WRITE (*,"(' Enter an real number using 11 characters (or less).')") WRITE (*,"(' Please try again:')") finished = .FALSE. ELSE answer = trial END IF ! problem with string, or not? END IF ! some bytes were entered END DO ! until finished END SUBROUTINE DPrompt_for_Real SUBROUTINE DPrompt_for_String (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a character-string value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! The same happens IF (mt_flashby), without waiting for the user. ! Note that prompt_text should usually end with '?'. ! It can be more than 70 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text CHARACTER*(*), INTENT(IN) :: default CHARACTER*(*), INTENT(OUT) :: answer CHARACTER*80 trial INTEGER :: blank_at, default_bytes, leftover, & & prompt_bytes, written prompt_bytes = LEN_TRIM(prompt_text) default_bytes = LEN_TRIM(default) written = 0 leftover = 79 - prompt_bytes - 4 ! unless changed below DO WHILE ((prompt_bytes - written) > 70) blank_at = written + INDEX(prompt_text((written+1):(written+70)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 70 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) leftover = 79 - (blank_at - (written+1) + 1) - 4 written = blank_at END DO IF (leftover >= default_bytes) THEN WRITE (*,"(' ',A,' [',A,']')") prompt_text(written+1:prompt_bytes), TRIM(default) ELSE WRITE (*,"(' ',A)") prompt_text(written+1:prompt_bytes) WRITE (*,"(' [',A,']')") TRIM(default) END IF WRITE (*,"(' ?: '\)") READ (*,"(A)") trial IF (LEN_TRIM(trial) == 0) THEN answer = TRIM(default) ELSE answer = TRIM(trial) END IF END SUBROUTINE DPrompt_for_String SUBROUTINE Pause() IMPLICIT NONE WRITE (*,"(' Press [Enter]...'\)") READ (*,*) END SUBROUTINE Pause END PROGRAM SHP_toFrom_xyDIG