What is BlueTopo?

BlueTopo is curated and compiled bathymetric source data that provide a definitive U.S. nationwide model of the seafloor and navigational lakes. It isn’t complete at the time of this writing, but NOAA Office of Coast Survey is adding to it regularly. While it isn’t meant for navigation, it is targeted at ocean science and engineering. One of the neat aspects compared to our usual topographic elevations is that it includes spatially variable information about the elevation uncertainty and spatial metadata about the sources of data used in the compilation. You can also view and query the data through nowCOAST.

What is it adding to DAV

BlueTopo adds a lot of bathymetry access that had been missing in the Data Access Viewer. The orange areas in the image below show where we had elevation data before adding BlueTopo, at least near the mainland U.S.

After adding the current BlueTopo data, we see a large part of the waters off the eastern side of the U.S. are covered in the image below.

As more data is added to BlueTopo, it will also be available through the DAV.

Roadblocks for DAV

You might think that BlueTopo was like any other digital elevation model (DEM) and wouldn’t be a problem in DAV. However, as alluded to above, instead of a single band of elevation the files are three bands with elevation, uncertainty, and source information respectively. All in floating point. While that isn’t inherently a problem for people to use the data, extra care has to be taken when you start reprojecting.

Source information band

The biggest issue is with the source information (or metadata) band. It contains values that are the index for a raster attribute table (RAT), such that for a given pixel you can get the information regarding which survey or other source was used to derive the elevation and uncertainty. For example, a pixel might have a value of 65091, for which the raster attribute table would tell you that in was interpolated from Office of Coast Survey’s survey H10674 in 1996, along with a few other details about the survey.

If we were simply clipping areas from the original tiles, there wouldn’t be an issue. When we reproject, there is always some sort of interpolation involved. For elevation data, we commonly use bilinear interpolation. Using something like nearest neighbor will often leave noticeable artifacts in the results. However, using bilinear interpolation on the metadata band would be a disaster. You might interpolate between a pixel derived from survey A and one derived from survey B and end up showing it was from survey C. Most likely you’d end up with a number that didn’t match anything in your RAT.

We have an additional issue if we change to a projection in feet, such as many of the state plane projections. This isn’t an issue once you’re beyond the bounds of the state plane system, but if you’re nearshore, DAV’s default projection for elevation is the state plane projections. Similar to the issue of interpolating above, applying a meters to feet conversion to the index won’t result in a number you can find in the RAT.

While a minor point, it’s worth noting that the source information/metadata band is only using integers for the indexing, but it is stored in floating point. Though it appears to work just fine in GIS software, having it as an integer band would make more sense.

Lastly, we’ll need to do something with the RAT, especially if multiple BlueTopo tiles are involved in one of our output tiles. Each input tile will have it’s own RAT that only contains the source information for that tile. Since we’re combining tiles, we need to make a RAT that has a combination of the sources, but without any that aren’t in our output tile.

Solutions

Our raster processing within DAV is done using GDAL. Fortunately, GDAL provides us with tools we’ll need to handle that metadata band differently. When reprojecting, we can easily specify the a list of the resampling methods we want to use, one per band. In this case, that’s [‘BILINEAR’,’BILINEAR’,’NEAREST’] if you’re using the GDAL python API for gdal.Warp.

We do a similar operation when we use gdal.Translate to convert the units and provide scaling for each band. In this case, we could convert the first two bands to international feet and leave the third (metadata) band alone by specifying the python API scaleParams option as

[ [0,0.3048,0,1], [0,0.3048,0,1], [0,1,0,1] ]

If the user is also changing the vertical datum, for instance to mean lower low water (MLLW), we have to take care that we only modify the elevation and not the uncertainty and metadata bands. That generally means separating the bands into separate files and doing raster math with gdal_calc as appropriate since the conversions are spatially variable, then recombining using gdal.Translate again.

For the RATs, we need to do a couple of steps. First, we need to build a RAT that contains all the unique rows from the, potentially, multiple input RATs. We do this by starting with the first RAT and then looping over the rows of the other RATs and adding any row with an index we haven’t already seen. The second step is to find all the unique values in the new tile we’ve created. There isn’t a way in the GDAL API to remove rows from a RAT, but we can clone the structure of a RAT and build the data by only adding the rows we want to keep, the ones that match our unique values in the raster.

Other Changes

Three-band to Single-band

We did some other changes to the files while we were at it. Even though it was a minor point, I wanted to change the floating point values in the metadata band to integer since all the RAT index values are integers. However, GDAL doesn’t let you create a tiff file with mixed band types. Combining that with my perception that most people are looking for their DEMs to be a single floating point band and will be a little confused by the three-band rasters, we decided to split the bands into individual rasters. Instead of getting our usual DAV naming that would look like bluetopo_J{jobid}tR0_C0.tif as a three-band floating point raster, you’ll bluetopo_J{jobid}tR0_C0_elev.tif, bluetopo_J{jobid}tR0_C0_uncert.tif, and bluetopo_J{jobid}tR0_C0_meta.tif, which are single-band floating-point, floating-point, and integer rasters respectively.

Resolution

The BlueTopo data has multiple horizontal resolutions, varying from 2 meter cells nearshore to 16 meter cells offshore. The DAV system creates a new set of tiles based on the bounds selected by the user. To keep things reasonably consistent, it outputs those new tiles using the highest resolution of the input data within the area. For selections that cross resolution boundaries, that means that some areas will have a higher resolution than they should, effectively up-sampled, or as I like to call it – geospatial fantasizing.

Loose Ends

There are a few things that probably still need work. The DAV size estimates that restrict how large an area you can grab are coming in much larger than the actual outputs. There are a couple of reasons behind this. First, the data compresses well and the estimates don’t take that into account. The metadata band compresses really well as there are probably only a few to a few dozen unique values in it. Second, the resolution mentioned above. The DAV system is assuming 2-meter data everywhere for BlueTopo. If you’re in the 16-meter resolution area, it’s already 64 times too big with the estimate. Don’t be surprised if you see DAV estimate a 10 Gb output but the zip file you receive is only 10’s of Mb.

The metadata band doesn’t have a nodata value, but it may have a legitimate value of zero (the code for National Bathymetric Source Generalization). When reprojected, you end up with nodata spaces along the edges filled with zeros.

The RATs work well in QGIS and are set to use the source_survey_id field as the legend label. They also work in ArcPro, though the field chosen as the legend label seems a bit random to me. Global Mapper doesn’t appear to read the *.aux.xml file containing the RAT at all, so we may need to find a workaround for that. If you have other software packages, feel free to let us know about the RAT functionality in the comments.

Feedback?

We’re very curious to know if splitting the bands into separate files was a good thing or a bad thing. While they can certainly be merged back together, that might require making the metadata band floating-point again. If you’ve got some thoughts on better ways to serve the data, let us know in the comments. Be sure to check out the BlueTopo web site for other ways to access the data too.

(Image credit to Climate Central)

Why don’t we provide maps / map layers up to this amount in our Sea Level Rise Viewer?

Our viewer focuses on the next 75 years (out to 2100) and uses the latest Sea Level science that was published in the 5th National Climate Assessment and the 2022 Interagency Sea Level Rise Technical Report.  This video provides a good summary of that report.  https://coast.noaa.gov/elearning/video/slr/takeaways/

The maximum amount of sea level rise (SLR) projected in the U.S is 8.86 feet (Grand Isle, LA) under a high (worse case) scenario which has less than a 1% chance of happening by 2100. So we feel our current mapping to 10 feet Mean Higher High Water (MHHW) by 2100 is sufficient. And, we now map in half-foot increments, which provides increased accruacy based on the increased vertical accuracy of the lidar-based elevation data now being collected and used.

If we decided to project out to 2150 based on the Tech Report, the range (2.6 -12.8 feet) between SLR scenarios would be almost 10 feet, which we feel is too uncertain for planning use.

If users want to see inundation impacts greater than 10 feet (MHHW) for events like hurricane storm surge or tsunamis, it is recommended they use other tools and data, such as the National Storm Surge Risk Maps provided by the National Hurricane Center, or Tsunami Aware tool which shows tsunami evacuation zones.

Perhaps the most compelling reason not to map up to 216 feet is that we currently map at half-foot increments, therefore provide more than 20 inundation map layers and services for download. If we mapped to 216 feet we would have to maintain and provide 432 map layers. This is not practical given our current 3 meter horizonal mapping resolution. We also update DEMs as new data becomes available and remap the layers. That amount of maintenance effort would also be impractical.

