Skip to content


30 Mar 17
Docs and Handoff — WA18

Prediction Run — Pass II

* rebuilt tiling
– there was a bug in the tiling code that caused gaps, due to a floating point truncation; using %f format for the float, as input to the transform, was insufficient resolution and the tiles were slightly malformed
– after retiling Inglewood and Humboldt Arcata, re-run training and prediction
– take advantage of the change, and clip Humboldt to a water layer

* Made a quick spatial join of Predicted set to Polys htarg2
– sent to garlynn for preliminary QC
– created a few QGis visualizations

Humboldt Data

Inglewood / Training
ctr_targ T/F; ctr_seg T/F; perc_overlap (0-1);

Humboldt / Prediction
coverage2 value (0-1);

MGRS Grid Reference

MGRS 10t 10s 11s start

buildings.buildings final

5563240 from all tables
5457201 removing duplicates


Google Earth Enterprise on Github

Developer See is using SEGNET based ML

24 Mar 17
Recognition Trials — WA18

Machine Learning — Tuning the Training Process

The BIS2 segmentation engine has three paramaters. For this series of tests evaltest, the imagery zoom level is always “full image resolution” and only one NAIP DOQQ is segmented (a portion of Inglewood). Iterating through input parameters, segmentation is performed and the result polygons are saved. This is repeated many times, and for several ‘modes’ of the imagery. Valid imagery ‘mode’ starts with true color ‘tc’, false color ‘fc’ and infrared ‘ir’; new products such as a Sobel Filter result can be added. (20mar17 – first sobel-enhanced segmentation added)

