Historical Shorelines for Puerto Rico from 1901 to 1987

Original shoreline shapefiles and source data (including ungeoreferenced topographic sheets [T- or TP-sheets], aerial photos, and USGS quadrangles) from OFR 93-574, Historical Shoreline Changes in Puerto Rico, 1901-1987 (Thieler and Danforth, 1993), were acquired from the archives of internal USGS data. The shapefiles received were separated by year based on the corresponding source information.
The first step of this QA/QC process was to make preliminary edits and additions to the original shoreline shapefiles. One of the first noticeable errors was the prevalence of duplicate features across the shorelines. For instance, identical features existed within the same shoreline shapefile, but each feature had its own unique attributes within the attribute table. To identify and delete these duplicate features, a "Length" field was added to each of the shapefiles. Using the Calculate Geometry tool within each shapefile's attribute table, the length of each shoreline feature was derived. With this new length field, a visual inspection the shapefiles' attribute tables with these identical shoreline lengths helped identify the duplicate shorelines.
There were also instances where there were identical features across multiple shoreline feature classes. In this case, a visual inspection was made in ArcMap and the Select by Location tool was used to identify these duplicate shorelines. If the target layer (meaning the feature class that contained the incorrectly duplicated shoreline) was spatially identical to the source layer (the feature class with the correct shoreline), then those selected shorelines were reviewed and deleted.
Each of these edited shapefiles was then exported as a feature class to a personal geodatabase for two reasons: for organizational purposes; and DSAS requires shorelines to be stored in a personal geodatabase (Himmelstoss and others, 2018). Fields that are necessary for the shoreline change analysis through DSAS as well as fields pertaining to other relevant information were added to each of the feature classes. These are described in the Attribute section of the metadata.
The historical shorelines were then assessed based on their accuracy by georeferencing the source data from which they originated. Two different shoreline sources are referenced: T-sheets and aerial photos.

