Connecting...

W1siziisimnvbxbpbgvkx3rozw1lx2fzc2v0cy9zawduawz5lxrly2hub2xvz3kvanbnl2jhbm5lci1kzwzhdwx0lmpwzyjdxq

Creating maps from GRIB2 data in Scala (and Java) by Freddie Poser

W1siziisijiwmtkvmdqvmtavmtuvndevmjuvnduvymfsbc1zagfwzwqtymx1cnjlzc1iywnrz3jvdw5klwnsb3nllxvwltewotg1mtuuanbnil0swyjwiiwidgh1bwiilci5mdb4otawxhuwmdnlil1d

Scala lets you do so much, even what you might not first consider using it for!

Developer, Freddie Poser used Scala for creating maps of predicted temperature, have you ever used it for this function before? Check out how Freddie did it.

 

'Featuring GeoTools, Natural Earth Data, SBT Headaches and more.

In my last post I explained how to go about reading data from GRIB2 files using the netCDF library. Now we are going to go about creating a map of predicted temperature, complete with country borders.

This is going to be a high level overview of how I achieved this, the full source is available on github here.

 

Colour Interpolation

Ok lets start with a bit of code unrelated to grib2, weather or really maps, colour interpolation. We want to be able to define a colour scale and then smoothly interpolate so that we can get a colour for each data point. This isn’t too hard but if we want some nice, reusable code, we get something like this:

case class ColourScale (colours: (Double, Double, Double)*) {

  def mapValue(minVal:Double = 1, maxVal:Double = 1)(value: Double,interpolate: Boolean = true): (Int, Int, Int) = {

    val fi = ((value - minVal) / (maxVal - minVal)) * (colours.length - 1)
    val index = fi.toInt
    val step = fi - index

    if (!interpolate || step < 0.0001) {
      val (r,g,b) = colours(index)
      return (r.toInt, g.toInt, b.toInt)
    }

    val (r1, g1, b1) = colours(index)
    val (r2, g2, b2) = colours(index + 1)

    (
      (r1 + step * (r2 - r1)).toInt,
      (g1 + step * (g2 - g1)).toInt,
      (b1 + step * (b2 - b1)).toInt
    )

  }

  def reverse: ColourScale = ColourScale(colours.reverse: _*)

}

This class allows us to define a colour scale as a list of RGB tuples, in order from lowest to highest. The main function, 'mapValue' is curried, so that we only need to pass our min and max values once to create an anonymous function and then we can use that many times. This code is adapted from a stack overflow answer (which I can’t find now) but the way it works is to figure out the nearest colour in our scale, and how far we need to be from that colour. It then uses linear interpolation to create (acceptably) smooth transitions between colours. This isn’t the best way of doing this, but its certainly good enough for now.

An example colour scale might be as follows:

val CS_STANDARD = ColourScale(
    (0xad, 0x05, 0x05),
    (0xFF, 0x00, 0x00),
    (0xf7, 0x85, 0x38),
    (0xed, 0xf7, 0x37),
    (0x42, 0xe2, 0xf4),
    (0x00, 0x00, 0xFF)
  )

In my code I have defined a number of these on the companion object.

 

Making a Map

The next step is to actually plot every point in the grib file to an image. For this we are going to convert our grib2 dataset into an array and use that. In this code you can see a number of helper classes that I have created to make this process cleaner. The main one is the 'GhridSpec'. This gives us information about how the gridded data relates to real-world coordinates. This function is also designed to create images of any size, by scaling in either the X (longitude) or Y (latitude) direction.

def generateImage(data: UcarVariableGridWrapper, resolutionX: Option[Int] = None, resolutionY: Option[Int] = None, interpolate: Boolean = true, scaleRange: Option[(Double, Double)] = None): BufferedImage = {

    import data.gridFilter.gridSpec._
    
    val (width, height) =
      if (resolutionX.isDefined) {
        val w = resolutionX.get
        val h = (w.toDouble / longSize) * latSize
        (w.toInt,h.toInt)
      }
      else if (resolutionY.isDefined) {
        val h = resolutionY.get
        val w = (h.toDouble / latSize) * longSize
        (w.toInt,h.toInt)
      }
      else (longSize, latSize)
    
    val (minValue, maxValue): (Double, Double) = scaleRange.getOrElse((data.javaArray.map(_.min).min,data.javaArray.map(_.max).max))
    
    val xScale = width.toDouble / longSize
    val yScale = height.toDouble / latSize
    
    val colourMapper = (colourScale.mapValue _)(minValue, maxValue)
    
    println(s"xScale = $xScale")
    println(s"yScale = $yScale")
    
    val image = new BufferedImage(width, height, BufferedImage.TYPE_INT_RGB)
    
    for (x <- Range(0, width)) {
      for (y <- Range(0, height)) {
    
        val I_X = (x / xScale).floor.toInt
        val I_Y = (y / yScale).floor.toInt
    
        val d = colourMapper(
          data.getIndex(I_Y,I_X),
          interpolate
        )
    
        val (r,g,b) = d
        var p = (0 << 24) | (r << 16) | (g << 8) | b
        image.setRGB(x,y,p)
    
      }
    }
    
    return image
}