Finally, we only provide mapping for the U.S. States and Territories as our NOAA mandate does not extend to international locations. Climate Central provides maps and tools that show up to 30 feet of inundation and provide it internationally. We provide them with DEMs for mapping in the U.S., but if international people ask us for mapping assistance, we usually refer them to Climate Central.

Edit July 18, 2025: If you’re having trouble since July 16, 2025, make sure you clear your cache or otherwise do a hard reload of the Data Access Viewer application. The back-end changes are such that the old code can’t search properly.

Front-end Changes

The main change you’ll see on the front-end is during the checkout process. There is now a little box near the top of the options page to enter a ‘label’. It looks like this:

What’s it for? It lets you add a prefix to all the filenames that come with your data. So, if you add a label like ‘my house’ and the output from your job would have been sc2017_cst_counties_Job1215633.zip, it will now be my_house_sc2017_cst_counties_J1215633.zip. Note that the space in the label was turned into an underscore. That will happen to any characters that wouldn’t be particularly filename friendly. That same label will also be applied to the files inside the zip. Providing a label is completely optional and is intended to help you with your file management.

There are a few other front-end changes that may be noticeable. The colored footprints you see at the start or after doing a search are now done as vector tiles and other parts of a search are also done with a new system. You may notice a speedier response. However, you may also notice that the colors aren’t getting more transparent as you zoom in. That should be a temporary issue. The whole front-end UI is going to be rewritten now that the back-end is done and there will be several other enhancements in it.

Admin Interface

The bigger changes right now are in what we call the ‘admin interface’ and it’s all about job management. This is where you’ll go to see the status of the data you’ve requested. Or to download the data because you deleted that email with the link. Or to rerun a job because you forgot to download it before the output was cleaned off.

First, some nomenclature. If you think of the cart in the front-end, you probably realize that you can put multiple things in that cart, each of those things being a dataset with a bit of geography specified. We call everything in the cart a request and each of the things in the cart a job. You’ll see that in the admin interface.

As usual, you’ll get an email when you make a request that confirms that we received it. There is a link in that email and it will take you to the admin interface. You’ll need to verify your email first. That’s just one of those typical widgets that sends an email to the address you enter and you click the link in that email. It opens a new tab in your browser with the admin interface page. You only need to verify once every 30 days. The page you come to looks like the image below.

The particular request shown is my test suite of fairly random things, but it works to illustrate a few things. At the top, you see overall information such as the request number, a map with the area selected for each job in the request highlighted, and a progress bar. The progress bar is counting the fraction of the jobs that have completed. It will turn red if there is a failed job. In this example, two of my seven jobs have already finished.

Below that summary is a table of the jobs within the request, including job identifier number, the type of data, the title of the data, the label you might have given it, and its current status. You’ll also see three horizontal dots at the end that give a menu for that job. We’ll get to that in a bit. First, let’s look at more detail about a job. Those job identifiers are hyperlinked, so you just select the one you want. That will bring you to a screen like the image below.

This is a job still under processing and along the top we can see a bit of information about what it has done so far. The stage where we setup the processing has a checkmark, so that’s done. It’s currently working on something labeled as ‘Clip Raster’ and it’s on step 4 of 11 of that. Behind the scenes, we’ve cut the job up into manageable tiles and spread the processing to multiple machines to handle in parallel. In this case, there are 11 tiles and 4 of them are already done. When all 11 are done, it will move onto the next step and we would see that step shown. The next step would be assembling all those tiles together in a zip file and sending you the notice that it is ready for you to download.

Under the status information there is a set of tabs that let you change the information shown. The options tab is open in the screenshot above and shows the options you have chosen. If you see that you’ve chosen something wrong, you can change some of those options and restart the job. It will start a new job associated with request and mark the old job as superseded.

The next tab is the AOI or area of interest, as shown in the screenshot below. It gives you a map with the area you selected highlighted in red. Currently, that’s always a rectangle. The new front-end will also allow polygon selection. You can’t change the AOI, but there is a link in the lower right corner that will take you back to the DAV with this dataset and AOI so you modify it and make a new request.

You may notice a couple of other tabs in the screenshots. The one for logs won’t show up for users, so you can ignore that one. The last one is for issues. If a job has a failure for some reason, we send you an email right away. The old system didn’t send anything until a human took a look. The issues tab will have the same information that was sent to you for a failed job. It does try to figure out what went wrong, though that system has limits. It may still require a human to take a look.

Lastly, lets take a look at that little menu for each job on the requests page. Yours won’t look exactly like this, but it’s close. In the image below, I’ve clicked on the three dots for the job I labeled ‘to_wgs84_ellip’ and it has opened a little menu of job actions.

The ‘View’ option simply takes you to the job details, which is the same thing as clicking the Job ID link on the other side of the table. The ‘Restart Job’ allows you to restart a job in the case that the output data has been deleted. The current cycle is to delete after 15 days. The ‘Download’ selection begins a download of the data, which is the same thing as clicking the link in the email you received when the job finished. It does also tell you the size of the zip file so you have some idea how much your downloading first. Lastly, the ‘Open in Data Access Viewer’ selection takes you to the DAV with the dataset and AOI set.

I hope you’ve found this quick tour useful. I expect there will be changes as we start to hear from real users like yourself.

Motivation

Let’s start with our motivation for this change. It’s not money. The current system runs on old cast-off machines that are long out of maintenance. There are costs to keep things patched and replace stuff when it dies, but it’s generally not too expensive. However, the server room is also on the first floor about 100 yards from, and maybe five feet above, a river by the coast in a hurricane zone. Moving the compute system and the data to somewhere else makes sense from a continuity of operations perspective. That is what drove the decision. There may have been other operating options besides the cloud, but NOAA has a “cloud first” directive, so that’s where we’re going.

Lift and shift or cloud native?

As we already have a system that works, we could just do a “lift and shift”. We’d reserve some machines in the cloud to replicate what we have on-premises and simply install the system there. Certainly feasible, but it wouldn’t take advantage of what the cloud has to offer. It’s also the case that the current DAV system has grown organically (haphazardly?) over 25 years and has a single retirement age maintainer that knows the system. Rebuilding the system as cloud native and more fully incorporating newer software offerings is the approach that made the most sense to us.

Cloud native can probably mean a number of things, but there are a couple that I’m particularly focused on. First, the system can expand and contract resources such that I only use, and pay for, the resources that I need at any given time. Second, the system is containerized such that it can easily be redeployed on different hardware. In our case, we’ll be using Docker containers that are managed by Kubernetes to do the scaling of resources.

Beta Testing

Before getting into more details of how the new system works, I should mention that we’ve set up a beta test so you can try it out. We’re mirroring all the requests that come into the current system such that they also go to the new system. This has let us spot and solve lots of problems with real requests for things I would never do. We’ve added an opportunity to opt-in during the checkout process of the current front-end so that you’ll get your outputs from both the old and the new system. The email from the new system will have a link to a feedback form so you can tell us what you think. It will also have a link to opt-out.

What has changed?

Parallel processing

The current system devotes one CPU to one job. It splits bigger jobs into tiles, but the job still only gets one CPU and does the tiles one at a time. The machines have multiple cores, so a given machine might be doing multiple jobs at once. This setup prevents a big job from clogging the system and getting in the way of small jobs.

The new system embraces the cloud advantages and expands the system as needed. The individual virtual machines are relatively small, but we spin-up more machines as needed and do multiple tiles at the same time. If a job had 100 tiles, there could be 100 virtual machines working on those tiles for a few minutes and then the machines would be scaled down. The next step in the process might require only one machine, to zip everything up for example, so we reduce the resources to only what we need (see Figure 1 below). The scaling up and down is being done using Kubernetes with Docker containers. We do have a limit to have far it will scale. This is to reduce our risk of something running wild and so we don’t cause our own networking bottlenecks. In addition, we’ve added custom scaling rules with Keda to control the resources based on the job requests instead of CPU or memory usage.

How do we handle the big job getting in the way of small jobs? One way is by letting the system scale to much larger than our current system, but that’s not really sufficient. In addition, we’re using BullMQ queue groups with a large group and a small group. As resources become available, tiles to process are chosen by alternating between the two groups. The large group has jobs with several tiles. The small group has jobs with only a few tiles. Even if a job with a lot of tiles came in first, a tile for a small job can jump in front of the those tiles.

Toolkit

