DAV: New Products

The Digital Coast Data Access Viewer (DAV) has some new lidar product outputs that I think you might find useful. They are still experimental. That means there might be some glitches, but it also means this is a good opportunity for you to give feedback and make them more useful to you. The products are: 1) a shaded relief TIFF image; 2) a Geospatial PDF; and 3) and KMZ.

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. Recently there has been a push for greater equity in government services. Combined 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”.

An example of the checkout screen where the Visual Representation output product has been selected and the output format options populated.

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.

Example of shaded relief TIFF output.

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.

PDF output for a neighborhood. I used the average value to interpolation method generate the raster and that results in holes where the buildings and ponds are. The contours use splines that interpolate over those holes.

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.

Zooming in to the PDF lets you see the labels on the contours.

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

The different layers can be seen in the Layers menu.

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.

Opening view in Google Earth for our Mount Pleasant neighborhood.

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.

Google Earth after zooming in. The buildings and ponds that had no data in the PDF now show the underlying imagery. Clicking on a contour line shows a table with the line’s properties, including the vertical units in the line’s title.

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.

KMZ view with the shaded relief image turned off.

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.

Example colorbar using the Turbo ramp.

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.