In evaltest, 675 segmentation run combinations of (t,s,c,mode) at full resolution on NAIP 2016 imagery were executed and the results are stored as shapefiles on disk. A small python program ( loads the shapefiles into a spatial table.

The 675 tables test tables:
 9 values for threshold (t) 10,20,30,40,50,60,70,80,90
 5 values for shape (a) 0.1,0.3,0.5,0.7,0.9
 5 values for compact (c) 0.1,0.3,0.5,0.7,0.9
 3 values for mode ir,fc,tc

    675 = 9 * 5 * 5 * 3

In the next step, statistics are calculated on the previous results (code in These stats are saved in tables in the ma_buildings relevance schema. Tables have parameter values embedded in their names in a structured way (see below).

Completion of this step is to choose the best parameter set of a search objective. In this experiment, we are matching shape-to-shape. Python generates SQL which contains a regular expression which filters the 675 tables of run result stats by encoded name, performs summary calculations on all rows in those tables (connected in the output by UNION ALL) and sorts all of that by a desired stats output column to discover outcomes in the test runs. Built-in Postgres math functions min() max() avg() stddev() are executed on three columns from each table, and counts of two more booleans. A SQL UNION ALL construct lists the chosen stats for all considered tables, one line per table, which then can be sorted in any way that is convenient to pick winners. Using this generalized format, statistics for any model run can quickly be compared. On the testing hardware, 450 tables of a few thousand rows each, are summarized and sorted for comparison in about 1.5 seconds.

Visual display of segmentation in a Training Area of Inglewood, California.

-- Five Spatial Tests are Performed on Every Segment Intersecting a Training Building --

CREATE TABLE relevance.evaltest_evaltest1_MODE_T_S_C
  gid integer PRIMARY KEY,
  class integer,      -- 0 not considered;  2 commercial bldg
  pctoverlap double precision,
  coverage1 double precision,
  coverage2 double precision,
  centr_seg boolean,
  centr_trg boolean

-- Structure of a model run stats table  23mar17
ma_buildings=# select * from "5bandt1_test2_5b_50_03_03" limit 5;
  gid  | class |    pctoverlap     |     coverage1     |    coverage2     | centr_seg | centr_trg 
 64311 |     2 | 0.547329611192823 | 0.970785850687081 | 35.9228431691447 | t         | f
 64559 |     2 | 0.137682284191657 | 0.836077333010945 | 0.54216208047878 | f         | f
 64107 |     2 |                 1 | 0.972281692288727 |   70.15447713594 | t         | f
 64698 |     2 | 0.728239585461471 | 0.924875366545267 | 17.6592649379458 | t         | f
 64721 |     2 | 0.923712040345278 | 0.943590470848252 | 30.8265007068801 | t         | f
(5 rows)

Relevance Query:
--    segmented poly (a)   < =  polys from a segmentation run w/ given params  
--    target poly (b)      <=  training polys, what we are matching  
--    raw areas (c)        <=  all known bldg polys in the area

        | insert into relevance.evaltest_evaltest1_fc_50_09_05                                                                                                                 
        |   select a.gid,                                                                                                                                                      
        |          1::integer,                                                                                                                                                 
        |          st_area(st_intersection(a.geom,c.geom))/st_area(a.geom),                                                                                                    
        |          st_area(c.geom) - st_area(st_intersection(a.geom,c.geom))                                                                                                   
        |            + abs(st_area(a.geom) - st_area(c.geom)),                                                                                                                 
        |          (st_area(c.geom) - st_area(st_intersection(a.geom,c.geom)))                                                                                                 
        |            / st_area(a.geom),                                                                                                                                        
        |          st_intersects(st_centroid(a.geom), c.geom),                                                                                                                 
        |          st_intersects(st_centroid(c.geom), a.geom)                                                                                                                  
        |     from seg_polys a                                                                                                                                                 
        |          join raw.areas c on st_intersects(a.geom, c.geom)                                                                                                           
        |           left outer join 
                       relevance.evaltest_evaltest1_fc_50_09_05 b on a.gid=b.gid                                                                                   
        |    where b.class is null  and  
               jobid='evaltest'  and  project='evaltest1'  and  
               mode='fc'  and  t=50.000000  and  s=9.000000  and  c=5.000000  
                 on conflict do nothing 

Relevance Tests

  pctoverlap => the area of the intersection between the target training poly and a segmented poly, divided by the total area of the segmented poly; meaning, the percentage of the segmented poly that is intersecting the target, converges to 1.0 on a segmented poly that is entirely intersecting the target

  cov1 => the area of a building that is not covered by the intersection with the segmentation poly plus the absolute value of the difference between the segmented poly and the building poly; meaning, a measure of non-fit; converges to zero on a perfect fit (actually, we changed this again, see code)

  cov2 => the area of the segmented poly outside of the intersection with the target, divided by the area of the segmented poly; meaning, the proportion of the poly that is not part of the target fit; converges to 0 on a perfect fit -m tc -f cov1 -o avg -g

  search all result tables in schema 'relevance' for test run summary tables matching
imagery mode 'tc' True Color with any parameters; from those tables, 
calc min() max() avg() and stddev() on all measures, but emit only the totals for the chosen column 'cov1', 
one line per table; sort those results on column 'average' with the default sort order  
of smallest first. In this example, segmentation parameter 20-07-07 gave the lowest
average value of cov1 for all test runs, with 20-07-03 closely following. 

  table                         class count minimum   maximum   average   stddev  centr_seg centr_trg
evaltest_evaltest1_tc_20_07_07   1   7237   0.000324  0.801199  0.052418  0.083101   3420   1414
evaltest_evaltest1_tc_20_07_03   1   6694   0.000476  0.800666  0.052979  0.083563   3158   1336
evaltest_evaltest1_tc_30_05_09   1   5404   0.000354  0.800875  0.053037  0.084356   2406   1232
evaltest_evaltest1_tc_30_03_09   1   6731   0.000324  0.800875  0.053139  0.082947   3202   1362
evaltest_evaltest1_tc_40_01_07   1   5310   0.000468  0.899413  0.053150  0.084180   2394   1223
evaltest_evaltest1_tc_30_05_03   1   5252   0.000320  0.801222  0.053332  0.082887   2362   1209
evaltest_evaltest1_tc_30_01_01   1   8033   0.000331  0.801431  0.053378  0.082763   3932   1425
evaltest_evaltest1_tc_30_01_05   1   8069   0.000331  0.800875  0.053403  0.081330   3936   1451
evaltest_evaltest1_tc_30_01_07   1   8212   0.000394  0.801083  0.053454  0.082504   4053   1459
evaltest_evaltest1_tc_20_05_03   1   9560   0.000320  0.801222  0.053536  0.081324   4769   1530
evaltest_evaltest1_tc_30_05_05   1   5379   0.000531  0.800875  0.053862  0.083052   2403   1210
evaltest_evaltest1_tc_20_05_05   1   9607   0.000476  0.801410  0.053940  0.083886   4795   1542
evaltest_evaltest1_tc_30_03_05   1   6658   0.000320  0.800875  0.054063  0.083892   3125   1324
-- Classified Segments with a Given Band and Segmentation --

ma_buildings=# select class, count(*) from evaltest_evaltest1_ir_50_01_03 group by class;
 class | count 
     1 |  4708
     2 |   324
(2 rows)

Topics in Training Target Selection

Process Steps

  • Generate 2000x2000 pixel tiles in the Training Area and the Search Area
    done -- 80 in Inglewood and 130 in Humboldt
  • Segment those tiles done -- chose 50_03_03
  • Load that segment data into the database done -- table sd_data.seg_polys_5b
  • Feed the segment data to ML tools (either train or predict) and
  • Write the ML results back to the database.

Building Size
  Very small buildings are more likely to be confused with other kinds of visible features, therefore remove training target buildings smaller than 90 m-sq

Next, what is a reasonable cutoff for buildings between 90 m-sq and MAX_SIZE
This is a stats distribution question.. We know the buildings of interest (commercial) look very much like other buildings. What is are identifying characteristics of commercial buildings that will help distinguish them from all buildings? Looking at distributions of features on all buildings is some curve, and the features of buildings we are looking for in the training set is a different curve.. Non-commercial buildings will have more occurrences in the low sq-m sizes, making false matches to commercial buildings more likely.

Let's start simple.

A quick look at the full training set shows that roughly 90,000 of the 250,000 bldgs are less than 150 m-sq, which means that 160,000 or 64% are >= 150 m-sq . My first reaction is ... look for >= 150 m-sq


20Mar17 Next Up

Segmentation is a general-purpose topic in imaging right now, so inevitably there will be other toolkits and advanced researchers that use it.. Some parts of the setup now would have to be generalized to handle other segmentation engines, but that is "out of scope" for this pass. For example, masking, poly splitting and poly merging, are well know in image segmentation. We have many options to explore, but must prioritize as we make the best use of the BIS2 tool we started with..

We have successfully added a fifth band to the raster layer 'mode'='sobel'.
Using this five-band input, combined with segmentation trials, definitely
improves the result. The toolchain did not have to change other than add the new layer.
Documenting the details of this is better left at the code level.

revised relevance query:

insert into relevance."inglewood_run1_5b_50_03_03" 
     -- pctoverlap                                                                                                                                                                                          
     -- coverage1                                                                                                                                                                                           
        (st_area(c.geom) + st_area(a.geom) -                                                                                                                                                                   
         2.0*st_area(st_intersection(a.geom,c.geom))) /                                                                                                                                                      
          (st_area(c.geom) + st_area(a.geom)),                                                                                                                                                                
     -- coverage2                                                                                                                                                                                           
        (st_area(c.geom) - st_area(st_intersection(a.geom,c.geom)))                                                                                                                                            
          / st_area(a.geom),                                                                                                                                                                                   
     -- centr_seg                                                                                                                                                                                           
        st_intersects(st_centroid(a.geom), c.geom),                                                                                                                                                            
     -- centr_trg                                                                                                                                                                                           
        st_intersects(st_centroid(c.geom), a.geom)                                                                                                                                                          
   from seg_polys_5b a                                 -- (3) segment polygons                                                                                                                            
    join la14.lariac4_buildings_2014 c on st_intersects(a.geom, c.geom) -- (4) building polygons                                                                                                            
      left outer join relevance."inglewood_run1_5b_50_03_03" b on a.gid=b.gid                                                                                                                                 
       where b.class is null  and  jobid='inglewood'  and  project='run1'  and  
         mode='5b'  and  t=50.000000  and  s=3.000000  and  c=3

 -- previous
 select a.gid,                                                                                       
    -- pctoverlap                                                                                                                                          
    -- coverage1                                                                                                                                           
    (st_area(c.geom) + st_area(a.geom) -                                                                                                                   
      2.0*st_area(st_intersection(a.geom,c.geom))) /                                                                                                      
        (st_area(c.geom) + st_area(a.geom)),                                                                                                                
    -- coverage2                                                                                                                                           
    (st_area(c.geom) - st_area(st_intersection(a.geom,c.geom)))                                                                                            
      / st_area(a.geom),                                                                                                                                   
    -- centr_seg                                                                                                                                           
      st_intersects(st_centroid(a.geom), c.geom),                                                                                                            
    -- centr_trg                                                                                                                                           
      st_intersects(st_centroid(c.geom), a.geom)                                                                                                             
 from seg_polys_5b a                                                                                                                                         
        join raw.areas c on st_intersects(a.geom, c.geom)                                                                                                      
          left outer join relevance."5bandt1_test2_5b_50_03_04" b on a.gid=b.gid                                                                                 
    b.class is null  and  jobid='5bandt1'  and  project='test2'  and  mode='5b'  and  
    t=50.000000  and  s=3.000000  and  c=4.000000  
      on conflict do nothing 

  from sklearn.neighbors import KNeighborsClassifier
  from sklearn.ensemble  import GradientBoostingClassifier

ma_buildings=# \d inglewood_run1_5b_50_03_03
Table "relevance.inglewood_run1_5b_50_03_03"
   Column   |       Type       | Modifiers 
 gid        | integer          | not null
 class      | integer          | 
 pctoverlap | double precision | 
 coverage1  | double precision | 
 coverage2  | double precision | 
 centr_seg  | boolean          | 
 centr_trg  | boolean          | 
    "inglewood_run1_5b_50_03_03_pkey" PRIMARY KEY, btree (gid)

We selected this parameter set to use for the test area.. It is good at selecting a lot of area from a large, single roof. Complex roofs that have colored parts, using the training tests we have, do not get joined using these settings.

This measure seems important, yet it is a substantial cost to calculate, at least the way we are doing it now.

 5415961.887 ms  statement: 
   update sd_data.seg_polys_5b a 
    set dist2road=(select st_distance(a.geom, b.geom) from roads b 
        order by b.geom < ->a.geom limit 1) 
     where dist2road is null and a.gid between 400001 and 500001

Completed Prediction Run

NAIP 2016 60cm segmented; predicted and classified based on Inglewood_5b model, sklearn Gradient Boost Tree (GBT).


Inglewood_5b training polygons, alongside classified segmented polygons

## Complete Machine-Learning Pipeline pass
##  23March17 -dbb

# tile the cousub for inglewood
./ -o inglewood -f 0603791400 -s 2000 -m 5b

# segment all the tiles in the inglewood directory
./ -m 5b -t 50 -s 0.3 -c 0.3 inglewood/

# load all the segment polygons into the database, compute dist2road
# with jobid=inglewood, project=run1
# segments go into sd_data.seg_polys_5b
./ -p run1 -m 5b inglewood inglewood/

# compare the segment polygons to the training polygons
# filtering any training polygons less than 150 sq-m in area
# save the results in relevance.inglewood_run1_5b_50_03_03
./ -b inglewood_run1_5b_50_03_03 -j inglewood -p run1
-m 5b -t 50 -s 0.3 -c 0.3 -z 150

# Using scikit learn and Gradient Boost Tree ML module
# train it using the segments in the database
# and save the trained object in inglewood_run1_5b_50_03_03-r0.5.pkl.gz
./ -j inglewood -p run1 -m 5b -r 0.5 -b
inglewood_run1_5b_50_03_03 inglewood_run1_5b_50_03_03-r0.5.pkl.gz

# tile an area for Humboldt based on a bounding box
# and place the tiles in directory humboldt
./ -o humboldt -b -124.17664,40.85411,-124.05923,41.00326 -s
2000 -m 5b

# segment all the tiles in the humboldt directory
./ -m 5b -t 50 -s 0.3 -c 0.3 humboldt/

# load all the segment polygons into the database, compute dist2road
# with jobid=humboldt, project=target
# segments go into sd_data.seg_polys_5b
./ -p target -m 5b humboldt humboldt/

# run the trained classifier inglewood_run1_5b_50_03_03-r0.5.pkl.gz
# against the humboldt segments and save the results to
# predict.humboldt_target_5b_50_03_03_r005
./ -b humboldt_target_5b_50_03_03_r005 -j humboldt -p
target -m 5b inglewood_run1_5b_50_03_03-r0.5.pkl.gz

Visualized output per parameter set:

Other Topics

ECN OSM Import





Jochen Topf reports that the multi-polygon fixing effort is going well and provides an issue thread to track the evolution.

On Mar 22, 2017 9:59 AM, "Brad Neuhauser" wrote:

> I think this is it?
> On Wed, Mar 22, 2017 at 11:54 AM, Rihards wrote:
>> On 2017.03.22. 18:37, Clifford Snow wrote:
>> > I am happy to announce that Microsoft has made available approximately
>> > 9.8 million building footprints including building heights in key
>> > metropolitan areas. These footprints are licensed ODbL to allow
>> > importing into OSM. These footprints where manually created using high
>> > resolution imagery. The data contains no funny field names such as
>> > tiger:cfcc or gnis:featureid or fcode=46003, just building height.
>> >
>> >
>> > Please remember to follow the import guidelines.
>> >
>> > The wiki [1] has more information on these footprints as well as links
>> > to download.
>> the link seems to be a copypaste mistake 🙂
>> > [1]
>> >
>> > Enjoy,
>> > Clifford Snow
>> >
>> >
>> > --
>> > @osm_seattle
>> >
>> > OpenStreetMap: Maps with a human touch
>> >
>> >


Development Seed


17 Mar 17
Recognition Trials -- WA18

object-based approaches to segmentation ... are becoming increasingly important as
   the spatial resolution of images gets smaller relative to the size of the objects of interest.

Andy Lyons, PhD
University of California, Division of Agriculture and Natural Resources
Informatics in GIS -LINK-

Topics for Detecting Buildings in Aerial Imagery:
- Object-based Image Analysis versus Pixel-based Image Analysis redux
- Change-detection using multi-temporal imagery
- Inherent Difficulty of Identifying Buildings in Imagery, even for a human
- LIDAR can increase detection dramatically, but is difficult to scale

Topics in Segmentation
- Object-based Image Analysis relies on Segmentation
- removing roads before segmentation
- removing vegetation before segmentation using IR layer
- using height information and shadow-detection

Measuring Success
- positive match rate to ground truth
- miss rate to ground truth
- false positives
- false negatives

Although many image segmentation methods have been developed,
there is a strong need for new and more sophisticated segmentation
methods to produce more accurate segmentation results for urban object
identification on a fine scale (Li et al., 2011).
Segmentation can be especially problematic in areas with low contrast
or where different appearance does not imply different meaning. In this
case the outcomes are represented as wrongly delineated image objects
(Kanjir et al., 2010).

excerpted from:

research with examples of Object-based Image Analysis (OBIA):

An object-oriented approach to urban land cover change detection. Doxani, G.; Siachalou, S.; Tsakiri-Strati, M.; Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci 2008, 37, 1655–1660. [Google Scholar]

Monitoring urban changes based on scale-space filtering and object-oriented classification. Doxani, G.; Karantzalos, K.; Strati, M.T.; Int. J. Appl. Earth Obs. Geoinf 2012, 15, 38–48. [Google Scholar]

Building Change Detection from Historical Aerial Photographs
Using Dense Image Matching and Object-Based Image Analysis -LINK-
Nebiker S., Lack N., Deuber M.; Institute of Geomatics Engineering, FHNW University of Applied Sciences and Arts Northwestern Switzerland, Gründenstrasse 40, 4132 Muttenz, CH

  more research -LINK-

A Few Definitions:
GSD -- ground sampling distance; real distance between pixel centers in imagery
DSM -- Digital Surface Model; similar to Digital Elevation Model
   eDSM -- extracted Digital Surface Model; post-process result, varying definitions
   nDSM -- normalized Digital Surface Model; varying definitions
FBM -- Feature-based mapping; store only features, not empty spaces
SGM -- Semi-global matching; computationally efficient pixel analysis
morphological filtering -- computationally efficient edge methods
automated shadow detection -- in OBIA
Hough Transform -LINK-
Representative Public Test Process:
  ISPRS Commission III, WG III/4
  ISPRS Test Project on Urban Classification and 3D Building Reconstruction


"... you dont want to be fooling with a lot of parameters.. you want paramater-free functions whenever you can.. thats one thing I think is exciting about neural networks (and the like), is that you just throw some images in there, and out come some magic answers. I think we will be seeing more of that."
scikit-image: Image Analysis in Python / Intermediate
SciPy 2016 Tutorial * Stefan van der Walt

Machine Learning Pipeline Components

Training buildings available in Inglewood, Calif.

Introduction: Source imagery is fed to a segmentation engine using parameter combinations of threshold (t), shape count (s) and compactness/smoothness (c). The results are saved. For every chosen entry in a table of known good training polygons, intersection is calculated and statistics are stored. Search is performed to find optimal "goodness of fit" by some criteria. The results parameters are then propagated to subsequent segmentation of new imagery, those results are stored, and searched for likely matches.

Known-good Building Polygons all have the following attributes available:

ma_buildings=# \d train_bldgs_08mar17
                 Table "sd_data.train_bldgs_08mar17"
          Column           |            Type             | Modifiers 
 gid                       | integer                     | not null
 building_type_id          | integer                     | 
 building_type_name        | text                        | 
 situscity                 | text                        | 
 shp_area_m                | integer                     | 
 geom                      | geometry(MultiPolygon,4326) | 
 naip_meta_1               | integer                     | 
 naip_meta_2               | text                        | 
 urban_ldc                 | integer                     | 
 compact_ldc               | integer                     | 
 standard_ldc              | integer                     | 
 intersection_density_sqmi | double precision            | 
 acres_grid_gf             | double precision            | 
 acres_grid_con            | double precision            | 
 use_du_dens               | double precision            | 

Theory of Determination of Relevance -- Clinton_Scarborough_2010_PERS

from Accuracy Assessment Measures for Object-based Image Segmentation Goodness: Clinton, et al

Pipeline Tools to date:  >options<
  create a table associated with subset of seg_polys and assign
  to a class based on mutual percent overlap with a truth polygon
    -i|--info  - report available jobid, project, t, s, c parameters and exit
    -b|--table tablename - table to create with relevance and class code (required)
    -j|--jobid name          -+
    -p|--project name         |
    -m|--mode ir|tc|fc|4b     |- filters for selecting specific records
    -t|--threshold n.n        |
    -s|--shape n.n            |
    -c|--compact n.n         -+

Segmentation Trials -LINEUP-
Usage: >options< trained-object.pkl.gz
    -i|--info  - report available jobid, project, t, s, c parameters
    -b|--table tablename - table with relevance and class code
    -j|--jobid name    -+
    -p|--project name   |
    -t|--threshold n.n  |- filter dataset to be extracted
    -s|--shape n.n      |
    -c|--compact n.n   -+
    -r|--relevance n.n - min relevance value to be part of class, default: 0.5
    -x|--test - Use 50% of data to train and 50% to predict and report
    -a|--algorithm - KNN - KNearestNeighbors|
                     GBT - Gradient Boost Tree, default: GBT

   conn = psycopg2.connect("dbname=ma_buildings")

predict-segments >options< trained-object.pkl.gz
    -i|--info  - report available jobid, project, t, s, c parameters
    -b|--table tablename - table to create with relevance and class code
    -j|--jobid name    -+
    -p|--project name   |
    -t|--threshold n.n  |- filter dataset to be extracted
    -s|--shape n.n      |
    -c|--compact n.n   -+
    -x|--test - do the prediction but don't write to the database
    scan intermediate relevance tables and compute descriptive statistics
  evaluate multiple t, s, c parameters of the segmentation process
  and try to pick the that gives up the best results

  for a given training set:
    * take 40% of the polygons and use them for training
    * take 40% of the polygons and use them for evaluation
    * save 20% of the polygons as a hold back set for final evaluation

  there are currently two proceses planned to do this:
  1. a short path:
     compute the relevance for all the polygons
     and try to increase the average mutual coverage
     and decrease the stddev on the average mutual coverage
     take the best 2-3 sets of parameters and run them through 2.

  2. a long path: (NOT IMPLEMENTED YET)
     compute the relevance for all the polygons
     create a trained classifier using the training polygons
     run the test polygons through the trained classifier
     evaluate the goodness of the fit

  In either case, save the parameters that gave the best fit.

  Relevance tables will be generated based on >jobid< >project< t s c

------------------------------------------------------- [-c] jobid dir
  -c    - optionally drop and tables
  jobid - string to uniquely identify this data
  dir   - segmentizer output directory to walk and load

For a given input file like project/tcase/test-fc.vrt
the segmentizer generates files in path project/tcase/:


We can generate files given this matrix:

  mode |  t  |  S  |  c  |  z  |
   fc  | 30  | 05  | 03  | 19  |
   ir  | 50  | ... | 07  | ... |
   tc  | ... | ... | ... | ... |

We can generate any number of values for t, s, c, and z.
Segmented results are stored in table seg_polys,

create table sd_data.seg_polys (
  gid serial not null primary key,
  project text,               -- from file path
  tcase text,                 -- from file path

segment-dir options dir
    Call the segmentation engine for each file in a directory, with the given settings
    Produce an HTML overview for human interpretation of results

      [-t|--threshold T]   - T1[,T2,...] thresholds, default: 60
      [-s|--shape S]       - S1[,S2,...] shape rate, default: 0.9
      [-c|--compact C]     - C1[,C2,...] compact/smoothness, default: 0.5
      [-v]                 - be verbose
      [--nostats]          - don't generate stats.csv file
      [--nolineup]         - don't generate the lineup image

The following database table enables reference of any dataset from a mapfile:

create table sd_data.seg_tileindex (
  gid serial not null primary key,
  project text,               -- from file path
  tcase text,                 -- from file path

Other Imaging Filters Dept.
 A Sobel Filter from sckit-image may be worth adding to the chain
The filter detects edges well (though they are not squared and various weights).
A target area could be passed through the filter and added as a segmentation clue.


08 Mar 17
Recognition Trials -- WA18

Project Big Picture
  "tripartite solution" -- Authoritative Sources of 2D building footprints, with varying attributes (dbname=auth_buildings), updates and ingestion are semi-manual; Openstreet Buildings with minimal attributes (dbname=osm_buildings), updates are automated; Machine-Learning generated building locations (dbname=ma_buildings), updates are machine-generated with supplementary attributes.

ML Theory
  Unsupervised Classification with pixel-based analysis (see below)
  Supervised Classification with Object-based analysis (what we are using, newer)

Project Details
  Database ma_buildings schemas now include: building_types, a classification system for buildings in all use_codes; US Census TIGER, road network, block groups, county subportions; Grid 150m Los Angeles County, vast, high-resolution landuse and urban metrics data; LARIAC4 Buildings 2014, every building in LA from SCAG; sd_data training_buildings, a table-driven source of ML training targets (see below), picked from 249,927 classified geometry records of commercial and MF5+ buildings in LA County.

Building a Training Set: filter a set of 2D building polygons of interest from LA County (table=sd_data.train_bldgs); run a segmentation on remote sensing data layers NAIP 2016, store the results (table=seg_polys).

Recognition: Segment a target search area using remote sensing layers NAIP 2016; use ML engine and training product to recognize likely matches; characterize the results.

Documentation & Handoff


* Training Set, continued:

- Garlynn's crosswalk JOIN table gives a LOT of category information on the buildings themselves. In effect, takes 300,000 or so non-residential buildings and puts a very fine grain classification on them.

"Each and Every Building in LA, provided by the County"
ma_buildings  la14.lariac_buildings_2014
code, bld_id, height, elev,
lariac_buildings_2014_area, source, date_, ain, shape_leng, status, old_bld_id,

"Building Categorization from Vision_California release"
  building_type_category text,
  building_type_id integer NOT NULL,
  building_type_name text,
  building_height_floors integer,
  total_far integer

12,Commercial,Store Combination,39,,1
13,Commercial,Department Stores,39,,1
15,Commercial,"Shopping Centers (Neighborhood, community)",41,,
16,Commercial,Shopping Centers (Regional),39,,1
17,Commercial,Office Buildings,32,,1
18,Commercial,Hotel & Motels,38,,1
19,Commercial,Professional Buildings,32,,1

First attempt at a JOIN with a VIEW.. did not work for QGis and had permissions problems within the database. Second attempt simply defines a TABLE in the sd_data schema. Success !

First pass at a finely categorized non-residential building set in LA County as ML Training Set candidates.

Misc URLs

.. and, Machine Learning problematic data

Characterize Training Set Buildings
- residential multifamily five units and up - are included
- how many buildings ? what sizes ? no filter on "small bldgs" yet

Attributes and Sources For Training Set (gw):
- Should record imagery attributes (direction of sun angle, year, etc). This will come from the NAIP metadata. WHY: It obviously matters from which direction the sun is shining, as shadows will be found on the opposite sides of objects from the sun. Time of year will tell about the presence or absence of leaves on deciduous vegetation. The year allows for time-series comparisons between different vintages of the same dataset.
- City (name) — this can come from a number of sources, including the LARIAC buildings shapes themselves (SitusCity). WHY: Different cities may have different attributes, citywide, in a meaningful manner.
- The area of the polygon, from LARIAC buildings training set (Shape_Area). WHY: We want to be able to filter by this attribute for only large buildings.
- urban_ldc from the UF base grid. WHY: These are super-urban areas, more than 150 walkable intersections/sqmi, most buildings here are skyscrapers.
- compact_ldc from the UF base grid. WHY: These are compact development types, more than 150 walkable intersections/sqmi, with a walkable street grid.
- standard_ldc from the UF base grid. WHY: These are suburban development types, less than 150 walkable intersections/sqmi,
- intersection_density_sqmi from the UF base grid. WHY: Intersection density is a critical variable in urbanism. It could be that the machine finds a strong correlation between it and other things. Let’s find out.
- acres_grid_gf from the UF base grid. WHY: Greenfield acres tell us how many non-urbanized acres are within the grid. If this is above some threshold, say 85% of the grid cell, we can probably skip that grid cell for ML analysis.
- acres_grid_con from the UF base grid. WHY: Constrained acres also tell us how many non-urbanized acres are within the grid. If this is above some threshold, say 85% of the grid cell, we can probably skip that grid cell for ML analysis.
- use_du_dens from the UF base grid. WHY: We may be able to do something clever, like filter out the residential records from the buildings dataset by only those grid cells from the UF base load with a DU use density above some threshold (I’m thinking 10-20 or so), which should get rid of some of the noise in the data.
- Building Size -- Filter all records by our threshold size: 2,000? 10,000? ...of 2D footprint. Maybe don’t try this during the first pass? Keep it in your back pocket for later? Not sure.


Server Layout

    bin        # all automation, plus small utils
        add-auth-building-layer  do-backup            json-to-postgres    osm_building_stats
        create-auth-buildings     load-osm-buildings  reload-osm-buildings
        createdb-dee1            fetch-osm-ca-latest  load-roads          wget-roads
    bis2       # the Image Segmentation engine
    cgi-bin    # mapserver(tm) delivers the view of data
    data       # links to large data stores, raw inputs, csv files
        auth_buildings  dbdump  misc  naip  sd_data  tiger2016-roads
    html       # mapserver site 
    maps       # mapserver layer definitions
    src        # misc code
        code_dbb  ECN_osm_import  naip_fetch  naip-fetch2  naip-process  osm-buildings  sql


California Department of Conservation
new maps site:
"... For further information or suggestions regarding the data on this site, please contact the DOGGR Technical Services Unit at 801 K Street, MS 20-20, Sacramento, CA 95814 or email "

also: < - referred from a wikipedia page, no other info  

Some Breadth on (pre)Machine Learning Topics
(see gw/rdir/snippets_misc for more like this)
Pixel Matching -- Notes on Maximum Likelyhood Classification Analysis
"... is a statistical decision criterion to assist in the classification of overlapping signatures; pixels are assigned to the class of highest probability." -LINK-
"one approach is to use Bayes' Rule and maximize the agreement between a model and the observed data"

  Chapter 5 - Model Estimation, page 116
The Nature of Mathematical Modeling; Gershenfeld, Neil, Cambridge University Press 1999

Maximum Likelyhood Estimation -WIKIPEDIA-
Bayes' Rule -WIKIPEDIA-


01 Mar 17
Recognition Trials -- WA18

* prep for segmentation

$ -o  /sand480/extractions_work/t2.tif -z 2 40.87427 -124.08024

  - web visible results -LINK-



from Multiresolution Segmentation: an optimization approach for high quality multi-scale image segmentation -- Martin BAATZ und Arno SCHÄPE -LINK-
Design Goals: (for a generalized analysis method ed.)

Multiresolution: Objects of interest typically appear on different scales in an image simultaneously. The extraction of meaningful image objects needs to take into account the scale of the problem to be solved. Therefore the scale of resulting image objects should be free adaptable to fit to the scale of task.

Similar Resolution: Almost all attributes of image objects – color, texture or form -
are more or less scale-dependent. Only structures of similar size respectively scale
are of comparable quality. Hence the size of all resulting image objects should be of comparable scale.

Reproducibility: Segmentation results should be reproducible

Universality: Applicable to arbitrary types of data and problems

Speed: Reasonable performance even on large image data sets

commentary: Let us focus here on the specific segmentation/extraction parameters (the "scale of the problem" is ambiguous in English, as it could refer to the total ambitions of the result set size). That is, the content of a scene by number of objects to be recognized, and the information in the scene, by measurable detail in the recorded bits. Both number of targets in a scene, and the detail available on each target, vary directly by scale -- aka zoom level.

The NAIP CA 2016 imagery has both IR and RGB layers, and with both layers, any reasonable desired "zoom level" is available, greater and lesser than the native resolution. Some questions here include:

* what combination of zoom levels + segmenting + analysis, within reasonable compute bounds given existing software and hardware, will yield accurate and consistent recognition for the size of the object we are searching for ?

* what segmentation detail will be most useful to identify likely targets?

* can we combine segmentation results from RGB and IR/Gray, in some way with our setup, to increase recognition? Do we need two training sets for ML, or can the two layers be combined ?

* if this search is focused on "big buildings", can we disregard completely "small buildings" footprints, like a typical residential, single-family home for example? Or, do we need the similarly shaped, but smaller footprint building in the result. Generally speaking, one building is more similar looking to other buildings, than to other kinds of features in the scene... except maybe parking lots !

* Can negative examples be used in a training set ? In other words, a segmented parking lot is "not a building".

We can think of "how big is a typical target, at a given scale" and "how many segments typically compose one target". Some explanation of these topics follows.

Consider a test area that is only 20% larger than a typical building of a certain class of interest. There is a huge amount of detail to match on, but, how many image analysis operations would it take to accomplish a survey ?

A counter-example is, analysis of an image at a zoom level such that there are one-hundred targets of interest. Each target is no longer a detailed polygon, but instead could be identified by the pattern of segments including the surrounding group, e.g. pixels for the building, its shadow, the road leading to it, cleared land around it, etc.

Can information from analysis at different scales, be combined in some rigorous way ? Does this overlap with combining non-segmentation data with segmentation data ? for example the presence of OSM polygons in a search area..


In case of boredom, a boot disk dropping off of its RAID is a remedy...

$ cat /proc/mdstat
Personalities : [raid1] [linear] [multipath] [raid0] [raid6] [raid5] 
[raid4] [raid10]
md0 : active raid1 sdh[2] sda1[0]
244007744 blocks super 1.2 [2/1] [U_]
[====>................] recovery = 23.5% (57386432/244007744) 
finish=53.0min speed=58672K/sec

unused devices: 

Mask Files ?
- it seems there may be something missing in the extraction chain. There is an option in the gdal_translate code to "generate internal masks" but the way the processing chain is right now, the dot-msk files were left as a side product. This may need to change.

Map URLs as the Basis of Linking
The short-form of the map urls has zoom information and a center LatLon. This is useful at all steps in the analysis, from human to machine internals. Adapt the tool to use those zoom levels when reading NAIP imagery.

Search Notes:
* we think that "distance to a road" should be included in a preliminary search as a weight,
and that making effort eliminate whole areas from the search is worthwhile.

24 Feb 17
Recognition Trials -- California AB 802 Support

James S. writes:

There is probably no documentation you haven't seen. The Python API source code is documentation
 of sorts. There is also a brief FAQ on

Honestly, different customers will either use or not use different parts for their specific
workflow. We lay it all out and let the integration happen at the customer end. As such there 
are different pieces and you are welcome to use, disregard, or innovate on top of them.

You're certainly putting more software engineering effort into it than most.

The is optional for doing the shape comparison to determine the "goodness" of the
segmentation match to a ground truth shapefie. That embodies the gist of the Clinton paper.
Another potential post processing is supervised classification on the image stats.

Note and modules that wrap the suggested workflow.

You are certainly welcome to build your own front end in Jupyter if that works for you. That
sounds like a good iterative/exploratory take on it. Some other customer in Africa just wanted 
a GUI with four boxes. The land cover mapping customer in Mexico just did the whole country 
at 10,0.5,0.5 and moved on. As I say, fit for purpose.


  • Threshold | 1-255 ? | close to 0 is "many colors in final"; 255 is "few colors in final"
  • Shape Rate | [0.0,1.0] | 0 -> "many shapes"; "few shapes" -> 1.0
  • Compact/Smoothness | [0.0,1.0] | 0 -> "compact; "spread" -> 1.0

Base Example -IMG- Arcata_Ex0.png
  A harsh first test; much detail, many similar objects.. ..
  run in grayscale (IR) only
    roughly 70 variations are produced (see exArcata dir)

Example Vectorization using BIS:
  Here is an IR layer image of Arcata, close up and shown at 50% intensity. Using the triple-inputs to BIS (threshhold, shape rate, smoothness) get a vector layer back. Do this twice, varying only the first parameter (threshhold). Notice how the same edges are picked in both cases, only more of them with a lower threshold.


  • ArcataEx0
  • ArcataEx0 50_04_04
  • ArcataEx0 30_04_04

Some observations: these vector result edges do follow the edge of change in pixel values, exactly. So the odd shapes of shadow and lighting, do carry through to the vector results. There appears to be no way for the vectorizing engine to estimate and resolve implied straight edges through shadows, for example. note: internally, BIS calls gdal.Polygonize() -LINK- to perform the vectorization.

-- Some buildings in Chico, at zoom level 14 in Infrared (IR)
-- RGB Color, same scene in 2014
-- osm outlines, more to do !

Another form of output from the BIS engine is statistics available in a dot-csv file, each and every geometry is summed, including area in square units, mean luminance per pixel, and other statistical measures.. These statistics can be used to calibrate engine response in a feedback loop, using larger frameworks.

Questioning Machine Learning:
  What is building and what is not building, in this image.

One question a ML process can answer is, is this (noisy thing) similar to this other noisy thing. Not the same to make a clean line in vector geometry! We may be asking several questions at once.. Where on the image is "building pixels" and where are not, is different than, provide a simple geometry representing each building. What is building, what is not building, in an image... is a different question than, show me a good representation of a building.

* developing a training set for supervised ML is work - a training set shows the learning system what a "right answer" is.. more than one hundred thousand, well-chosen examples is not uncommon in a training set. Developing the training set is not to be taken for granted.

* we have targets that do not change in life, but the view of them changes, due to weather, lighting, focus, image processing.

Example BIS Automation Run from the Command Line using IPython

$ cd /home/dbb/CEC_i7d/bis-2/bisapi-hamlin/bis-workd
$ ipython
Python 2.7.6 (default, Oct 26 2016, 20:30:19) 
Type "copyright", "credits" or "license" for more information.

IPython 5.2.2 -- An enhanced Interactive Python.
In [3]: from automate import automate
In [4]: 
    ...: files = automate('/wd4m/exArcata/Arcata_Ex0.png', t=[42,52,62,72,82], 
    ...:   s=[0.5], c=[0.3,0.7],do_stats=True, do_colorize=False,
    ...:   do_train=False, do_classify=False, do_outline=True, do_lineup=True)

## -- BIS virtualenv Fresh Install --
ipython / jupyter / scipy / pyside  virtualenv re-install

    source /home/dbb/CEC_i7d/bis-2/bisenv/bin/activate bisenv
    pip install ipykernel
    python  -m ipykernel install --user --name bisenv --display-name "bisenv"
    pip install shapely fiona
    pip install scikit-image
    pip install numpy rasterio
    pip install pyproj cartopy
    pip install psycopg2
    pip install gdal   < - gives 2.1.3 which appears to work ok with 2.2dev

## - BIS working set of packages in virtualenv-
(bisenv) dbb@i7d:~/CEC_i7d/bis-2/bisapi-hamlin/bis-workd$ pip freeze
affine==2.0.0.post1  Cartopy==0.14.2   Fiona==1.7.3    GDAL==2.1.3
ipykernel==4.5.2    ipython==5.2.2    ipython-genutils==0.1.0
jupyter-client==4.4.0    jupyter-core==4.3.0
matplotlib==2.0.0    numpy==1.12.0    Pillow==4.0.0    psycopg2==2.6.2    ...
pyproj==    pyshp==1.2.10    PySide==1.2.4    ...
rasterio==0.36.0    scipy==0.18.1    Shapely==1.5.17    ...
six==1.10.0    subprocess32==3.2.7    tornado==4.4.2    ...

## Fiona CLI tool - (no PG driver)
$fio env
  ESRI Shapefile (modes 'r', 'a', 'w')
  GPKG (modes 'r', 'w')
  GPSTrackMaker (modes 'r', 'a', 'w')
  GPX (modes 'r', 'a', 'w')
  GeoJSON (modes 'r', 'w')
  OpenFileGDB (modes 'r')

i7d Disk Updates:
 * naip-fetch2/
   - add dir /sand480/tmp


Google Search \ filetype:pdf inurl:asprs object

BIS Test 2

visually picked an area in IR; important variables: zoom level, color/gray; composition

/osmb/?zoom=17&lat=40.8835&lon=-124.08918&layers=0000BTFFFFFF -LINK-

Use the DOQQ labels layer in the map to find the filename of the right image on disk; picked out a subset using rio tool

rio info final-ir/m_4012408_sw_10_1_20140607.tif

take one-hundredth of a degree in each direction, from some buildings

$ rio clip input.tif output.tif --bounds xmin ymin xmax ymax
--bounds #(-124.087-0.01) (40.878-0.01) (-124.087+0.01) (40.878+0.01)
rio clip final-ir/m_4012408_sw_10_1_20140607.tif /tmp/work1.tif --bounds -124.097 40.868 -124.077 40.888

try expanding the file with gdal to 2x resolution

gdal_translate -outsize 200% 0 /tmp/work1.tif /tmp/work2.tif
-rw-rw-r-- 1 dbb dbb 24M Feb 20 19:09 work2.tif
-rw-rw-r-- 1 dbb dbb 6.0M Feb 20 19:00 /tmp/work1.tif

but, the image result looks exactly the same, just 4x the size on disk.
onward -- go to BIS work dir

dbb@i7d:~/CEC_i7d/bis-2/bisapi-hamlin/bis-workd$ source /home/dbb/CEC_i7d/bis-2/bisenv/bin/activate bisenv
(bisenv) dbb@i7d:~/CEC_i7d/bis-2/bisapi-hamlin/bis-workd$ ipython

In [1]: from automate import automate

files = automate('/tmp/work1.tif', t=[42,48,54,60,72], s=[0.1,0.2,0.4,0.6,0.9], c=[0.5],do_stats=True, do_colorize=False, do_train=False, do_classify=False, do_outline=True, do_lineup=True)

CRASH in outline burn()
repeat with do_outline=False
also note, turn off stats if they are not being used

files = automate('/tmp/work1.tif', t=[42,48,54,60,72], s=[0.1,0.2,0.4,0.6,0.9], c=[0.5],do_stats=False, do_colorize=False, do_train=False, do_classify=False, do_outline=False, do_lineup=True)

automate() invokes the engine five times, each time thresholds are evaluated from one to the largest threshhold in t, and len(t) number of outputs are emitted. A run in this example allocates 12GB of RAM.... about 40 seconds per run.. 5 runs with 5 outputs each, in this example.

example output:



New Data - NAIP CA 2016 0.6m
copying ...

$ /naip_16/d2/naip/ca_60cm_2016/40124$ ls -lh m_4012408_sw_10_h_20160528.tif
-rwxrwxrwx 2 root root 439M Sep 29 04:22 m_4012408_sw_10_h_20160528.tif

##-- test area in Arcata, CA
rio info /naip_16/d2/naip/ca_60cm_2016/40124/m_4012408_sw_10_h_20160528.tif
{"count": 4, "crs": "EPSG:26910", "colorinterp": ["red", "green", "blue", "undefined"], "interleave": "pixel", "dtype": "uint8", "driver": "GTiff", "transform": [0.6, 0.0, 405054.0, 0.0, -0.6, 4532580.0], "lnglat": [-124.0937491743777, 40.90625190492096], "height": 12180, "width": 9430, "shape": [12180, 9430], "tiled": false, "res": [0.6, 0.6], "nodata": null, "bounds": [405054.0, 4525272.0, 410712.0, 4532580.0]}

NAIP California 2016, 0.6 meter per pixel, compressed JPEG YCbCr in 4326


Prioritizing LA County in NAIP 2016

  update doqq_processing as a 
    set priority=6 from tl_2016_us_county b, naip_3_16_1_1_ca c 
        b.statefp='06' and b.countyfp='037' and 
        st_intersects(b.geom, c.geom) and a.doqqid=c.gid;


Feb 15 -ARCHIVE-

Jan 17 -ARCHIVE-

2016 ARCHIVE Main -LINK-