The current system has grown bit by bit over the 25 years and the code-base reflects choices made at the beginning. It uses Perl as an orchestrator language that calls programs written in C or C++. Over time, open source options became available and replaced C code we had written, though not all of it. The primary packages pulled in include GDAL, PDAL, and LAStools (the open part). In approximately 2010, the job control system was rewritten in Java. It has always run on a *nix system, sometimes on more than one flavor at the same time, including IRIX, SunOS, Tru64 UNIX, and Linux.

The new system is using Python as the orchestrator language for processing and JavaScript as the orchestrator language for job control. Everything is going into Linux-based containers and we’re avoiding adding code we need to compile. We rely on PDAL, GDAL, and their dependencies such as PROJ, either calling them as a sub-process or with Python hooks. Those programs are doing all the heavy lifting. This means that all the outputs you see in the DAV can be created with software available to anyone without cost.

Contours

One of the tentative changes is the process for making contours. Contour generation has gone through a few iterations over time. I believe our first method was to create a raster from the average value per cell and use GMT program pscontour to generate the contours. At some point we switched to the commercial LAStools program las2iso, which makes the contours from a TIN surface of the points. This allowed us to make contours over void areas, but also makes contours that are very accurate to the data, which means very jagged. Taking the view that contours are generally for human interpretation, we did want smoother contours. We switched to use tensioned spline interpolation (using ANUDEM) to make a raster and then GDAL’s gdal_contour to generate the contours. Spline interpolation isn’t recommended for dense data like lidar as it tends to smear things a bit, but that gave us the nice smooth lines that make visual sense.

Unfortunately, the latest version of ANUDEM did not want to run in our Docker container. Looking for another way to make a smoother raster, we are trying inverse distance weighting (IDW) as our interpolation method and a coarser grid. At this point, feedback from the beta testing is conflicted. Some people absolutely want the smoother lines back. Others are liking that newer approach that can show features that the splines smoothed over. We may end up adding a choice to the new DAV related to how smooth you want your contours. An older version of ANUDEM does work in the container, though it is a licensed product.

I put together a little comparison of some of the possible options. The only difference is the method used to make a raster. In all cases gdal_contour was used to extract the contours. The three cases are:

• 1) IDW interpolation 3-foot resolution raster;
• 2) IDW interpolation 6-foot resolution raster; and
• 3) Spline interpolation 6-foot resolution raster.

If we look at one-foot contours on a small hilly area in South Carolina, we can see some difference in these methods.

Several features stand out to me for these contours. While all get the general features of the landscape, the least smoothed contours do a better job capturing the ridge in the middle. However, they also look unnaturally jagged in the flatter areas, such as in the upper right quadrant of the images. In those areas, the spline looks most like what we expect in nature. The coarse IDW grid is in between these two, possibly making nobody happy. That’s the one the beta system is using as I type this, but maybe that needs to change. Giving some options to the user may be a good idea, though every option adds more confusion for novice users.

What about the user interface?

We haven’t started the work on updating the front-end user interface that you interact with. We usually get positive feedback about it and the new back-end obviously works with the current front-end, so why would we change it? Well, there are definitely a few things that could use attention on the front-end, but the bigger drivers are enabling the passing of some new features to the back-end and an update of the technology used to enhance maintainability.

The feature that I think will be appreciated by the most people will be the ability to draw a polygon instead of just a box. I anticipate that mechanisms similar to those available in the United States Interagency Elevation Inventory will employed. That includes drawing the polygon on the screen or providing a GeoJSON string in the search box. It probably won’t include uploading a shapefile.

Another feature that I think people will like is providing a custom label for each job. This would be optional, but it would let you add some text to all the filenames associated with a job so it will be easier to remember what it is. Instead of a name like sc2020_sav_pd_Job22286, you could turn that into Bobs_house_sc2020_sav_pd_Job22286.

There will be some other changes coming too, such as a better way to see and interact with all your data requests. That should allow you to restart an old job if the results have already been removed instead of redoing it or asking us, or simply see what you’ve requested in the past.

Please try out the beta system and provide some feedback if there’s something you want to let us know. The feedback form doesn’t point you in any direction, it simply has a box for feedback. Whatever you want to let us know is fine. If we need some clarification, we’ll email you, but otherwise you won’t hear from us unless you request contact.

Background on Our Cycle and Specifications

The Sea Level Rise Viewer (Viewer) is a web mapping tool on NOAA’s Digital Coast that shows the impact of coastal flooding and sea level rise at a community level. To do so, it uses high-accuracy elevation data from lidar surveys and sea surface information from NOAA’s tide gauges and VDatum software. The Viewer is continuously updated with new data as they become available to reflect changes in the dynamic coastal zone.

When selecting a region to update, we consider the age and resolution of the existing DEMs underlying the Viewer, along with the availability of new data. Most areas are mapped at a 3-meter horizontal resolution, but a few states are still at 5-meter, which dates back to the initial version of the tool (GA was at 5-meter until this update). We also work with the National Hurricane Center (NHC), as they use our DEMs to model storm surge in the various oceanic basins. With an unfortunately busy several years, Florida was in need of updated modeling. 

When incorporating new data, our approach is to refine, striving to make the data more accurate, or ‘less wrong’. In recent years, the quality of contractor-provided breaklines and DEMs has improved significantly, capturing and flattening even small lakes, rivers, and streams, which simplifies our work. However, we still thoroughly review the data and make necessary edits to ensure we capture the potential inundation pathways from sea level rise. These edits may involve ‘burning-in’ culverts that allow water to flow beneath roadways, smoothing interpolation artifacts after nearshore buildings have been removed, and incorporating flood control structures that may be too narrow to be accurately represented in a 3-meter DEM, or otherwise weren’t captured by the lidar. These were all relevant, particularly in South Florida, during this update.

New Collections Become Available

Peninsular:

In Florida, a significant lidar collection effort took place in 2018-19 across an expansive 34,873 square miles. ‘The Florida Peninsular’ collect was supported by FDEM collaborating with the USGS, the Florida Department of Transportation (FDOT) and representatives from all five Florida Water Management Districts.

The specifications for the project aimed to meet or surpass the USGS Quality Level 1 (QL1) standards, with a nominal density of 10 pulses per square meter along with correspondingly high vertical accuracy standards. To collect the data, 11 lidar equipped aircraft were deployed nearly simultaneously across the state. 

High-quality elevation data is vital for inundation mapping, enabling identification of subtle terrain changes and vulnerable regions. It allows for accurate delineation of potential inundation areas, considering factors like slopes, drainage patterns, and barriers.

In addition to meeting the QL1 standards, the project included extra features to make it more useful. Mapping potential flooding from rising sea levels requires understanding how water bodies and low-lying areas are connected. To address this, the project added “connectors” to the breakline vectors. These connectors show areas that may be linked hydrologically, even if they appear disconnected from the airborne lidar due to obstacles like roads. The connectors were designed to seamlessly align with other water features, such as ocean, river, and lake polygons. After review, we were able to assimilate them into the DEMs through semi-automated processes – capturing many more flooding pathways than would have been possible otherwise.

Everglades:

Accurate lidar surveys in dense marsh areas like the Florida Everglades pose unique challenges. Thick vegetation or water inhibit signal returns, leaving void areas. The marshy topography, with subtle elevation changes, further complicates data collection and processing – identifying real ground returns from dense vegetation is no small feat, especially with a dearth of traditional survey marks. This first-ever collection was therefore a noteworthy accomplishment. 

Led by the USGS and National Park Service, data collection took place during a period of historically low water levels within the park. The contract team used both a topographic and a topobathymetric lidar sensor to collect terrestrial near-infrared, and submerged green wavelength data. Extensive ground surveys were conducted prior to data collection for testing and calibration (see the Project Report for impressive photos from this campaign). There was even a comparison between the topobathy raster elevations and a plumb-bob survey conducted across the extents of the northern park reaches – underscoring the difficulty of capturing traditional ground control points in this environment.

Final Pieces:

Rounding out the southern reaches of the state, NGS collected topobathy lidar throughout the entire Keys as part of a post-Hurricane Irma supplemental. They also flew Biscayne Bay and the Dry Tortugas. The latter, like the Everglades, is making its first appearance in the Viewer with this round of updates.

Combining These, Then Breaking Them Apart

