Posted on by Joseph Ellsworth

DEM – Digital Elevation Model work with Scala

n51e001-london-ukDigital Elevation Models are the height maps of the Earth taken by Satellites.  This data was produced using the Shuttle Radar Topography Mission (SRTM)

The STRM3 data is available for most of the world while the SRTM1 data is available for most of the USA.   The SRTM1 data provides finer resolution than the SRTM3.

The basic quick start and topography docs are OK but working code is always easier to use.

I have included Sample Scala to download the entire world wide set of SRTM3 files along with the SRTM1 files for the USA.  There is also sample Scala code to read the .hgt files and convert them into grey scale renderings. 

Sample Images Rendered with Scala Code

Remember these are not actually maps.   They are DEM (Digital Elevation Models) visually expressed.    Darker areas indicate lower elevations.   Lighter areas indicate higher elevations.  Yellow dots are bad data.

Houston TX Area

n29w095-houston-tx

United Kingdom rendered using SRTM3 Elevation Data

n51e001-london-uk-france

Mauritius Island SRTM-3 Data

Mauritius Island SRTM3 DEM (Digital Elevation Map)

Los Angeles Valley all the way to Ocean

Los Angeles Valley DEM (Digital Elevation Model) using SRTM3 data from NASA

Mountains East of Los Angeles, CA USAn34w118-losangeles-ca

Puget Sound Near Seattle

Mountains Near Seattle WA USAn47w122-seattle-wa

Not sure where this was but was visually interestingn15w084

Rendering Notes

The NASA data is not ideal for rendering without a few tweaks.

  • Detect 0 foot altitude and render in blue to indicate sea.
  • Detect invalid data and attempt to compute based on surrounding land.  If that fails then render in Yellow.
  • Figure out how much altitude changes in this tile and use that to adjust the visual contrast.  This is the single most important aspect to obtain useful images.
  • Modify the default Java .jpg file save to save with higher quality.

Data Structure Notes

The hgt files are essentially binary files full of 16 bit integer in MSB format.   Java naturally reads integers in this format and scala uses the same libraries.   We have to read as shorts or Java tries to read 32 bit ints.

File names refer to the latitude and longitude of the lower left corner of the tile – e.g. N37W105 has its lower left corner at 37 degrees north latitude and 105 degrees west longitude. To be more exact, these coordinates refer to the geometric center of the lower left pixel, which in the case of SRTM3 data will be about 90 meters in extent.  You can use this information to compute which hgt tile a given lat, lon will be located in.

Disclaimer & Consulting

I don’t claim this is good Scala only that it works.    I am roughing together some ideas while testing certain Scala features.   I do offer consulting services but please remember I normally put a lot more effort into the code design than was invested here.
joeeatbayesanalytic phone number

 

Scala Source Code to Read and Render JPG from SRTM Data

import java.io.File
import java.io._
import javax.imageio._
import javax.imageio.stream._
import java.awt.image.BufferedImage
import java.util.Iterator
import scala.io.Source
 
// Default method for saving images does too much
// compression so this version saves the file with
// minimum compression 
// see also: ImageIO.write(out, "jpg", new File(destFiName))
def SaveHighQualityJpg( bimage : BufferedImage,  fiName : String) {
  val iter = ImageIO.getImageWritersByFormatName("jpeg");
  val writer = iter.next();
  val iwp = writer.getDefaultWriteParam();
  iwp.setCompressionMode(ImageWriteParam.MODE_EXPLICIT);
  iwp.setCompressionQuality(1);   // an integer between 0 and 1
  // 1 specifies minimum compression and maximum quality
  val file = new File(fiName);
  val output = new FileImageOutputStream(file);
  writer.setOutput(output);
  val image = new IIOImage(bimage, null, null);
  writer.write(null, image, iwp);
  writer.dispose();
}
 
 
class HGTElevation(val inFiName : String) {
  //TODO:  Move our Transform logic here
  val invalidElev = -32768
  val maxPossElev = 32768
 
  val xPix = 1201
  val yPix = 1201
  val elevations = Array.ofDim[Short](yPix, xPix)
  var maxElev = 0
  var minElev = maxElev
 
