Continuous and optimized 3-arcsecond elevation model for the United States west coast (32-bit GeoTiff, geographic, NAD83)

Step 1. Download CRM regions: Coastal Relief Model (CRM) data were downloaded in binary float format except for southern California, which was downloaded as a zip folder of .grd files. The southern California individual .grd files were merged in alphabetical order using the ArcGIS 10.3 Data Management\Raster\Raster Dataset\Mosaic to new raster tool to a 32-bit float GeoTiff using a mosaic method order of "LAST" and all other options as defaults. This process step, and all subsequent process steps, were performed by the same person: Kevin M. Befus.

Step 2. Download NED tiles: Individual 1-arcsecond National Elevation Dataset (NED) 1 degree by 1 degree tiles were downloaded within the extents of the CRM data. Each tile was downloaded and unzipped in Python 2.7.11 using the standard libraries zipfile and urllib2. For the U.S. east and west coast separately, the NED tiles were merged into a single NED DEM using the ArcGIS 10.3 Data Management\Raster\Raster Dataset\Mosaic to new raster tool to a 32-bit float GeoTiff using a mosaic method order of "LAST", an output cell size of 3-arcseconds, and all other options as defaults. No horizontal datum transformation was required for NED data.

Step 3. DEM inspection and artifact identification: Visual inspection of NED and CRM products revealed numerous artifacts associated with the creation of each dataset, primarily due to merging 1 degree x 1 degree subsets of the elevation data or due to interpolation in data-poor regions. Areas where grid merging created offsets in elevation were manually delineated with polygon features (DEM_artifacts.shp), indicating whether the merging problem was "terrestrial" (elevation greater than 0 meters), "offshore" (elevation of less than 0 meters), or "both". Other areas where data quality was questionable were also indicated with polygon features and assigned a type of "Merging", "Interp", or "Channel". Not all areas where artifacts were present were delineated. The primary focus of the artifact delineation was between -50 to 50 meters elevation. Only those interpolation artifacts that could be replaced by other data, were unavoidable/unfillable using other datasets, or were distinguishable from nearby bathymetry or topography were delineated and saved in DEM_artifacts.shp.

Step 4. Supplementary data download: Additional bathymetry data were downloaded to replace the cells where the NED and CRM data contained artifacts or had no bathymetry data. These data were downloaded from the Canadian Hydrographic Service. The source name of this step corresponds to the source identifier name and is the original downloaded dataset.

Step 5. Central and Northwestern Pacific region corrections: In the central_pacific_crm_v1.flt and nw_pacific_crm_v1.flt, coastal inlets and low elevation river channels with elevation artifacts were filled. First, the artifacts were identified as type = "channel" artifacts in DEM_artifacts.shp, although these artifacts were primarily associated with merging CRM 1 degree cells, and given a flag attribute to "fill" (Action="fill"). Then, the original CRM .flt's were individually filled using the ArcGIS 10.3 Spatial Analyst\Hydrology\Fill tool with default settings, yielding temporary rasters "nwp_filled" and "cp_filled". The features with the channel artifacts were then used to extract the filled elevation data by mask and were merged with the original CRM DEMs using the Extract by mask and Raster calculator merging procedure in Step 5 with the output extents set to the CRM .flt files.