So this sounds like just a few datasets that needed to be merged, right? Well, the Peninsular collection was released per county (>30), the Everglades were in at least three pieces, while the Miami-Dade area had been part of an altogether different collection, plus all of that NGS data. All told, there were 41 different source DEMs that were wrangled in order to paint one complete picture of the peninsula and Keys. This many disparate datasets creates challenges for data management, assimilation and subsequent mapping. Mapping, at this point, seemed like the easy part!  For, along with the data, we needed all of those breaklines, connectors and other ancillary metadata to do things like pull out and flatten all the ocean connected water to one value, and to establish and map confidence in the final products. 

We combined these data, ensuring consistent projections and geoids, then clipped and shipped these annoyingly large datasets into manageable chunks, as no one wants to download a single, 50 GB DEM.  

Determining where to split is a bit of a balancing act. For our internal, mapping purposes, bigger is generally better – to a point. Bigger input DEMs result in fewer output layers to manage through all of the subsequent steps, culminating in the final Viewer and downloadable products. But as the size increases, so do the amount of CPU cycles to make it through one area. Loathe to find a mistake after a 24 hour script run! 

The other downside of huge DEMs hinges on you, the consumer.  We try to keep DEMs to somewhere around 5 GB, that way it and any derived products, like inundation depth grids, don’t get unmanageable to pull down for further analysis. To keep things relatively clear, we also look to break things at county boundaries, and, if we are able, those that are demarcated by a river. That way any ocean connected water is easy to piece back together for a seamless look within the Viewer.

Mapping Products

While we have our mapping process published, I’ll summarize. We need two critical pieces to start from. The first, detailed here, is the underlying elevation data. The second, equally integral portion, is a variable surface reflecting local sea heights – a topic for another time. 

But there are a few new additions with this round of updates to share. First, we have started to map ½ foot increments, in addition to the whole feet currently in the Viewer. For now, these are only available for download in the regions updated with this cycle – that is, the majority of Florida, Georgia, CNMI, Guam and Catalina, CA. We are working hard to map these additional increments throughout the rest of the country and territories. Once the mapping is done, we will create the map services, modify the Viewer and roll it all out to production. Stay tuned.

We have also taken a cue from projects like the NCEI CUDEM and started to create spatial metadata files to go along with the DEMs. After all, with 41 different source datasets for the Florida peninsula alone, knowing where the data came from is not that straightforward. These are currently available through the Data Access Viewer for the aforementioned regions or by request. We will soon be updating our download page, where they will also have a home.

Other Relevant Links

• Sea Level Rise Viewer, Florida Peninsula Update
• SLR Data Download, including ½ foot increments for FL, GA, CNMI, Guam, Catalina Island
• Articles about the Peninsular Collect:
  • XYHT
  • Lidar Magazine

As federal agencies and others have tried to use the inventory to find gaps in coverage, the simplified footprints have become a problem. They can incorrectly show a sliver of a gap where none exists or vice-versa. This can waste a lot of time and effort as everything has to be double checked.

An additional problem is that the footprints can change during the process of data acquisition and there isn’t an indication of what the footprint represents. One footprint might be planned, but it might be changed during the contracting phase when you find out just how far your money can stretch. It might change again with the final delivered data if, for example, persistent cloud cover made it impossible to collect an area.

To address these issues, we’re making a change to what we store and serve with the inventory. The app will still use the simplified footprints for viewing and won’t look any different. However, we will have a second spatial table that has the detailed footprints. This table will be used when generating the GDB and GeoPackage download files. That means the download data will be a little different than the app. This change will happen slowly. The detailed footprint table is initialized with the simplified footprints, but we’ll start to change them to more detailed footprints during updates.

That leads us back to the issue of knowing what the footprint you’re looking at represents. The download GIS files will have an additional column in the NationalElevationInventory table for the FootprintType. This is just a text field and will contain entries like “Planned Footprint”, “Contracted Footprint”, “Completed Footprint”, “Simplified Web Footprint”, and “Unknown”. I’m hoping to never use “Unknown”. Other types may be added as we find the need. As you probably surmised, they are almost all “Simplified Web Footprint” as I write this.

Extra details/weeds

Now that the real news is out of the way, I wanted to provide some details on how all this came about and why it didn’t happen before. The addition of the detailed footprints was an unexpected opportunity that resulted from something else. While the inventory is held at NOAA, we send a replica of the database to USGS to get their updates. The database is managed by ArcGIS, which allows us to create the replica and sync the data changes made by USGS back into the primary enterprise geodatabase.

One issue we found is that updating a footprint in the database was not as simple as it should be. We couldn’t figure out a way in ArcPro or ArcMap to take the feature in a shapefile and use it to replace the feature in a specific row of the database using the GUI. Instead, staff would add a new row with the new shape, copy all the information from the old row, and delete the old row. As we expanded from a single table to a relational model to accommodate multiple entries for a few things, these updates became much more painful and time consuming. I’ll take the blame for causing that.

While there may be a way to do that update more efficiently in the GUI and we just missed it, it looked pretty straightforward to do it using arcpy and a toolbox. The basic approach is to reproject the new feature class if necessary, simplify it for the web, extract the shape with an arcpy.da.SearchCursor, and then use an arcpy.da.UpdateCursor to replace the shape for the correct row in the database. No more copying of all those attributes. There are some weird requirements though. In our setup, you have to be in an edit session for that UpdateCursor and if you have the attribute table open you open an edit session. The basics of the update code are shown below.

 # pull the shape we need from the feature class with the shape to update
simp_fc = "simplifed_dissolved_shape" # really our output from dissolve and simplify processes not shown here
shps=[] # a list, but there should only be one since we flattened
with arcpy.da.SearchCursor(simp_fc,["SHAPE@"]) as s_cursor:
    for row in s_cursor:
        shps.append(row[0])
shp = shps[0]

# whereclause has the condition to pick only the records we want
# database_table is the table to be updated.
# UpdateCursor is a wrapper class for arcpy.da.UpdateCursor
# adapted from the solution answer at https://community.esri.com/t5/python-questions/versioned-data-multiple-updatecursors-runtimeerror/td-p/1043525
# we have a versioned table, so we need to have an edit session. This makes sure we do.
with UpdateCursor(database_table,['SHAPE@'], whereclause) as u_cursor:
    for row in u_cursor:
        row[0] = shp
        u_cursor.updateRow(row)

Once I had a toolbox to do the updates, it made sense to add tools to create or delete a record in general. This gives us most of the pieces of CRUD. Now the tasks of simplifying and reprojecting could be done in the tool and make the process easier for staff. It could also provide more consistency to the process.

Detailed Footprints

How does this lead us to detailed footprints? We certainly could have made a separate spatial table for the detailed footprints before. However, we would have had to make sure everything stayed synchronized between the tables in a manual process. That appeared to be a recipe for disaster. Once we had a toolbox, and therefore a scripted consistent process, we could take the new or updated shape and put it in the detailed footprint table before simplifying it for the web app. Thus, the opportunity to finally capture the full detail was there and we’re going down that path. New entries will automatically get this treatment, but we’ll have to update the older records a bit at a time.

One last thing on the footprints. You might be wondering how the tables were joined so that the GDB and GeoPackage have the detailed footprints and the attributes from the original simplified table. While it would make sense to use the AddJoin function in arcpy, this resulted in field names that included the origin table. The solution was to do the join in the database to make a view with the desired fields and names, and then export (arcpy.conversion.FeatureClassToFeatureClass) into the GDB.

What’s a VRT?

I admit that I thought that VRT stood for virtual raster table, but I can’t find proof of that. It’s currently called the GDAL Virtual Format and can be used for more than rasters. VRT is the short name and the file extension. A VRT, is essentially a simple XML catalog for some raster files. It has information about the data projection, bands, and each file’s bounds. It could also have information about operations to be done, such as reprojection, but we won’t touch on that much here. VRTs are part of the Geospatial Abstraction Library (GDAL), which supports a great many GIS and remote sensing software packages. You can learn more about the details of a VRT on the GDAL Virtual Format page.

For our purposes, a VRT provides a way to reference a bunch of files as if they are one giant file. There are some some restrictions to what we can put in a VRT that follow logically if we think of it as one big image. for instance, you can’t put in files with different projections or a differing number of bands, just as you couldn’t do that with a single image. However, we can take thousands of image tiles that make up dataset, put them into a VRT, and treat it as a single image.

How do I use a VRT?

Let’s start with a simple example using QGIS (a free open source GIS package) to view some data. If you open this link to the VRT file I’ll use for some 2009 imagery of Kings Bay, GA, you’ll see lots of text referencing the individual image tiles. In QGIS, you select the Raster data type in the Data Source Manager, set the source type to a protocol, and paste in the URI for the VRT.