  // Control nature of greyscale 
  var elevOffset = 10
  val maxWhite = 230 // 210 //230
  val offset = ((255.0 - maxWhite) / 2.0).toInt
  var elevationStep = 1.0
 
  // Compute the pixel offset for a decimal decimal degree
  def computePixelOffset(decDegree : Float, numPixel : Integer) : Integer = {
    val fractPort = decDegree - decDegree.toInt
    if (decDegree > 0) {
        return (fractPort * numPixel.toFloat).toInt
    }
    else {
      return (numPixel.toFloat - (fractPort * numPixel.toFloat)).toInt
    }
  }
 
  /* Create a new elevation model as subset of the existing model
  to allow that area to be rendered in isolation */
  def subRegion(lat1 : Float, lon1 : Float, lat2 : Float, lon2 : Float) = {
    // TODO: Fill this in. 
    // Will ahve to make it possible to set xPix, YPix for this to work
    // correctly but
  }
 
  /* Save a retangular region specified between the 
   two geo points to ARC ASCII file format */
  def saveArcASCII(destFiName : String, lat1 : Float, lon1 : Float, lat2 : Float, lon2 : Float) = {
    // TODO: Fill this in
  }
 
  def readArcASCII(sourceFiName : String) = {
    // TODO; Fill this in 
  }
 
  //Adjust our colors ased on the elevations used in this
  //hgt file so we can clearly see the elevation changes then
  //write the .jpg file to disk.
  def makeJPG(destFiName : String) : BufferedImage = {
    val out = new BufferedImage(xPix, yPix, BufferedImage.TYPE_INT_RGB)  
    // Now Adjust our colors so we properly adjust them
    // to use the full color spectrum cased on the
    // maximum height.          
    println("Save " + destFiName);
    for (row < - 0 to yPix - 1) {
      for (col <- 0 to xPix -1 ) {
        val elev = elevations(row)(col);
 
        //println("row=" + row + " col=" + col + " elev=" + elev)
        if (elev == invalidElev) {
           // These error pixels mess up the map and show
           // as simple white so need to interpolate a reasonable
           // value based on surrounding land. 
           //
           // TODO: Spiral around these picking up non-error points in collection
           // then yeidl the average of the non error points.  Willing to spiral
           // out 8 or 10 pixesl to find a good location
           if ((row < 2) || (row >= yPix) || (col < 2) || (col >= xPix)) {
             out.setRGB(col, row, 0x0A3375) 
           }
           else {
             val estElv = (elevations(row-1)(col)  + elevations(row+1)(col) + elevations(row)(col -1) + elevations(row)(col+1)) / 4         
             if ((estElv > maxElev) || (estElv < 0)) {
              out.setRGB(col, row, 0xF4E842) // if still an error then color yellow
             }           
             else {
               val adjElev = offset + (estElv * elevationStep).toInt
               val gray = (adjElev << 16) + (adjElev << 8) + adjElev        
               out.setRGB(col, row,  gray)  
               println("invalid elev row=" + row, " col=" + col + " estElv=" + estElv)
             }           
           }
        }
        else if (elev == 0) {
          // unspecified as water normally means water so make blue
          out.setRGB(col, row, 0x0A3375) 
        }
        else {       
          val adjElev = offset + (elev * elevationStep).toInt
          val gray = (adjElev << 16) + (adjElev << 8) + adjElev
          out.setRGB(col, row,  gray)        
        }
      }  // for col
    } // for row  
    // save image to file
    //ImageIO.write(out, "jpg", new File(destFiName))
    SaveHighQualityJpg(out, destFiName)
    return out
  }
 
