sábado, 2 de febrero de 2013

Preparing a DEM and slope layer

In order to have more layers than can be used to create meaningful examples for the processes included in the OpenGeo Suite, I took the task of extending the Medford dataset that currently ships with the Suite and is used for most of the examples and use cases. In particular, I wanted to add a DEM and a slope layer, since the dataset has several vector layers, but lacks useful raster ones. I used quite a few SEXTANTE algorithms to do this, so I have decided to write down the process and publish it here as a case study.

The first thing I did was to look for good elevation data. The National Elevation Dataset was clearly my first choice. As it is divided in tiles, metadata shapefiles are provided, with polygons that represent the extent of each tile.

Data can be dowloaded from the National Map download website, selecting the region you want to download. The process is rather cumbersome (you have to give an email address and a download link is sent to you!!), but I ended up having two raster layers that covered the Medford area, as shown below.

These layers had two problems:

  • They covered an area that was too large for what I wanted (I was interested in a smaller region around the Medford city center, shown in brown in the image above)
  • They were in two different files. (The Medford limits fall into just one single raster layer, but, as I said, I wanted some extra area around it.

Both of them are easily solvable using SEXTANTE.

First, I create a rectangle defining the area that  want. To do it, I create a layer containing the bounding box of the layer with the limits of the Medford area, and then I buffer it, so as to have a raster layer that covers a bit more that the strictly necessary.

To calculate the bounding box , we can use the Polygon from layer extent algorithm

To buffer it, we have several alternatives. The GRASS v.buffer.distance algorithm has the option to leave square corners, which in this case it is interesting, since we are then going to use the resulting polygon to crop a raster layer, which has to be rectangular.

Here is the resulting bounding box obtained using the parameters shown above

With this layer that contains the bounding box of the raster layer I want to obtain, I can crop both of the raster layer that I downloaded, using the Clip Grid with Polygons algorithm.

Before clipping, the clipping vector layer and the raster layer to clip must be in the same CRS, since the algorithm assumes that. Originaly, they are not (the bounding box has the same CRS as the Medford limits layer, which is not the same as the downloaded DEM layers), although they appear to match, since QGIS is performing on-the-fly reprojecting before rendering, and also because both CRSs (EPSG 4326 and EPSG 4269) are rather "similar", so even without the on-the-fly reprojection everything would appear more or less "in place".

Reprojecting the clipping polygon can be done with the Reproject layer algorithm.

Now we can crop the raster layers.

Once the layers have been cropped, they can be merged using the Merge Raster Layers algorithm

A cellsize is needed for the merged layer. I used the same one of the input ones. You do not need to know it in advance before calling the algorithm. Just click on the button in the right-hand size of the text field and you will have a dialog to enter small mathematical formulas, and a list of frequently used values, among them the cellsizes and bounding coordinates of all available layers.

Note: You can save time merging first and then cropping, and you will avoid calling the Clip... algorith twice. However, if there are several layers to merge and they have a rather big size, you will end up with a large layer than it can later be difficult to process.

With that, we get the final DEM we want.

A slope layer can be computed with the Slope,Aspect,Curvature algorithm, but the DEM obtained in the last step is not suitable as input, since elevation values are in meters but cellsize is not expressed in meters (the layer uses a CRS with geographic coordinates). A new reprojection is needed. To reproject a raster layer, the GDAL Warp algorithm can be used. We reproject into a CRS with meters as units, so we can then correctly calculate the slope.

Here is the reprojected DEM.

With the new DEM, slope can now be computed.

The slope produced by the SAGA  Slope,Aspect,Curvature algorithm is expressed in radians, but degrees are a more practical and common unit. The Metric conversions algorithm will help us to do the conversion.

A similar conversion can be made using the raster calculator (whether the SAGA one or the GRASS one).
Converting to a different unit, such as percentage, can also be done with the calculator.

Reprojecting the converted slope layer back with the GDAL Warp algorithm, we get the final layer we wanted.

The reprojection processes have caused the final layer to contain data outside the bounding box we calculated in one of the first steps. This can be solved by clipping it again.

The process can be automated, as usual, using the modeler.

No hay comentarios:

Publicar un comentario en la entrada