After adding it, it will show up in the map and you can use it like any other image, except you can’t write to it for obvious reasons. This example isn’t a particularly large one as there’s only about two gigabytes of data there. However, it would have taken a lot longer to download a couple of gigabytes compared to the time it took to show up. Zooming in and out works relatively quickly. It’s half meter data, so there is only so far you can zoom in. Part of the reason it can zoom well is because the underlying data that the VRT references is in cloud optimized GeoTiff format, so it only pulls what it needs over the wire.

Where do I find these VRTs?

Both NOAA and USGS have been putting up VRTs lately. I’ll review them separately. You can be sure there are a lot of other sources of VRTs too.

NOAA

I’ll start with source I have some control over, the data on the NOAA Digital Coast. If you search for data on the Data Access Viewer, you’ll find a link to ‘Bulk Download’ for each data set. For the data in raster format (i.e. not the lidar point clouds), that bulk download page will have a link to the VRT.

In addition to the Digital Coast, the National Geodetic Survey has created VRT files for their emergency response imagery. These are primarily post-hurricane imagery. The data and the VRTs are under https://noaa-eri-pds.s3.amazonaws.com/index.html. It takes a little bit for the page to populate and then it will show links to the various datasets. Each dataset tends to be separated into days and the flights in a day, with a VRT for each flight. The VRTs are in with the rest of the data, so use the search in the upper right to find them by entering ‘vrt’. You could also search for vrt on the opening page to see all the VRTs. Note that only the more recent datasets (maybe 2020 and later) have VRTs.

USGS

The US Geological Survey currently has VRTs for the 1/3 arc-second, 1 arc-second, and 2 arc-second seamless elevation data sets. These datasets cover the nation, so that’s a lot of territory in one file. If there is demand, I wouldn’t be surprised if they add VRTs for more data too.

Beyond Viewing

While it’s great that you can use the VRT in QGIS to view the data, what if you want to do more than viewing? You can certainly do that within QGIS where you can treat it like any other dataset. You can also use the applications that come with GDAL to do things like reproject and subset the data. Suppose I have a vector with the outline of South Carolina and I want to extract data from the USGS 1/3 arc-second elevation for South Carolina and reprojected it to State Plane coordinates. I can do something like this:

gdalwarp -cutline sc_outline.shp -crop_to_cutline -t_srs EPSG:2273 /vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/USGS_Seamless_DEM_13.vrt sc_dem13.tif

Notice that I put /vsicurl/ in front of the link to the VRT file. That tells GDAL what driver to use to access the data. The -t_srs EPSG:2273 specifies that the target reference system is South Carolina State Plane feet. The vertical units are still in meters though. You could do a two step process where you make a VRT as output from the gdalwarp and then use gdal_translate to change the vertical units.

gdalwarp -cutline sc_outline.shp -crop_to_cutline -t_srs EPSG:2273 /vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/USGS_Seamless_DEM_13.vrt sc_dem13_vertmeters.vrt
gdal_translate -scale 0 0.3048 0 1 sc_dem13_vertmeters.vrt sc_dem13_vertfeet.tif

The intermediate VRT file is very small and fast to create. Naturally, it’s going to have to download all the raster cells that fall in South Carolina, but it doesn’t have to download the full dataset and I don’t have to figure out which tiles I need for South Carolina from The National Map.

What about ArcGIS?

Many people use some form of ArcGIS to work with geospatial data, so how well do VRTs work there? The ArcPro documentation shows that a VRT will work as it is listed in the raster file formats. However, it has to be a local file as far as I can tell. Trying to give it a URL, with or without the /vsicurl/, does not work. If I missed the trick to using a URL, please tell us in the comments. However, you could make a local VRT that just points to the web VRT. For instance:

gdalbuildvrt noaa.vrt /vsicurl/https://chs.coast.noaa.gov/htdata/raster2/imagery/KingsBayGA_RGB_2009_1107/KingsBayGA_RGB_2009.vrt

or

gdalbuildvrt usgs.vrt /vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/USGS_Seamless_DEM_13.vrt

Oddly, the NOAA one gives some warnings about files that don’t, and have no reason to, exist. The USGS one doesn’t. The warnings don’t matter, but the difference in behavior is puzzling.

Once you have the local VRT, you can add that to a map. Unfortunately, my tests for the Kings Bay, GA, data were much slower for ArcPro than for QGIS, using the same local VRT that references the VRT on the web. QGIS (version 3.24.1) took about 1 minute to load and display. ArcPro (version 2.9.2) took 25 minutes. I don’t know why. It’s possible that ArcPro is ignoring the overviews embedded in the cloud optimized GeoTiffs and downloading all the data.

Conclusion

The primary gain you get from a VRT is an aggregation of potentially lots of files that can be referred to as a single file. It’s not the only way to do something like this. The group of files could be managed as a map or image service, though that requires a lot more management and infrastructure on the part of the data owner. Something like a SpatioTemporal Asset Catalog could also serve this purpose and provide additional flexibility, though they might lack the ability to encapsulate operations like the reprojection as that would be outside their scope.

From the perspective of managing the Digital Coast Data Access Viewer, we have found it very useful to have other people provide VRT files. For example, we can use the VRTs that were created for the emergency response imagery to allow people to do custom processing without having to either download or reference all the individual TIFF files. We simply reference the VRTs and let GDAL pull the data as needed for a job request. That provides greater data access options with less work for us.

Why these products?

For over twenty years, our lidar processing system has focused on data products suitable for scientists, engineers, and GIS professionals. People with access to software for further data manipulation and a need to be able to extract the data values. The main data types were point clouds, floating-point raster grids, and contours. None of these were particularly useful to someone with no special software and looking for a map of their area. With the requests we’ve had over time for more accessible products, this seemed an opportune time to try something.

How do I get them?

The DAV’s usual process for searching, adding to cart, and so on, hasn’t changed. The change is in the checkout selections. Note that this only applies to elevation datasets that start as point clouds (have three dots by their name) and not ones that start as raster DEMs (have four squares). We’ll probably add the DEMs at some point too. Previously you could choose an output product of points, raster, or contour, with the default being raster. We’ve added a selection for “Visual Representation”. You can thank our communications department for that name. I was going to call it “Map”.

Once you select Visual Representation, the file output selections will change to show the shaded relief TIFF, PDF, and KMZ options. You’ll also see options for the interpolation method and contour interval. Those options might be removed in the future since this is intended for an audience that doesn’t know much about interpolation methods. The defaults should be fine if you don’t want to muck with the options.

There is a caution that will pop up for the KMZ and PDF outputs if you select a large enough area that you might have trouble with the output. It uses an index of the number of pixels expected for a raster of the area divided by the contour interval and then takes the logarithm in base 2. That means that to lower the index by one, you’d need to double the grid cell size or double the contour or select half the area. The caution won’t stop you from doing anything, it’s just there to give you a heads-up. It’s only an estimate and there are important variables the user interface doesn’t know, such as the height range of your area.

That’s it for the interface changes, so let’s talk about what the products are.

Shaded Relief TIFF

The shaded relief TIFF is the simplest selection, at least right now. We generate a raster digital elevation model and apply a color map to it. We also generate a hillshade and then merge the two images. This gives a pretty image with elevation increasing from cold to hot colors that also has some shading to reveal detail. Most importantly, it’s an image that the photo editor, paint program, or whatever that came with your computer can open and display. You’ll also get a little image with the colorbar so you can tell what heights the colors represent.

Geospatial PDF

The PDF can be opened in a PDF reader, such as Adobe Acrobat Reader. The colored shaded relief image above is in the PDF, but it also has contours added on top of it. It is important to note that a large area with lots of contours is going to really bog down your computer. It works great for small areas, but a large area is painfully slow. How large is large? There are too many variables for a set answer, but a neighborhood is probably just fine while a county is almost certainly too big. In between depends on your computer, how hilly the area is (increasing the number of contours), and the grid spacing you pick. Increasing the grid spacing will let you work with a larger area by reducing the size of the raster and smoothing the contours.

For illustration, I’ll use the Adobe reader. Below is an example output for a neighborhood in Mount Pleasant, SC. You can see that we’ve add a legend as a colorbar in the lower left to provide the height meaning for the colors. This is the same colorbar that was a separate file for the tiff.

While we can’t see them in the neighborhood image, there are labels on the contour lines. They aren’t perfect but they are there if you zoom in. I’ll describe how those were made in the “Weeds” section at the end.