  def readHGT() {      
    val afile = new File(inFiName);
    val fin = new FileInputStream(afile);
    val bin = new BufferedInputStream(fin);
    val din = new DataInputStream(bin);
 
    // TODO: Compute minimum non-zero elevation and use 
    // it to compute the white scale for grids that only
    // contain higher altitude points.
    for (row <- 0 to yPix - 1) {
      for (col <- 0 to xPix -1 ) {
        // must read as short convert to int
        // val elev = din.readShort() & 0xFFFF; // if  need to convert to unsigned int
        val elev = din.readShort()
        elevations(row)(col) = elev.toShort
        if ((elev != 0) && (minElev > elev) && (elev != invalidElev)) {
          minElev = elev
        }
        if ((maxElev < elev) && (maxElev != invalidElev)) {
          maxElev = elev
        }
        // println("row=" + row + " col=" + col + " elev=" + elev)
      }
    }
    fin.close()
    println("fiName=" + inFiName)
    println("max Elevation=" + maxElev) 
    println("min Elevation=" + minElev)
 
    //elevationStep = maxWhite.toFloat / (maxElev - minElev)
    elevationStep = maxWhite.toFloat / maxElev
    elevOffset = (maxElev.toFloat * 0.8).toInt
    println("Elevation step = " + elevationStep + " ElevOffset=" + elevOffset)     
 
 
  } // func
} // class
 
// Load a set of hgt files and combine them into
// a single virtual image to allow cross border
// pictures..
def loadAndMerge(InputFileName : String) {
}
 
/* Load the specified file and convert it into a JPG file
with elevations represented as grayscale. */
def transformHGT(inputFiName : String,  outFiName : String) {
  val ahgt = new HGTElevation(inputFiName);
  ahgt.readHGT();
  ahgt.makeJPG(outFiName)
  println("\r\n")  
}
 
println("readHGTElevation.scala")
transformHGT("E:/srtm3/n-america/N15W084.hgt", "N15W084.jpg")
transformHGT("E:/srtm3/n-america/N29W095.hgt", "N29W095.Houston-tx.jpg")
transformHGT("E:/srtm3/n-america/N47W123.hgt", "N47W123.Seattle-wa.jpg")
transformHGT("E:/srtm3/n-america/N34W119.hgt", "N34W119.LosAngeles-ca.jpg")
transformHGT("E:/srtm3/africa/S21E057.hgt",    "S21E057.Mauritius-island.jpg")
transformHGT("E:/srtm3/eurasia/N51E001.hgt", "N51E001.London-uk.jpg")
transformHGT("E:/srtm3/africa/N04E029.hgt", "N04E029.jpg")

Scala Source Code to Download worldwide SRTM Data

/* fetchSTRMData 
  Use a Base URI to fetch a list of files from a NASA server and then
  download each of those files to a target drive on local system. 
  Unzip those files and delete the orignal archives.
 
  The NASA server is a little wierd and instead of returning a simle text file
  list of files returns formatted HTML so we have to parse the actual sub files
  from the original path and then combine them with the base path to 
  come up with the actual URI where the file is stored. 
 
*/
import scala.io.Source
import scala.util.matching.Regex
import scala.language.postfixOps
import java.io._
import java.net.URL
import sys.process._
import java.nio.file.Files
 
 
object fetchSRTMData {
  // Returns true if the path exists false otherwise.
  def pathExists(apath : String) : Boolean = {
    val fDesc = new File(apath)
    return fDesc.exists()
  }
 
  // delete the file if it exists
  def deleteIfExist(apath : String) {
   val fDesc = new File(apath)
   if (fDesc.exists()) {
     fDesc.delete()
   }
  }
 
  // Download a directory listing from remote server parse out the 
  // files named in href and return a new list of URI to download
  // Base Path:  http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Africa/ 
  // Line from HTML Returned:  <li><a href="N00E006.hgt.zip"> N00E006.hgt.zip</a></li>
  // Actual URI to fetch: http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Africa/N00E006.hgt.zip
  def getURIDownloadList(fileListURI : String) : Array[String] = {
    try {
      val html = Source.fromURL(fileListURI)("ISO-8859-1")
      val s = html.mkString
      val lines = s.split("\n")
      val pattern =  """.*href\=\"(.*)\"\>(.*).*""".r
      val amat = pattern.findAllIn(s).matchData 
      val fnames = for  {
          m < - amat
          grp = m.group(1)
          if (grp endsWith "zip") // other files like zip directory will not have .zip suffice
      } yield {
         println(" grp=" + grp)
         grp 
      }          
      return fnames.toArray
    } catch {
      case e: Exception => println("fail fetch URI: " + fileListURI + " e=" + e); return Array.empty[String];
    }
  }
 
 
 