Step 6. Add Vancouver bathymetry: To extend the bathymetric data smoothly into Canada on the west coast, an additional bathymetric dataset was required. The Vancouver data, 092G.bag, was converted to GeoTiff using Hydroffice.bag 0.2.15 BAGFile command to load the data and elevation.Elevation2Gdal to export the GeoTiff file ("CHS_Vancouver_500m_proj.tif"). The cooridinate system for CHS_Vancouver_500m_proj.tif was assigned with the ArcGIS 10.3 Data Management\Projections and Transformations\Define Projection tool with the Coordinate System specified as World_Mercator. CHS_Vancouver_500m_proj.tif was projected to the WGS_1984 geographic coordinate system as a temporary raster ("CHS_Vancouver_500m.tif") using ArcGIS 10.3 Data Management\Projections and Transformations\Raster\Project Raster tool using an output cell size of 3-arcseconds, BILINEAR interpolation for the resampling technique, and otherwise default values. CHS_Vancouver_500m.tif was first set to null above 0 meters elevation and then merged with nw_pacific_crm_vc.tif. CHS_Vancouver_500m.tif values were used directly for cells where nw_pacific_crm_vc.tif was null. The merge used the ArcGIS 10.3 Raster Calculator command "nw_pacific_crm_vc2.tif = Con(IsNull(SetNull(CHS_Vancouver_500m.tif greater than 0,CHS_Vancouver_500m.tif)),Con(IsNull(nw_pacific_crm_vc.tif),CHS_Vancouver_500m.tif,nw_pacific_crm_vc.tif),CHS_Vancouver_500m.tif)" with the cell size set to 1-arcsecond, the output coordinate system to the same as nw_pacific_crm_vc.tif, and the output extent to the union of inputs in the environmental settings. Raster Calculator automatically resampled CHS_Vancouver_500m.tif to match the output cell size. Low data density at the southern end of Vancouver Island in the nw_pacific_crm_vc2.tif led to unavoidable artifacts that persist to the final raster product, west_cdem_v1.tif.

Step 7. Bathymetry merge: After these region-specific corrections, the CRM region-based rasters were merged into two separate coastal datasets for the U.S. east and west coasts. The ArcGIS 10.3 Data Management\Raster\Raster Dataset\Mosaic to new raster tool was used to merge the CRM region-based rasters to temporary rasters ("east_CRM_WGS.tif","west_CRM_WGS.tif") with an output cell size of 3-arcseconds and a merge order of "FIRST". The order of the east coast regions was geographically in order along the coast from the northeast region to the western Gulf of Mexico. The order of the west coast regions was from south to north. The east_CRM_WGS.tif and west_CRM_WGS.tif rasters were transformed into the same horizontal datum as the NED data (NAD 1983 from WGS 1984) using the ArcGIS 10.3 Data Management\Projections and Transformations\Raster\Project Raster tool with a cell size of 1-arcsecond and otherwise default values.

Step 8. Correct CRM merging: The final step for the CRM bathymetry data was to smooth any merging artifacts present in the merged CRM DEMs. For CRM grid cells within features from DEM_artifacts.shp that indicated a merging error in an offshore area where the dataset had an elevation of less than or equal to 0 meters (type = "Merging", Loc="offshore" or "both"), the SciPy 0.17.0 function griddata was used in Python 2.7.11 to bilinearly interpolate over the erroneous values (method="linear"). Grid cells with centers inside a stitching error feature had their values set to NaN (not a number), and the grid cell locations were then used as "unknown" value locations to receive the interpolated values. PyShp 1.2.1, Shapely 1.5.12, GDAL 1.11.1 were used in Python 2.7.11 to select the grid cells corresponding to the stitching artifacts. For the interpolation step, stitching errors above 0 meters elevation in the CRM data were not changed unless the stitching feature contained the additional flag of "higher" (DEM_artifacts.shp, Action="higher"), which set the maximum interpolation elevation to 5 meters instead of 0 meters. This "higher" option was only used for the seam between the northwest Pacific region and the Vancouver bathymetry data.