We can also turn on and off the different layers using the Layers menu on the left.

I should have the naming of the layers in better shape by the time I post this, but you can see that there is a layer for the color shaded relief raster, one for the labels, one for the legend and NOAA logo, and a weird one that should have been named contours. By clicking on the little eye in the box next to the names, you can turn the layers on and off.

KMZ (Keyhole Markup Language)

The KMZ format is generally used with Google Earth. This brings a lot to the table and is my favorite of the new products, but it will also be overwhelmed if you have too many contours. Below is the view in Google Earth after loading the KMZ version of the same area we saw with the PDF.

Google earth opens the view far enough away that the contours just make a big black box. However, it does let you see the imagery of the surrounding area and get oriented. We’ve also put the legend and the NOAA logo on the screen. We’re going to need to zoom in to see our data.

After zooming in, we can see the contours as individual lines. Google Earth is also providing labels for streets, businesses, and schools, that make it helpful to orient. However, there are no labels on the contours to show their heights. In the image above, I’ve clicked on a contour and you can see the table that pops up that includes the elevation. There are major and minor contour lines, with the major lines being thicker and are every 5 increments of your contour interval.

As with the PDF, you can turn the layers on and off. There are only two layers and the legend and logo go with the shaded relief image. Turning that shaded relief off does let you see the imagery with the contours overlain, which can be helpful to tell exactly where you are looking.

Another advantage of the KMZ is that you can add your own information to the view to augment it. It’s pretty easy to do and you’ll find lots of information on the web to add in Google Earth directly and to write your own KML files.

Next up is the details on creation. Or you can skip to the Conclusion where you’ll find some hints about your creation options.

Weeds – How we create the outputs

This section is for the folks that want to know some of the details of how the products were created. Most people probably don’t want to read this. You’ve been warned. All the products start with an elevation model (raster DEM) and contours created in the same way they would be if you requested them as data products from the DAV and there’s a previous post that describes how they were made. Open source software was leveraged to take those products and make the new ones, specifically GDAL and GMT.

Shaded relief image creation

As I write this, the shaded relief does not have the contours burned into it. It might in the future. The basic process steps are:

• Make a colored version of the DEM. This involves selecting a colormap to apply and using GDAL’s gdaldem color-relief.
• Create a hillshade from the DEM. This uses gdaldem hillshade.
• Darken the hillshade by applying a gamma of 2.0. You could do this with gdal_calc.py with a –calc argument like –calc=”uint8(((A/255.)**(2)) * 255)”, but I did it in python with the osgeo gdal module so I could use a lookup table and only do the calculation 256 times.
• Combine the colored DEM with the hillshade using gdal_calc.py and a –calc argument of “uint8((2*(A/255.)*(B/255.)*(A<128) + (1-2*(1-(A/255.))*(1-(B/255.)))*(A>=128))*255)” where A is the hillshade and B is the colored DEM.

There are many tutorials on the web for making a colored shaded relief map. Many of them also use the slope map derived from the DEM. That may be an enhancement for the future.

Colorbar or Legend

I’ve probably called it both a legend and a colorbar. A legend consisting of a scale with colors is usually called a colorbar. Without it, you don’t have any information about what the colors in the shaded relief image mean. To make it, we rely on the Generic Mapping Tools (GMT) suite of programs. We’ll need a couple pieces of information.

• The set of colors we want to use from low to high, often called a color ramp.
• The range of elevations those colors need to represent.

For our products, we’ve chosen two potential color ramps. One from Google called Turbo and the other from cpt-city called elevation. Turbo is used when the dataset has bathymetry points somewhere it it because it has blue at the bottom of the ramp. Elevation is used for datasets without bathymetry points and doesn’t have blue.

For the range of elevations, we can use our DEM file. You could get this a lot of different ways, including running gdalinfo -stats on the DEM and looking for the minimum and maximum values. I did it in python using NumPy and GDAL with this routine:

import numpy as np
import numpy.ma as ma
from osgeo import gdal
def get_dem_z_range(demfile):
    gdal.AllRegister()
    ds = gdal.Open(demfile, GA_ReadOnly)
    if ds is None:
        print("Could not open " + demfile)
        sys.exit(1)
    ndv=ds.GetRasterBand(1).GetNoDataValue()
    array=ds.GetRasterBand(1).ReadAsArray()
    masked_data = ma.masked_where(array == ndv, array)
    return np.min(masked_data),np.max(masked_data)

With those two pieces of information, we can make a color ramp with the elevation values we want with:

gmt makecpt -Cramp.cpt -Tminz/maxz > colorbar.cpt

Notice that arguments with GMT are a little different than most packages. There is no space between an option and its arguments. Now we want to make a graphic with a scale and numbers and all that. GMT’s outputs are often postscript files designed for use in publication. We’re going to make a graphic that is 0.5 cm in from the edge of the page in x and y, is 7.5 cm wide by 1.25 cm tall. It will be in portrait mode, so the width and height flip their usual sense. The labels for the scale will be at increments of 5. We’ll also tack on the units (in this case feet).

gmt psscale -Dx.5c/.5c+w7.5c/1.25c+e -Ccolorbar.cpt -Bx5+lElevation -By+lfeet -P > colorbar.ps

That gives us a pretty postscript file, but we need a raster graphic if we’re going to add it to some of these products. GMT has a psconvert command that leverages the ghostscript command. You can use it to convert to other outputs and specify how many dots/pixels per inch there are. For example, to make a PPM output, we can do:

gmt psconvert -A -Tm -E300 colorbar.ps

and it will produce colorbar.ppm rendered at 300dpi, cropping to only the content (the -A option). We render at 300 dpi for the PDF, but 72 dpi for the KMZ.

PDF creation

The first step for the PDF is to create the shaded relief map and colorbar above. To simply get the contours added, all you’d need to do is something like:

gdal_translate shaded_relief.tif outfile.pdf -co OGR_DATASOURCE=contours.shp -co OGR_DISPLAY_FIELD=Elevation -co LAYER_NAME='Color Shaded Relief'

Fancier contour lines

However, that’s going to put all the contour lines on without anything to differentiate them. We want to make the major contours, let’s say every 5th one, thicker and we’d like to label the elevation. That’s two different problems. While we can make an extra field in our contour shapefile to hold an OGR_STYLE that controls how thick the line is, GDAL/OGR’s pdf driver only understands the PEN option for lines and not the LABEL option. We do want to make that OGR_STYLE entry though so we can make thick major contours. You can see the result of that in the KMZ table for a contour line. You can see the options for styles in the OGR feature style specification. I added the column for OGR_STYLE with python and the osgeo module, but there are many ways you could do that.

Contour labels

The labels were a bigger problem and I can’t say my solution was elegant or anything like the results you get from something really designed to handle contours. While the GDAL pdf driver can’t put labels on lines, it can put labels on points. I made a new vector point layer, again with python and the osgeo module, by taking the middle vertex of each line and using that. I used the elevation for the line as the label for the point, again differentiating major and minor labels by making the major ones a little bigger with the OGR_STYLE field for those points. I also took the points on each side of the middle vertex to calculate an angle for my label at each point. The angle was atan2(dy,dx) where dx and dy are the difference in x and y for my two points on either side of my middle vertex and the atan2 output was converted to degrees.

Combining vector layers

The gdal_translate command expects only one data source and now we’ve got two. One for the contour lines and one for the labels as points. We need to combine them into one file. We can’t combine them into a shapefile, because that’s supposed to contain only one type of feature, not lines and points. I used ogrmerge.py to combine them into an sqlite file. That gives me a vector layer I can use in the gdal_translate command to output a pdf.

Colorbar and logo

We also want the colorbar and the logo and they need to sit in a fixed position on the page without having any georeferencing information. We do that with the EXTRA_IMAGES creation option for the pdf in gdal_translate. We can give that layer a name with the EXTRA_LAYER_NAME creation option.

KMZ creation

The KMZ needs a lot of the same things as the PDF. The contours and adding styles to the contours are the same, but we don’t need the labels because Google Earth will let us query the lines. In theory, the PDF lets us highlight lines by selecting them from the model tree, but all I’ve achieved by trying that is high CPU usage and the need to terminate Acrobat.

KMZ raster

For the raster, we’re going to make a KML super overlay with gdal2tiles.py. All we need to do is use gdalwarp to convert our shaded relief to web Mercator. For instance:

gdalwarp -t_srs EPSG:4326 shaderelief.tif webmercshade.tif

Then we can run gdal2tiles.py. I didn’t use many options, it was simply:

gdal2tiles.py webmercshade.tif -k -t 'Digital Coast is great' -x --processes=10 kml_out_dir

I could have skipped the gdalwarp step and used the -s option to set the input spatial reference frame, but that would have been another thing to track. It would probably be more efficient though. The gdal2tiles.py outputs everything into the kml_out_dir, including a kml file called doc.kml. Note that the KMZ is always going to be in Web Mercator, so the choice of projection during checkout doesn’t do much useful, however the reprojection to Web Mercator is not changing the vertical values, so the datum choices do matter.

The image KML is where we’re going to slip in the colorbar and the logo as images that stick in place in the window. We edit the doc.kml file using python and the xml.etree.ElementTree module. The idea is to read in the original doc.kml, add a ScreenOverlay element to the Document element with the appropriate positioning. We also need to copy our images to the kml_out_dir so they will be relative paths. The python to add an image looks like this:

import xml.etree.ElementTree as ET
tree = ET.parse("doc.kml")
root=tree.getroot()
doc=root[0]
scrover = ET.SubElement(doc,'ScreenOverlay')
name = ET.SubElement(scrover,'name')
name.text = "colorbar'
icon = ET.SubElement(scrover,'Icon')
href = ET.SubElement(icon,'href')
href.text = "colorbar.png"
overlay = ET.SubElement(scrover,'overlayXY')
# ovrx and ovry give the reference corner of the overlay. 0,0 is the bottom left corner, 1,1 is the upper right
overlay.attrib={"x":str(ovrx), "y":str(ovry), "xunits":"fraction", "yunits":"fraction" }
screen=ET.SubElement(scrover,'screenXY')
# scrx and scry give the screen location for the overlay's reference corner.  0,0 is the bottom left corner, 1,1 is the upper right
screen.attrib={"x":str(scrx), "y":str(scry), "xunits":"fraction", "yunits":"fraction" }
#write it back out again
with open ("doc.kml", "w") as file:
    file.write(ET.tostring(root, encoding='unicode')

To put the colorbar in the lower left corner, we use ovrx, ovry = 0,0 and srcx, srcy = 0,0. To put the logo in the upper right, we use ovrx, ovry = 0,1 and srcx, srcy = 0,1

Note that we’ve put the legend/colorbar and logo at the same level as the images. We could have put them in a different folder later to make them toggle separately.

KML contours

The KML contours are easily made with the ogr2ogr command, such as:

ogr2ogr -f KML contours.kml contours.shp

Combine KMLs

Now we have a problem. We have two KMLs. One with the imagery, colorbar, and logo, and the other with the contours. A KMZ file is just a KML and all its referenced pieces (e.g. images) zipped up. However, it can’t have more than one KML file in it. If it does, only the first one it finds is used and you can’t necessarily predict which one that will be. What we can do is combine the two KMLs into a single KML. Using the XML module again, my python routine looks like:

def kml_combine(rasterkml, vectorkml, outfile):
    """Combine the rasterkml and vectorkml files into a single kml"""
    rtree = ET.parse(rasterkml)
    rroot = rtree.getroot()
    vtree = ET.parse(vectorkml)
    vroot = vtree.getroot()
    rdoc = rroot[0]
    vdoc = vroot[0]
    rfold=ET.SubElement(vdoc,'Folder') # add a folder to the vector document to hold the raster contents
    for child in rdoc:
        rfold.append(child)
    with open(outfile, 'w') as file:
        file.write(ET.tostring(vroot,encoding='unicode'))

The vector (contour) KML already has a folder structure to it. The raster doc.kml doesn’t. We read in both files and add a folder to the vector as a sub-element of Document. We then put the raster contents in that sub-element and write out our new file.

KML to KMZ

All that’s left is to take our KML and make a KMZ. We have our output folder from the raster with the original doc.kml in it and it has all the images we need. We need to make sure our new combined KML is in there and we need to make sure the original raster and vector KMLs are not. Once we’ve done that, we can zip the entire directory contents into a zip file, giving it a kmz extension during creation or by renaming.

Conclusion

Thanks for sticking with me to the bitter end. Or maybe you skipped all the weeds. I hope you find the new products useful. I’d love to hear feedback on them. I’ll mention again that you want to keep the PDF and KMZ files relatively small. It’s not a problem to create a big file, but the usability is seriously impacted. You have some control when you create them. One of the most important is the grid size. Smaller grid size has more detail, but it’s that detail that is going to bog down your system. Larger grid size will make smoother contours and fewer vertices. Another is the contour interval. It defaults to the finest allowed interval based on the data accuracy. You may not need a contour every foot, especially if you have a lot of relief.

The products will initially show only for data sets that start as point clouds. Those are the ones with three dots by the dataset name and they won’t have DEM in the title. Since the new products all start from a raster DEM, there is no reason that we couldn’t make similar products for the datasets that start as rasters, though we will have to do something different to smooth the contours. If these turn out to be popular, I expect we will expand them to the raster DEMs.

What types of data are available?

The rockyweb site has a wealth of data available for the elevation enthusiast. This includes lidar points, original product resolution (OPR) DEMs, hydro breaklines, project and survey reports, USGS determined vertical accuracy values, and tile index and control point shapefiles. If you know the name of the data set you are looking for, you can probably just dive into the site. If you don’t know the name and just want to find out what data is available, a good place to start is the work unit extent spatial metadata (WESM) GeoPackage file.

Work Unit Extent Spatial Metadata (WESM) GeoPackage

In my frequent rockyweb stalking (er, researching), I have learned that new data is added to rockyweb almost daily. USGS is adding data to the site in work units (which are just subsets), not waiting until the entire data set has completed processing and validation, before making the data set available. To view the work units and all the data available, USGS provides the WESM. This is the spatial metadata that lists the data sets and their work units and is a good source of information for current lidar data availability and basic information about the lidar projects, including quality level, acquisition dates, and links to the project level metadata. This file is currently accessible from this rockyweb folder: https://rockyweb.usgs.gov/vdelivery/Datasets/Staged/Elevation/metadata. You’ll need software that understands an OGC GeoPackage, many of which are listed on the GeoPackage implementation page.

Lidar Points

The lidar points are the workhorse of the elevation products and the top need for many elevation data users. If you want to download the laz files in their delivered parameters, access is provided from this rockyweb folder: https://rockyweb.usgs.gov/vdelivery/Datasets/Staged/Elevation/LPC/Projects/.  The LPC in the folder name stands for lidar point cloud. As mentioned earlier, there is some variability in how the data sets are named. Some are listed with the state 2 letter abbreviation at the beginning of the folder name. Some are listed with ARRA or USGS_LPC at the beginning of the folder name, but the state two letter abbreviation follows. If you have downloaded the WESM.gpkg and found the name of the data set you are interested in, that will be a help in looking through the LPC/Projects folder to find your data.

Once you have found your desired data set, you can either download all the laz files in the data set using wget or uGet or, if you just need a single or few specific files, you will need to download the tile index.shp (discussed below) to determine which files you need to download. You might see a TL in the lidar point cloud data set folder. The TL is short for tribal lands and indicates that the files in that folder fall in a tribal land area.

If you want to download an entire data set, you’ll generally want some sort of download manager. An internet search will give you lots of options. Linux users probably already have a free command line tool called wget that works well.  Another free tool that might be more Windows friendly is uGet. The USGS has directions on how to download their data using uGet. 

Digital Elevation Models (DEMs)

If the lidar points are considered the workhorse product, the DEMs are the polished, easy to understand, regularly spaced grids that can be used for visualization and computing. When elevation data is contracted for collection, most of the time, DEMs are also delivered along with the lidar points. USGS calls these the Original Product Resolution or OPR DEMs and access to them is provided via this rockyweb folder: https://rockyweb.usgs.gov/vdelivery/Datasets/Staged/Elevation/OPR/. The organization is pretty much the same as the lidar points, with the state two letter abbreviation at the beginning of the folder name (with some exceptions), however, data that was added to rockyweb since December 2019 are located a folder down in the Projects folder. It is available here: https://rockyweb.usgs.gov/vdelivery/Datasets/Staged/Elevation/OPR/Projects/

If you are unsure of when the data was likely added to rockyweb, the best approach would be to just check both locations to determine if the OPR DEMs are available.

Hydro Breaklines

The hydro breaklines are also often a requested deliverable when lidar data is collected. These data are lines that show where water features or shorelines are located. These are used to help create the hydro-flattened DEMs. These datasets are usually provided in ESRI File GeoDatabase (gdb) format and can be accessed from the metadata/project/work unit/breaklines/ folder. 

Reports, Metadata, Tile Index

The last, but certainly not the least, important type of data available on rockyweb is all of the supporting information about the data. This accompanying information includes the metadata, reports, and shapefiles. As anyone who has worked with great data, the lack of any descriptive metadata about the data can limit its usability. Fortunately, USGS provides access to all this information so the user can make an informed decision about the suitability of the data for their project. All of the reference information is available in the metadata/project/work unit folder here:

https://rockyweb.usgs.gov/vdelivery/Datasets/Staged/Elevation/metadata/.

USGS has recently started adding their project (also known as work package, or WP) report which details their vertical accuracy value, as well as the dates of collection, work unit information, list of point classes, sensor used, las version, etc. This file can be identified by the WP_Report in the file name. 

Also, in the same location as the WP report is the vertical accuracy folder. In this folder you can find both the USGS and the contractor’s shapefiles for the control points and the checkpoints, as well as the USGS las checkpoint report which details the results of the Non-Vegetated Vertical Accuracy (NVA) and the Vegetated Vertical Accuracy (VVA) value determinations.

If the contractor collection, survey, or processing reports are needed, they can be found in the metadata/project/work unit/reports folder.

If the project, lidar point, DEM, or breakline metadata record from the original vendor is needed, those records can be found in the metadata/project/work unit/reports/vendor_provided_xml folder.

As mentioned earlier, if you want just a specific file or files, you will need to first download the USGS tile index shapefile. The LPC index is usually available as a shapefile and GeoPackage (gpkg) format.  These files are available in the metadata/project/work unit/spatial_metadata/USGS folder.

Cloud Access

As mentioned at the top, all of this data is also in the cloud and that is where USGS points you to for direct data access. That means you can do some data processing in the cloud without having to download the data. It’s in the AWS cloud, so you need to know the s3 bucket name to access it. The key is to recognize that web address https://prd-tnm.s3.amazonaws.com/index.html?prefix=StagedProducts/Elevation matches bucket location s3://prd-tnm/StagedProducts/Elevation. However, the point cloud data isn’t there. Instead it is in a requester pays bucket. You can find the addresses to use on the AWS open data registry page for 3DEP.

Summary

While the data available from the rockyweb site may be a bit of a challenge to navigate, the wealth of elevation information that USGS has made available for free download, from the lidar point data, to the DEMs, hydro breaklines, and the vital metadata, you too, may find yourself becoming a rockyweb creeper.

Contours

Perhaps the biggest change is with contours. DAV has been using a commercial tool that creates the contours from a triangulated irregular network (TIN) and then tries to simplify the lines and make them a little smoother. However, they were still pretty jagged and the simplify flag could result in output that was not topologically correct. There was the potential for contours to cross, which is generally considered “not good”. Although extracting the contours from a TIN will make them truer to the data, you don’t necessarily want every little dip and rise in someone’s lawn represented. Nor do you want the inherent noise in any point cloud to add irrelevant detail.

We’ve taken the approach that contours are meant for human interpretation and smoothness is desirable. Instead of a TIN, we’re going to generate a raster and then extract the contours from it. There are multiple methods of interpolation to create that raster, but we’ve noticed that the rasters from the ANUDEM program, using splines, are a bit smoother than some of the other methods we have (more on those methods below). We’re going to use that method for now, but it may change as we learn more. It’s possible that the interpolation method may become an option for contours, although it will probably need to be an advanced option.

To give you some idea of how different the contours are, below are two sets of contours over the same raster. The raster was generated as just the average value per cell with filling for small gaps. The first set of contours are from a TIN with simplification turned on. The second set is using splines to make a raster and contouring with gdal_contour.

Just the reduction of small polygons that probably have little meaning is enormous. The contours from the raster have roughly a third as many features as the ones directly from the TIN. It’s entirely possible that some of the features that were smoothed out of existence had significance for someone’s project. Unfortunately, there is no perfect answer for contours. They are already a simplification of the data, so if you need all the detail, you’re likely better off going to a raster or TIN version of the data.

Raster Interpolation

The DAV systems offers several options to generate a raster from a point cloud. Several of them used commercial tools requiring regular acquisition of a new license, which put the system at risk if we couldn’t get the paperwork done in time. These methods included TIN, average/mean, minimum, and maximum. Luckily, the open source PDAL and GDAL packages were able to step in for all of them. PDAL was already in use for inverse distance weighting (IDW).

Min, Max, and Mean

The page for writers.gdal in the PDAL documentation gives a good example of how PDAL and GDAL are combined to generate a raster from a lidar point cloud and that’s the pattern we followed. For the case where you want to fill small gaps, we use the GDAL command gdal_fillnodata.py with a maximum distance of 5 pixels. If you would like to fill further than that, just take the outputs and run gdal_fillnodata.py with a larger maximum distance.

TIN

The triangulated irregular network was a bit trickier. There is an example in the PDAL documentation, but it isn’t obvious. If you look at the documentation for the writers.raster, the first example shows how to do a TIN. However, it uses a filter called filters.faceraster that doesn’t currently show up in the list of filters on the PDAL page. That filter does have some needed options for setting the bounds of the output. However, we still have a problem with the TIN output, especially when we’re working at the coastline. We want to be able to limit how far it extrapolates so we aren’t making triangles out into a harbor if we don’t have bathy data. At the coast, we’re very likely to have a jagged edge of valid data and the seaward or lakeward areas should be empty instead of having a TIN interpolated value.

I don’t believe there is a simple automated solution to removing the areas you don’t want interpolated while retaining the ones you do using only the data itself. We do want to fill in areas like lakes and the holes left by buildings, but don’t want to extrapolate out to sea. Our previous solution would remove triangles that had a leg longer than a threshold. This kept it from going out to sea, but the threshold small enough to do that also left a lot of holes.

PDAL does not appear to have a feature to remove long triangles, so a different approach was needed. Our current solution is to take the point cloud and derive the concave boundary with PDAL’s hexbin filter. This boundary has polygons describing the areas with data, including having all the holes. A python script using GDAL’s OGR routines deletes interior holes in the hexbin output as we want to allow interpolation over the holes left by buildings. We use the boundary as a cutline for the raster in GDAL’s gdalwarp command, thus clipping out all the areas where we didn’t want interpolation.

The result is by no means perfect, but it’s pretty good for only using the data to drive it. The example above shows a large job in the Winyah Bay, SC area. Little lakes and building footprints have been filled in, but the bay has not, though it retains the islands. This is what we were looking for. However, you’ll note at the bottom that there are a couple of rivers that aren’t filled for a little bit because their gaps extend to the edge of the area I selected. That’s the trade-off. Gaps that extend to the edge will have their interpolated values clipped out.

I think the next step that would be required to improve the output is to use breaklines to determine the clipping or filling. While those exist for many datasets, the variations in their construction make using them in an automated system challenging. A lot of manual labor will need to go into cataloging and normalizing what we have.

Other Changes

Tile Sizes

Some of the changes above ended up enabling other beneficial changes. Previous products that used a TIN, including the contours, were using a 32-bit program and that had us limited to a tile layout that would keep the number of input points below about 15 million per tile. Using 64-bit PDAL, we aren’t hitting that limit and can make significantly bigger tiles. We haven’t tested the limits, but 120 million points is not a problem and we may go larger than that. That means a lot less tiles for you and fewer potential issues at tile boundaries. Tiles are now targeted to have 10,000 pixels per side.

DXF Format

We originally created DXF format from the shapefiles, created by the contour program, and then modified to add the z-component in each vertex. The converter was old, home-grown, circa 2000 software that never seemed to work well for AutoCAD users. While gdal_contour is still outputting shapefiles (maybe geopackage in the future), we’re now using ogr2ogr to create the DXF files and I’m hoping that will be a better product. Time, and users, will tell.

Future

Incorporating PDAL as more fundamental to the processing stream could allow some additional products to be generated. Already in the plans are potential products that can be used without a GIS system, though they will lack the ability to extract data values. PDAL can also output a mesh in PLY format and that could be translated to formats used by 3D printers, making it an interesting candidate. If you know of other outputs you’d like to see from the system, please let us know. Just recognize that it is automated, so anything that requires someone to put eyes on isn’t feasible.