Following on from Andy’s blog about how visualising photogrammetric data can be completed easily and freely; Andy has written another piece on how to take it to the next level. Over to Andy…
In my last blog I posted a short(ish) video showing some of the simple ways standard aerial photogrammetric data could be used. Well, what happens if you want to get your hands dirty and take it to the next level? If you asked that question the following blog shows just one of the many possible ways to make your data work even harder.
Using the various image datasets created in the previous blog can allow for some fairly complex analysis. But for large areas data management becomes a consideration and the facility to perform quick, simple and flexible analysis and visualisation becomes less easy. However, there are a myriad of ways and a number of software solutions to help this situation; I have picked just one for this blog to provide an illustration of what is possible.
Through a little bit of scripting I stored all the data sets in a single volume efficient file; using the las format. Those of you familiar with the LiDAR side of aerial data will already be used to las files, but they are far less common on the photogrammetric side of the discipline. The images below are taken from the same sample data as used in the previous blog, all stored in one single las file.
Left: RGB colourised point cloud created from the 1m 25km2 DSM and 12.5cm orthophoto.
Right: Same RGB point cloud blended with the DSM elevation.
Top: Oblique viewed, point cloud colourised by RGB.
Middle: Oblique viewed, point cloud coloured by elevation.
Bottom: Oblique viewed, point cloud RGB blended with elevation.
Top: Oblique viewed, triangulated point cloud colourised by RGB.
Middle: Oblique viewed, triangulated point cloud coloured by elevation.
Bottom: Oblique viewed, triangulated point cloud coloured by class.
So where to begin? Let me introduce the real triumvirate of geospatial data manipulation – GDAL/OGR, LibLAS and Python. GDAL/OGR and LibLAS are a set of Open Source Creative Commons licensed libraries that allow tasks to be carried out on geospatial data without needing to be able to code at the ‘uber’ level. GDAL for working on rasters, OGR for vectors and LibLAS for working with LiDAR data. Python is just one of many languages which can be used to bring it all together. The programming language you use is very much a case of ‘choose your poison’ but as all three can be installed from one package, OSGEO4W, it makes a good starter option. Links to the installation download and some setting up instructions are at the end of this blog.
Okay, so to address the elephant in the room, this does involve a little programming. But trust me; programming is not some big scary monster that only strange pasty creatures understand, speaking in arcane tongues while squinting at the sunlight. It is certainly easier to understand than a least squares bundle adjustment for aerial triangulation (and more fun he whispers!) – A contention vigorously disputed by our photogrammetric manager (Ed.). Also, the online community is huge, enthusiastically helpful and provides so much material that anyone can rapidly develop useful outputs.
So how do we put it all together? The LibLAS libraries provide access to the las format, an efficient and simple format for storing millions of points of data and associated georeferenced information. Let me say that again AND ASSOCIATED GEOREFERENCED INFORMATION! So all those other bits of data from our rasters: colour, NDVI and classification can be stored efficiently on a single XYZ point.
What, you do not believe it is that easy? “Proof!” I hear muttered from somewhere in the stalls. Okay, below are some example code snippets for using the libraries to carry out the storage of a raster DSM as points, attach the land-use classification and then apply an RGB colourisation to them. This formed the core of the script used to generate the example images earlier in this blog.
###Open required files – in this case the ESRI ASCII Grid file dsm1m.asc dsm = gdal.Open(r"C:\Users\An\Documents\Work\las\dsm1m.asc", GA_ReadOnly) if dsm is None: print "Fail" sys.exit(1)
###Create DSM raster transform dsmRows = dsm.RasterYSize dsmCols = dsm.RasterXSize dsmBands = dsm.RasterCount dsmTransform = dsm.GetGeoTransform() dsmX0 = dsmTransform dsmY0 = dsmTransform dsmPixelWidth = dsmTransform dsmPixelHeight = (dsmTransform * -1) ###Used a negative value as ASCII grid stored from lower left so need to count up the column
###Reads in dsm heights from Esri ASCII grid for i in range(dsmBands): dsmBand = dsm.GetRasterBand(i+1) dsmData = dsmBand.ReadAsArray(0, 0, dsmCols, dsmRows) dsmXOffset = dsmX0 / dsmPixelWidth dsmYOffset = dsmY0 / dsmPixelHeight
###Begins populating las file with points pt = point.Point() c = pt.color for k in range(dsmCols): for j in range(dsmRows): ###Adds a xyz value to a point pt.x = 1000 * (MinX + (j * dsmPixelWidth)) pt.y = 1000 * ((MaxY - MinY) + (MinY + (k * (dsmPixelHeight * -1)))) pt.z = 1000 * dsmData[k, j] ###add RGB values to colourise point cloud b = int(((pt.x / 1000) - rgbX0) / rgbPixelWidth) a = int((rgbY0 - (pt.y / 1000)) / rgbPixelHeight) c.red = rgbDataRed[a,b] c.green = rgbDataGreen[a,b] c.blue = rgbDataBlue[a,b] pt.color = c ###add classes -> class 20 and above to avoid reserved ASPRS class values f = int(((pt.x / 1000) - classesX0) / classesPixelWidth) g = int((classesY0 - (pt.y / 1000)) / classesPixelHeight) pt.classification = 20 + classesData[g,f] ###Writes the point to the las file opF.write(pt)
So there you have it, a next step in using standard photogrammetric data and if this has whet your appetites consider adding in R. A statistical library set – again Open Source and then welcome to the huge world of Geoanalysis, Geovisualisation and some really cool fun.
The OSGEO4W installers can be found at: http://trac.osgeo.org/osgeo4w/
Installation instructions for LibLAS libraries at: http://www.LibLAS.org/osgeo4w.html
This also installs Quantum GIS an Open Source GIS package to rival some of the commercially available packages.