The T-sheet shorelines were the first to be assessed. The T-sheets from OFR 93-574 (Thieler and Danforth, 1993) were cross-referenced with NOAA georeferenced T-sheets available through the NOAA Shoreline website (https://shoreline.noaa.gov/data/datasheets/t-sheets.html). Upon further investigation of NOAA's georeferenced T-sheets, original shoreline data, and Esri's World Imagery basemap all viewed together in ArcMap, there was a roughly 200-meter offset that existed in the georeferenced T-sheets from NOAA. This offset is likely due to a datum error. Most historical T-sheet surveys were conducted in the continental United States and reference the North American Datum of 1927. T-sheets in Puerto Rico, however, reference the Puerto Rico Datum of 1940. Therefore, we needed to georeference all T-sheets by using the non-georeferenced T-sheets.
While most of these T-sheets were already acquired from the internal data transfer, a few T-sheets were searched for and replaced based on their quality (e.g. part of the T-sheet may have been clipped) using the Non-georeferenced NOAA Shoreline Survey Scans website (https://nosimagery.noaa.gov/images/shoreline_surveys/survey_scans/NOAA_Shoreline_Survey_Scans.html). A total of 31 T-sheets were included from the original OFR and used to validate shorelines.

After all T-sheets were acquired, additional T-sheets adjacent to the original T-sheets were also acquired. To locate these additional T-sheets, an index to adjoining T-sheets on the bottom left of each T-sheet was used. The T-sheets were downloaded from the Non-georeferenced NOAA Shoreline Survey Scans website. 27 additional T-sheets were used to add to the existing shoreline data.
Some of these additional T-sheets from 1966 were only available as PiCture eXchange (PCX) image files, which are not compatible with ArcMap. These files were converted from PCX to TIF files using GIMP software (version 2.10.22, https://www.gimp.org/), an open-source graphics editor.

T-sheets dating from 1959 and beyond in this dataset were georeferenced using the NGS Coordinate Conversion and Transformation Tool (NCAT, https://geodesy.noaa.gov/NCAT/) and ArcMap. After opening up an ArcMap document, the data frame's coordinate system was set to NAD 1983 StatePlane Puerto Rico Virgin Isl FIPS 5200 (Meters). Each non-georeferenced T-sheet was added to the ArcMap document.
Eight well distributed control points were selected using the coordinates indicated along the edges of the T-sheets. These coordinates were transformed using the NCAT Tool from the Puerto Rico datum (PR40) to the North American Datum of 1983 (NAD83) Puerto Rico State Plane Coordinate System. After the eight control points were georegistered, adding or removing control points from the link table was considered to obtain an optimal root mean square error (RMSE). No less than six control points were used to georeference each of these T-sheets. It is important to note that the RMSE conforms to the National Map Accuracy Standards.

In this dataset, T-sheets from 1901 needed to be georeferenced using a different method than the relatively newer T-sheets. These T-sheets reference an older Puerto Rico datum that was not available for conversion through the NCAT tool. These older T-sheets include triangulation stations which can also be used in the georeferencing process.
Locations of these triangulation stations were determined by using a shapefile of NGS known geodetic marks in Puerto Rico (https://geodesy.noaa.gov/cgi-bin/sf_archive.prl). This shapefile was projected to NAD 1983 StatePlane Puerto Rico Virgin Isl FIPS 5200 (Meters). Then, the non-georeferenced T-sheets were georeferenced by matching the name of the triangulation station to the name of the geodetic mark ("NAME" field) in the shapefile's attribute table. If there were multiple geodetic marks by the same name, the fields "FIRST_RECV" and "LAST_RECV," which indicate the years that the geodetic mark was first observed and last recovered, respectively, were used to select the appropriate mark based on the time frame.

After all T-sheets were georeferenced, a visual inspection in ArcMap determined further edits that needed to be made to the shoreline features.
Some of the features contained errors as a result of the software that was used to digitize the original shorelines (Thieler and Danforth, 1993). There were a series of clustered vertices and rough, jagged edges within some shoreline features that required deletion. These were likely artifacts of the digitizing track mode that may have inadvertently added erroneous vertices not reflective of the shoreline.
There were also errors regarding the alignment of the shoreline features to the georeferenced T-sheets. In some cases, the digitized shoreline did not line up with the shoreline delineated on the T-sheets. In other cases, there were T-sheet shoreline features that were omitted from the digitized shoreline as they may have been cut off when stitching together overlapping T-sheets.
Based on this assessment, it was determined that the majority of these shorelines would be re-digitized. In a few instances, only the vertices from the digitized shorelines were adjusted to match the shoreline delineated on the T-sheet as there was little spatial discrepancy between the digitized shoreline and the georeferenced T-sheet.

The next shoreline features to be QA/QC'd were the aerial photo shorelines. Unlike the T-sheets, the control points for the aerial photos were established using orthorectified imagery downloaded from the NOAA Data Access Viewer (DAV, https://coast.noaa.gov/dataviewer/#/imagery/search/where:ID=394). In an effort to reduce the uncertainty values for the aerial photos, natural color orthophotos were used as opposed to the quadrangles mentioned in the OFR (Thieler and Danforth, 1993).
Like the T-sheets, the aerial photos were acquired through an internal data transfer. These aerial photos were not georeferenced but were separated by municipality and year. The aerial photos were georeferenced using ArcMap. In the ArcMap document, the data frame's coordinate system was set to match the shoreline data (NAD 1983 StatePlane Puerto Rico Virgin Isl FIPS 5200 (Meters)). In order to retrieve the corresponding natural color orthophoto, a few data sources were needed.
First, even though the quadrangles were not used directly in the georeferencing process, they were used to help locate the general vicinity of the aerial photos. The quadrangles were georeferenced using the same method as the T-sheets as the quadrangles contained latitude and longitude points in the Puerto Rico datum of 1940 (see Step 4).
Then, the tile index shapefile downloaded from the bulk download page on the NOAA DAV (https://chs.coast.noaa.gov/htdata/raster1/imagery/PuertoRico_RGB_2007_394/) was brought into the ArcMap document. The tiles from this shapefile contained the associated orthophoto direct download link. These along with the Esri basemap imagery helped determine the approximate location of the non-georeferenced aerial photos.
The aerial photos were georeferenced using consistent ground points available in the orthophotos (e.g. road intersections, man-made structures, rocks along the shoreline). At least four ground control points for a first order polynomial transformation and at least six points for a second order polynomial transformation were selected for each of the aerial photos to achieve an acceptable RMSE. It is important to note that the RMSE conforms to the National Map Accuracy Standards.

After all available aerial photos were georeferenced, a visual inspection in ArcMap determined further edits that needed to be made to the shoreline features.
Some of the features contained errors as a result of the software that was used to digitize these original shorelines (Thieler and Danforth, 1993). There were a series of clustered vertices and rough, jagged edges within some shoreline features that would require deletion. These were likely artifacts of the digitizing track mode that may have inadvertently added erroneous vertices not reflective of the shoreline.
Based on this assessment, it was determined that many of these shorelines would be re-digitized. In some instances, only the vertices from the digitized shorelines were adjusted to match the shoreline delineated on the aerial photo as there was little spatial discrepancy between the digitized shoreline and the georeferenced aerial photo.

Additional fields that were added to the attribute table needed to be calculated. All fields except "UNCY" were calculated using information provided in the OFR or directly from the source itself. The uncertainty calculations were largely adapted from the uncertainty equations created by Fletcher and others 2012.
The following summation in quadrature equation was used to calculate the uncertainty for the T-sheet shorelines:
Sqr(([U_SURVEY]^2) + (([U_CELL_SIZE] * [U_RMSE])^2) + ([U_CELL_SIZE]^2) + ([U_DIGITIZE]^2) + ([U_SHORE]^2))
where U_SURVEY = error associated with the T-sheet survey (10 meters for T-sheets before the year 1960, 3 meters after 1960; see Shalowitz, 1964), U_CELL_SIZE = the pixel cell size of each rectified T-sheet, U_RMSE = the unitless root mean square error value of each T-sheet recorded from the georeferencing link table, U_DIGITIZE = the maximum error associated with digitizing as represented in past studies (1 meter; see Hapke and others, 2010), U_SHORE = the estimated uncertainty in the shoreline position due to the fact that the high water line varies with tide. U_SHORE was estimated to be 3 meters and was derived by taking the average great diurnal tide range (or the difference between mean higher high water and mean lower low water) and dividing it by the average beach slope. The tide range values were taken from six well-distributed tide stations from NOAA Tides and Currents (https://tidesandcurrents.noaa.gov/). The average beach slope was calculated using the average slope value from the 2018 Lidar data. Dividing the great diurnal tide range (an approximation of the water level's vertical range) by the average beach slope (represented as rise/run) yields the horizontal range of where the water meets the beach. It is presumed that the shoreline proxies derived from historical sources will fall within this horizontal range, thus an appropriate value to apply to the uncertainty equation.
The following summation in quadrature equation was used to calculate the uncertainty for the aerial photo shorelines:
Sqr(([U_AIR_PHOTO] ^2) + ([U_BASEMAP] ^2) + (([U_CELL_SIZE]* [U_RMSE])^2) + ([U_CELL_SIZE]^2) + ([U_DIGITIZE]^2) + ([U_SHORE]^2))
where U_AIR_PHOTO = the average aerial photo error due to distortion and cartographic error (3 meters; see Crowell and others, 1991), U_BASEMAP = the horizontal accuracy from the orthophotos' metadata (5 meters), U_CELL_SIZE = the pixel cell size of each rectified aerial photo, U_RMSE = the unitless root mean square error value of each aerial photo recorded from the georeferencing link table, U_DIGITIZE = the maximum error associated with digitizing as represented in past studies (1 meter; see Hapke and others, 2010), U_SHORE = the estimated uncertainty in the shoreline position due to the fact that the wet/dry line varies with tide. U_SHORE was estimated to be 3 meters.
While all T-sheets were available to cross-reference with the original digitized shorelines, only a sample of the aerial photos were available for the QA/QC. In the event that there was not an aerial photo available for a digitized shoreline, the originally reported uncertainty value from the OFR (9.25 meters) was left as the uncertainty value.

After all the individual shoreline feature classes within the personal geodatabase were evaluated a final time, they were exported as shapefiles and merged together into one shapefile (ArcToolbox >> Data Management >> General >> Merge).

The data were projected to WGS 84 using ArcToolbox >> Data Management Tools >> Projections and Transformations >> Feature >> Project (Parameters: Input Coordinate System = NAD_1983_StatePlane_Puerto_Rico_Virgin_Isl_FIPS_5200; Output Coordinate System = GCS_WGS_1984; Geographic Transformation = Puerto_Rico_To_NAD_1983 + Puerto_Rico_To_WGS_1984_4).