This code starts out by calculating how much it needs to scale each direction by inorder to match the desired size. Next it finds the maximum and minimum values of the data to use for the scale, unless we have defined an absolute scale. Next it iterates over every pixel in the image. 'I_X' and 'I_Y' represent the index in the grib2 data nearest to the pixel. This is a very dumb method with no interpolation in “blown up” images but is good enough. We then use our colour mapper function, defined from our colour scale to get the correct RGB values for the pixel. Finally we use some bit-shift magic to get the correct final value and set that in the image.

This function run on a temperature scale will look something like this. You can just about see the outline of major landmasses but nothing more.

 

Adding Borders

The final step is to add the outlines of landmasses and countries. To do this I first downloaded a map from the fantastic Natural Earth Data. We are interested in the shapefile data (the '.shp' and '.shx' files). To read this we will use GeoTools, a java library with a huge number of features but we will use 'gt-shapefile'. Currently an issue with SBT makes it quite hard to download this library. The work around is to turn of “gigahorse”, which seems to solve the issue. Our 'build.sbt' now looks like this:
updateOptions := updateOptions.value.withGigahorse(false)

resolvers += "Boundless" at "http://repo.boundlessgeo.com/main"

resolvers += "OSGEO" at "http://download.osgeo.org/webdav/geotools/"

// https://mvnrepository.com/artifact/org.geotools/gt-shapefile
libraryDependencies += "org.geotools" % "gt-shapefile" % "19.0"

// https://mvnrepository.com/artifact/org.apache.commons/commons-compress
libraryDependencies += "org.apache.commons" % "commons-compress" % "1.16.1"

// https://mvnrepository.com/artifact/commons-io/commons-io
libraryDependencies += "commons-io" % "commons-io" % "2.6"

Now we can open the map file, read it as a list of polygons, and extract their points. The code for doing this is as follows. The convertToIDXOption method simply converts the point from (lat,long) to an index in the grib array. If the point is outside of the boundaries of the array None is returned and the point is ignored. In future I want to have it be able to deal with this more gracefully but for now this is it.

class MapData(val data: List[List[(Double, Double)]]) {

  def interpret(filter: GridFilter): List[List[(Int, Int)]] = {

    data.map(
      _.map(
        P => filter.convertToIDXOption(P._1, P._2)
      ).filter(_.isDefined).map(_.get)
    )

  }

}

object MapData {

  def readFromShapefile(file: File): MapData = {

    import scala.collection.JavaConverters._

    val map: Map[String, String] = Map ("url" -> file.toURI.toString)
    val dataStore = DataStoreFinder.getDataStore(map.asJava)

    val typeName = dataStore.getTypeNames.head
    println(s"Reading $typeName")

    val iterator = dataStore.getFeatureSource(typeName).getFeatures().features()

    val polygons = ArrayBuffer[List[(Double,Double)]]()
    while (iterator.hasNext) {

      val feature = iterator.next()
      val geom = feature.getDefaultGeometryProperty
      val poly = geom.getValue.asInstanceOf[MultiPolygon]

      for (i <- Range(0, poly.getNumGeometries)) {
        val geom = poly.getGeometryN(i)
        val coords = geom.getCoordinates.map(T => (T.y, T.x)).toList
        polygons.append(coords)
      }

    }

    iterator.close()

    new MapData(polygons.toList)
  }

}

To add these points to the map we are going to amend our generate image method so that it takes in each polygon and draws lines between all the points. There are many ways of doing this, my implementation is below:


for (poly <- mapData.get.interpret(gridFilter)) {

  var lastPoint: (Int, Int) = null

  for (i <- poly) {

    val point = gridFilter.filterIDX(i)

    if (lastPoint != null) {

      val (startX, startY) = (lastPoint._1.toDouble, lastPoint._2.toDouble)
      val (endX, endY) = (point._1.toDouble, point._2.toDouble)

      val STEP = if (startX > endX) -0.05 else 0.05
      var xSTEP = STEP
      val GRAD =
        if (Math.abs(endX - startX)< 0.001) {
          xSTEP = 0
          if (endY > startY) 0.4 else -0.4
        } else (endY - startY) / (endX - startX)

      var posX = startX
      var posY = startY

      //println(s"From $startX,$startY to $endX,$endY with grad $GRAD and xSTEP $xSTEP")

      while (!(posX.round == endX.round && posY.round == endY.round)) {

        posX += xSTEP
        posY += GRAD * STEP

        image.setRGB((posY * xScale).round.toInt, (posX * yScale).round.toInt, 0)

      }

    }

    image.setRGB((point._2 * xScale).round.toInt, (point._1 * yScale).round.toInt, 0)

    lastPoint = point

  }

}

My line interpretation for loop is over-complicated and inefficient but it works, drawing black lines between the points. The lines are calculated at the scale of the grib data and then scaled like before with xScale and yScale.'

 

This article was written by Freddie Poser and originally posted on blog.vogonjeltz.com