Step 9. Merge elevation products: Merging the CRM and NED data required multiple steps to extract the highest quality data from each source while ensuring a smooth boundary across the continuum between bathymetric and topographic elevations, chosen for visual purposes as -10 to 10 meters. The individual steps for the merge are expanded as substeps. Two differences between the CRM and NED datasets led to the development of this multiple step merging procedure. First, NED data are always based on a vertical datum of NAVD 88, but NED data are nominally assigned a value of 0.0 meters elevation for anywhere with an elevation less than or equal to 0 meters. However, many anomalous, non-zero values persist for marine or other coastal waters due to the source datasets as well as other processing used in NED creation. Thus, only using a cutoff of 0 meters elevation to snap CRM bathymetry data into NED will result in anomolous "islands". Additionally, in some coastal waterbodies, NED may have positive elevations, whereas CRM contains bathymetry with elevations below sea level. Secondly, CRM data above 0 meters elevation are only whole numbers with no non-zero decimals and also appear smoothed, potentially due to a processing step in CRM creation. This rounding in CRM smears the interface between bathymetric and topographic data (i.e., 0 to 1 meters) and creates the problem that coastal areas with elevations below 0.5 meters may have been assigned to 0 meters. Thus, no single value in neither NED nor CRM allowed the extraction and insertion of one dataset into the other. Instead, CRM values equal to and above 0 meters elevation were assigned values from NED. Then, NED elevations less than 8 meters above sea level (0 meters, NAVD88) were replaced with the edited CRM values. This assigned bathymetric data in coastal areas where it was available in CRM without relying on only one of the datasets to specify the coastal boundary (e.g., a coastal pond with a water elevation of 3 meters in NED can be supplied with CRM bathymetry data so long as CRM elevation is less than 0 meters).

Step 9a. Fill NED data with CRM: CRM values were assigned where the NED data were null (i.e., in areas far offshore) using ArcGIS 10.3 and the Raster Calculator command "west_NED_filled.tif = Con(IsNull(west_NED.tif),west_CRM.tif,west_NED.tif)" with the output extent set to the union of inputs in Environment Settings.

Step 9b. Set CRM topography to null: CRM values were set to null values if they equaled or were above 0 meters of elevation using ArcGIS 10.3 and the Raster Calculator command "west_CRM_null.tif = SetNull(west_CRM.tif > 0 ,west_CRM.tif)" with default Environment Settings.

Step 9c. Fill CRM_null data with NED_filled: The filled NED elevations were used to assign elevation values to the nulled CRM resulting in a temporary raster with original CRM bathymetry below 0 meters elevation and NED elevations otherwise. This step used ArcGIS 10.3 and the Raster Calculator command "west_CRM_null_filled.tif = Con(IsNull(west_CRM_null.tif),west_NED_filled.tif,west_CRM_null.tif)" with the output extent set to the union of inputs in Environment Settings.

Step 9d. Merge adjusted CRM and NED datasets: The adjusted CRM and NED datasets were merged by extracting low elevation areas from the NED-filled CRM raster and higher elevations from the CRM-filled NED raster using ArcGIS 10.3 and the Raster Calculator command "west_cdem.tif = Con(west_NED_filled.tif < 8,west_CRM_null_filled.tif,west_NED_filled.tif)" with default Environment Settings.

Step 10. Extract coastal elevations: A polygonal mask was manually drawn around the extent of the planned coastal groundwater models ("dem_mask.shp"). The mask was extended offshore a considerable distance to allow the dataset to be used for other purposes. No consistent offshore extent was used (i.e., did not use a bathymetric contour or a buffer from shore). Instead, the mask was drawn a relatively consistent distance offshore that still contained bathymetric data. The terrestrial extent was chosen solely based on the potential groundwater model extents. The extraction was achieved by using ArcGIS 10.3 Spatial Analyst\Extraction\Extract by mask tool with default settings using dem_mask.shp as the mask. No data values in west_cdem_v1.tif and east_cdem_v1.tif represent both where the mask did not extract elevations and where the underlying datasets contained no data values or were not downloaded for a particular area. Both NED and CRM data outside of the mask can be joined to west_cdem_v1.tif and east_cdem_v1.tif using either the Raster Calculator conditional merging procedure in Step 5 or using the Mosaic tools in ArcGIS or software and environments capable of handling GeoTiff data.

Added keywords section with USGS persistent identifier as theme keyword.