2013-12-29

NullAttributeMapper Use Case: Translate Mesh Code to Geometry

(FME 2014 Beta build 14223)

There is a standard mesh system to divide all Japan area into uniform rectangular areas in geographic coordinate system.
mesh typemesh widthmesh height
Primary1 degree2/3 degrees (40 minutes)
Secondary7.5 minutes (Primary width / 8)5 minutes (Primary height / 8)
Tertiary45 seconds (Secondary width / 10)30 seconds (Secondary height / 10)

Every mesh has unique code to identify its location; the code is defined by the left-bottom corner coordinate (latitude, longitude) of the mesh.
mesh typecode formatdefinition of code elements
PrimaryAABBAA = North latitude * 1.5, BB = East longitude - 100
SecondaryAABB-CDC = Row index, D = Column index in the Primary Mesh
TertiaryAABB-CD-EFE = Row index, F = Column index in the Secondary Mesh
* Row / Column index are 0-based; start from south-east corner in the higher rank mesh.
* Hyphens may be omitted.

For example, the left-bottom coordinate of "5438-23-45" tertiary mesh is:
Lat. = 54 / 1.5 (deg.) + 2 * 5 (min.) + 4 * 30 (sec.) = 36.2 (deg.) = N 36d 12m 00s
Lon. = (38 + 100) (deg.) + 3 * 7.5 (min.) + 5 * 45 (sec.) = 138.4375 (deg.) = E 138d 26m 15s

The mesh system is a national standard of Japan. It's used frequently in various geographic data processing scenarios, I often translate mesh codes to geometries with FME workspace.
The following workflow example creates a rectangular polygon representing a mesh area based on the mesh code. Assume every input feature has a primary, secondary or tertiary mesh code as its attribute named "_mesh_code".

Bookmark 1
The three StringSearchers validate format of the mesh code and split it into code elements - AA BB [C D [E F]], the elements will be stored in a list attribute named "_mcode{}". Then, the AttributeCreator calculates mesh width and height depending on mesh type. Mesh type can be determined based on the number of code elements, i.e. a primary mesh has 2, a secondary mesh has 4 and a tertiary mesh has 6 elements.











Bookmark 2
The AttributeCreator calculates coordinate (latitude, longitude) of the mesh left-bottom corner based on the code elements. Here, the number of list elements from Bookmark 1 is variable - 2, 4, or 6. So I've used the NullAttributeMapper beforehand so that the number of elements would be 6 and also values of the ex-missing-elements would be 0. Thereby, the AttributeCreator can always calculate without any error.










The workflow above is a use case of the NullAttributeMapper. I think it can be also used in many cases other than handling <null>.
And I guess many countries have standards about geographic data. Hope the information will be exchanged through various channels.

====
2014-01-02: Number of StringSearchers (Bookmark 1) can be reduced to just only one with this regular expression.
-----
^([0-9]{2})([0-9]{2})(-?([0-7])([0-7])(-?([0-9])([0-9]))?)?$
-----
But in that case, since the number of list elements will be always 8, mesh type cannot be determined based on it. Other approach will be necessary instead of the ListElementCounter and the AttributeCreator. And, the NullAttributeMapper (Bookmark 2) settings should be:
- Map: Selected Attributes
- Selected Attributes: _mcode{3} _mcode{4} _mcode{6} _mcode{7}
- If Attribute Value Is: Empty
- Map To: New Value
- New Value: 0
=====

=====
Yes, I like scripting :-)
-----
# Python Example (FME 2014 Beta build 14223)
# 2013-12-30 Updated
import fmeobjects, re

def createMeshPolygon(feature):
    mcode = feature.getAttribute('_mesh_code')
    if mcode == None:
        return

    mcode = str(mcode)
    mtype, matched = None, None
    matchers = [(1, re.compile('^([0-9]{2})([0-9]{2})$')),
                (2, re.compile('^([0-9]{2})([0-9]{2})-?([0-7])([0-7])$')),
                (3, re.compile('^([0-9]{2})([0-9]{2})-?([0-7])([0-7])-?([0-9])([0-9])$'))]
    for t, m in matchers:
        matched = m.match(mcode)
        if matched:
            mtype = t
            break
    if matched:      
        codes = [int(n) for n in matched.groups()]
        for i in range(6 - len(codes)):
            codes.append(0)
         
        w = {1: 1.0, 2: 7.5 / 60.0, 3: 45.0 / 3600.0}
        h = {1: 40.0 / 60.0, 2: 5.0 / 60.0, 3: 30.0 / 3600.0}
        xmin = codes[1] + 100.0 + codes[3] * w[2] + codes[5] * w[3]
        ymin = codes[0] / 1.5 + codes[2] * h[2] + codes[4] * h[3]
        xmax, ymax = xmin + w[mtype], ymin + h[mtype]
     
        bndry = fmeobjects.FMELine([(xmin, ymin), (xmin, ymax), (xmax, ymax), (xmax, ymin)])
        feature.setGeometry(fmeobjects.FMEPolygon(bndry))
        feature.setAttribute('_mesh_type', mtype)
-----
=====
# Tcl Example (FME 2014 Beta build 14223)
# 2013-12-30 Updated
proc createMeshPolygon {} {
  array set formats [list  \
      1 {^([0-9]{2})([0-9]{2})$}  \
      2 {^([0-9]{2})([0-9]{2})-?([0-7])([0-7])$}  \
      3 {^([0-9]{2})([0-9]{2})-?([0-7])([0-7])-?([0-9])([0-9])$}]
  set mcode [FME_GetAttribute "_mesh_code"]
  set mtype 0
  foreach {t} {1 2 3} {
    if {[regexp $formats($t) $mcode m aa bb c d e f]} {
      set mtype $t
      break
    }
  }
  if {0 < $mtype} {
    if {[string length $e] < 1} {
      set e [set f 0]
      if {[string length $c] < 1} {
        set c [set d 0]
      }
    }
    array set w [list 1 {1.0} 2 [expr 7.5 / 60.0] 3 [expr 45 / 3600.0]]
    array set h [list 1 [expr 2 / 3.0] 2 [expr 5 / 60.0] 3 [expr 30 / 3600.0]]
    set xmin [expr $bb + 100.0 + $d * $w(2) + $f * $w(3)]
    set ymin [expr $aa / 1.5 + $c * $h(2) + $e * $h(3)]
    set xmax [expr $xmin + $w($mtype)]
    set ymax [expr $ymin + $h($mtype)]
    foreach {x y} [list $xmin $ymin $xmin $ymax $xmax $ymax $xmax $ymin] {
      FME_Coordinates addCoord $x $y
    }
    FME_Coordinates geomType fme_polygon
    return $mtype
  }
}
-----

No comments:

Post a Comment