  // Download File,   save in target path,  unzip and 
  // then remove original archive.  If the file has already
  // been downloaded then skip that file 
  def downloadFile(targetPrefix : String,  uriPrefix : String, fileName : String) {
    val uri = uriPrefix  + fileName
    val outName = targetPrefix + "/" + fileName
    val outhgtName = outName.replace(".zip", "")
    if (pathExists(outhgtName)) {
      println("STRM data already exists " + outhgtName + " skip")
      deleteIfExist(outName)
    }
    else {
      try {
        if (pathExists(outName)) {
          println("Zip file alread exists skip download " +  outName)
        }
        else {    
          println("Fetch " + uri)
          new URL(uri) #> new File(outName) !!
        }
        if (pathExists(outName)) {
          // Now Try to unzip the existing or new downloaded zip file 
          val cmd = "unzip -o  -d " + targetPrefix + " " + outName
          val result = cmd.!!
          println(" cmd=" + cmd + " result=" + result)
          if (result.contains("inflating")) {
            // if we get the inflating message we can assume the 
            // unzip worked OK so we can delte the zip file.
            deleteIfExist(outName)
          }
        }
      } catch {
        case e: Exception => println("fail fetch URI: " + uri + " e=" + e); 
      } 
    }
  }
 
 
  def main(args: Array[String]){
     System.out.printf("HTTPGet.scala")     
     if (args.length != 2) {
        println("fetchSRTMData.scala usage=   scala fetchSTRMData saveDrive ParentDirURI");
        println("""example scala fetchSTRMData.scala d:\datax\africa http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Africa/""");
     } else {
       val savePathPref = args(0);
       val uri = args(1);             
       println("savePathPref=" + savePathPref);
       println("uri=" + uri);
       val fileDesc = new File(savePathPref)
       if (fileDesc.exists == false) {    
          try {      
              fileDesc.mkdirs  
          } catch {
            case e: Exception => println("fail create directory " + savePathPref + " e=" + e); 
            return
         } 
       }
       val fnames = getURIDownloadList(uri)
       fnames.foreach( downloadFile(savePathPref,  uri, _))
     }
  }
}

Batch File to Call the Download Utility for Each Region

You may need to modify this for a hard disk where you have a lot of space available. For those running on linux you will need to convert the syntax to a compatible shell.

call scala fetchSRTMData.scala e:\srtm3\africa http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Africa/
call scala fetchSRTMData.scala e:\srtm3\austrailia http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Australia/
call scala fetchSRTMData.scala e:\srtm3\eurasia http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Eurasia/
call scala fetchSRTMData.scala e:\srtm3\islands http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Islands/
call scala fetchSRTMData.scala e:\srtm3\n-america http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/North_America/
call scala fetchSRTMData.scala e:\srtm3\s-america http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/South_America/
 
call scala fetchSRTMData.scala e:\srtm1\region_01 http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_01/
call scala fetchSRTMData.scala e:\srtm1\region_02 http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_02/
call scala fetchSRTMData.scala e:\srtm1\region_03 http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_03/
call scala fetchSRTMData.scala e:\srtm1\region_04 http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_04/
call scala fetchSRTMData.scala e:\srtm1\region_05 http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_05/
call scala fetchSRTMData.scala e:\srtm1\region_06 http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_06/
call scala fetchSRTMData.scala e:\srtm1\region_07 http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_07/

Reference Links:

Beginning Digital Elevation Model work with Python

4 Replies to “DEM – Digital Elevation Model work with Scala”

  2. I have the ArcAscii version working but it is too large to post in-line. ArcAscii is the the easiest form to use when downloading SRTM1 data from the USGS servers. The SRTM1 version 3 final provides much higher resolution data and renders nicer DEM than the hgt files we parsed here. If you send me a note then I can make it available on github.

    Factored out complexity to several classes
    Read ArcAscii
    Grab a rectangular sub region from a tile.
    Save as ArcAscii or .hgt
    Render sub regions of the hgt
    Better handeling of min / max elevation for shade coloring

  3. Pingback: Modeling Water flow using DEM (Digital Elevation Models) | Air Solar Water
  4. Pingback: Scala Resource Links Getting started | Predictive Engines, Search, Browse, Discovery High performance

Leave a Reply