In [24]:

    

%pwd


    

    
Out[24]:





u'/home/liux13/Downloads'

In [25]:

    

%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [26]:

    

import os
import laspy
import fiona
import shapely
import shapely.geometry
import rasterio
from pprint import pprint
from tqdm import tnrange


    

In [27]:

    

import rasterio.crs


    

In [28]:

    

from tqdm import tnrange, tqdm
from glob import glob
from time import sleep
import gc


    

In [29]:

    

import traceback


    

In [30]:

    

from tqdm import tqdm_notebook
from multiprocessing import Pool, cpu_count

print cpu_count(), "cpu cores on this machine"


    

    




40 cpu cores on this machine

In [31]:

    

VECTOR_DATA_FN = r'/home/liux13/Desktop/tmp/boxes_section11.shp'
assert os.path.isfile(VECTOR_DATA_FN)


    

In [32]:

    

GROUND_FN = r'/mnt/wwn-0x5000c5008aadf02c-part2/ASU_LiDAR_Processing/section11/merged_ground.asc'
assert os.path.isfile(GROUND_FN)


    

In [33]:

    

OUTPUT_METERS_PER_PIXEL = .05 # five centimeters, about 2 inches
OUTPUT_PADDING = 5 # expand output region 5  meters (to fit a local neighborhood into the image at the edges)


    

In [34]:

    

SLICE_THICKNESS = 0.313


    

In [35]:

    

OUTPUT_FN = './stack.vrt'


    

In [36]:

    

vector_file = fiona.open(VECTOR_DATA_FN)


    

In [37]:

    

minx, miny, maxx, maxy = vector_file.bounds

# Expand the region so that a local windo centered at points on the edges
# does not goi outside the grid
minx -= OUTPUT_PADDING
miny -= OUTPUT_PADDING
maxx += OUTPUT_PADDING
maxy += OUTPUT_PADDING

print "Region of interest (based on vector data)"
print minx, miny, "--", maxx, maxy


    

    




Region of interest (based on vector data)
231638.410976 251008.999357 -- 233216.25982 252568.834735

In [15]:

    

width_in_meters = maxx-minx
height_in_meters = maxy-miny
print 'width (m):', width_in_meters
print 'height (m):', height_in_meters


    

    




width (m): 1577.84884394
height (m): 1559.83537839

In [16]:

    

grid_width = int(round(width_in_meters/OUTPUT_METERS_PER_PIXEL))
grid_height = int(round(height_in_meters/OUTPUT_METERS_PER_PIXEL))
print 'grid width (cols):', grid_width
print 'grid height (rows):', grid_height


    

    




grid width (cols): 31557
grid height (rows): 31197

In [38]:

    

pattern = os.path.join(r'/mnt/wwn-0x5000c5008aadf02c-part2/ASU_LiDAR_Processing/section11/tiled/50/*.laz')
las_files = glob(pattern)
print len(las_files)


    

    




1680

In [18]:

    

def save_stack(stack, las_fn, min_point, max_point, las_width, las_height):
    output_transform = rasterio.transform.from_bounds(min_point[0], min_point[1],
                                                      max_point[0], max_point[1], 
                                                      las_width, las_height)
    output_crs = rasterio.crs.from_epsg(26949)

    output_fn = os.path.split(las_fn)[1][:-3] + 'tif'
    if os.path.isfile(output_fn):
        os.remove(output_fn)
        
    output_file = rasterio.open(output_fn, 'w', 
                            driver=u'GTiff', 
                            crs = output_crs,
                            transform=output_transform,
                            dtype=rasterio.uint16, 
                            count=3, 
                            width=las_width, 
                            height=las_height)
    
    output_file.write(indexes=(1,2,3), src=stack.astype(np.uint16))
    
    output_file.close()
    return output_fn


    

In [39]:

    

def process_las_file(las_fn):
    """Processes a single LAS/LAZ file
    This function either returns the name of a file, or an exception object. 
    It is intended to be used in multiprocwesing.Pool.imap
    """
    output_fn = ''
    gc.collect()
    try:
        las_file = laspy.file.File(las_fn)
        X = las_file.x
        Y = las_file.y
        Z = las_file.z
        min_point = las_file.header.min
        max_point = las_file.header.max
        las_file.close()

        las_min_x = int(floor((min_point[0]-minx) / OUTPUT_METERS_PER_PIXEL))
        las_min_y = int(floor((min_point[1]-miny) / OUTPUT_METERS_PER_PIXEL))
        las_max_x = int(floor((max_point[0]-minx) / OUTPUT_METERS_PER_PIXEL))
        las_max_y = int(floor((max_point[1]-miny) / OUTPUT_METERS_PER_PIXEL))
        las_width = las_max_x - las_min_x+1
        las_height = las_max_y - las_min_y+1

        if max_point[0] < minx or min_point[0] > maxx or max_point[1] < miny or min_point[1] > maxy:
            return "Skipping {}, out-of-bounds file".format(os.path.split(las_fn)[1])

        ground_file = rasterio.open(GROUND_FN)
        ground_min_x, ground_min_y = ~(ground_file.affine)*(min_point[0], min_point[1])
        ground_max_x, ground_max_y = ~(ground_file.affine)*(max_point[0], max_point[1])
        ground_min_x = int(round(ground_min_x))
        ground_min_y = int(round(ground_min_y))
        ground_max_x = int(round(ground_max_x))
        ground_max_y = int(round(ground_max_y))

        flipped = False
        if ground_max_y <  ground_min_y:
            ground_min_y, ground_max_y = ground_max_y, ground_min_y
            flipped = True

        ground_width = ground_max_x-ground_min_x + 1
        ground_height = ground_max_y - ground_min_y + 1

        ground = ground_file.read(1, window=((ground_min_y, ground_max_y+1), (ground_min_x, ground_max_x+1)), 
                                  boundless=True)
        if flipped:
            ground = ground[::-1,:] 
        ground_file.close()

        C = np.floor(((X - minx)/OUTPUT_METERS_PER_PIXEL)).astype(int) - las_min_x
        R = np.floor(((Y - miny)/OUTPUT_METERS_PER_PIXEL)).astype(int) - las_min_y

        height_above_ground = Z - ground[(R*ground_height)/las_height, (C*ground_width)/las_width]      

        level = np.floor((height_above_ground)/SLICE_THICKNESS).astype(int)
        
        def _generate_level(lvl):
            bc = np.bincount(R[level==lvl]*las_width + C[level==lvl], minlength=las_width*las_height)
            output = bc.reshape((las_height, las_width)).astype(float)
            return output

        stack = np.stack([_generate_level(lvl) for lvl in 1,2,3])

        # I seem to be working upside-down...
        output_fn = save_stack(stack[:,::-1,:], las_fn,  min_point, max_point, las_width, las_height)
    except e:
        return (e, traceback.format_exc())
    
    return output_fn


    

In [40]:

    

print process_las_file(las_files[41])
print process_las_file(las_files[-1])


    

    




Skipping tile_231350_251400.laz, out-of-bounds file
Skipping tile_233600_252200.laz, out-of-bounds file

In [41]:

    

jobs = las_files 

for i in 0, 1, 2:
    print '{:04}'.format(i), jobs[i]
print "..."    
for i in len(jobs)-3, len(jobs)-2, len(jobs)-1:
    print '{:04}'.format(i), jobs[i]


    

    




0000 /mnt/wwn-0x5000c5008aadf02c-part2/ASU_LiDAR_Processing/section11/tiled/50/tile_231200_252550.laz
0001 /mnt/wwn-0x5000c5008aadf02c-part2/ASU_LiDAR_Processing/section11/tiled/50/tile_231250_251350.laz
0002 /mnt/wwn-0x5000c5008aadf02c-part2/ASU_LiDAR_Processing/section11/tiled/50/tile_231250_251400.laz
...
1677 /mnt/wwn-0x5000c5008aadf02c-part2/ASU_LiDAR_Processing/section11/tiled/50/tile_233450_252250.laz
1678 /mnt/wwn-0x5000c5008aadf02c-part2/ASU_LiDAR_Processing/section11/tiled/50/tile_233500_252100.laz
1679 /mnt/wwn-0x5000c5008aadf02c-part2/ASU_LiDAR_Processing/section11/tiled/50/tile_233600_252200.laz

In [44]:

    

import gc
gc.collect()


    

    
Out[44]:





0

In [45]:

    

if not os.path.isfile('./try2/stack.vrt'):
    workers = Pool(processes=20)
    try:
        pbar = tqdm_notebook(total=len(jobs))
        completed = []
        for output_fn in workers.imap(process_las_file, jobs):
            print output_fn
            completed.append(output_fn)
            pbar.update(1)
    finally:
        workers.terminate()
        workers.close()
        pbar.close()


    

    




Skipping tile_231200_252550.laz, out-of-bounds file
Skipping tile_231250_251350.laz, out-of-bounds file
Skipping tile_231250_251400.laz, out-of-bounds file
Skipping tile_231250_251750.laz, out-of-bounds file
Skipping tile_231250_251850.laz, out-of-bounds file
Skipping tile_231250_252000.laz, out-of-bounds file
Skipping tile_231250_252050.laz, out-of-bounds file
Skipping tile_231250_252550.laz, out-of-bounds file
Skipping tile_231300_250750.laz, out-of-bounds file
Skipping tile_231300_250800.laz, out-of-bounds file
Skipping tile_231300_250850.laz, out-of-bounds file
Skipping tile_231300_250900.laz, out-of-bounds file
Skipping tile_231300_250950.laz, out-of-bounds file
Skipping tile_231300_251050.laz, out-of-bounds file
Skipping tile_231300_251250.laz, out-of-bounds file
Skipping tile_231300_251350.laz, out-of-bounds file
Skipping tile_231300_251400.laz, out-of-bounds file
Skipping tile_231300_251450.laz, out-of-bounds file
Skipping tile_231300_251500.laz, out-of-bounds file
Skipping tile_231300_251550.laz, out-of-bounds file
Skipping tile_231300_251600.laz, out-of-bounds file
Skipping tile_231300_251700.laz, out-of-bounds file
Skipping tile_231300_251750.laz, out-of-bounds file
Skipping tile_231300_251800.laz, out-of-bounds file
Skipping tile_231300_251850.laz, out-of-bounds file
Skipping tile_231300_252000.laz, out-of-bounds file
Skipping tile_231300_252050.laz, out-of-bounds file
Skipping tile_231300_252100.laz, out-of-bounds file
Skipping tile_231300_252500.laz, out-of-bounds file
Skipping tile_231300_252550.laz, out-of-bounds file
Skipping tile_231350_250700.laz, out-of-bounds file
Skipping tile_231350_250750.laz, out-of-bounds file
Skipping tile_231350_250800.laz, out-of-bounds file
Skipping tile_231350_250850.laz, out-of-bounds file
Skipping tile_231350_250900.laz, out-of-bounds file
Skipping tile_231350_250950.laz, out-of-bounds file
Skipping tile_231350_251000.laz, out-of-bounds file
Skipping tile_231350_251050.laz, out-of-bounds file
Skipping tile_231350_251200.laz, out-of-bounds file
Skipping tile_231350_251250.laz, out-of-bounds file
Skipping tile_231350_251350.laz, out-of-bounds file
Skipping tile_231350_251400.laz, out-of-bounds file
Skipping tile_231350_251450.laz, out-of-bounds file
Skipping tile_231350_251500.laz, out-of-bounds file
Skipping tile_231350_251550.laz, out-of-bounds file
Skipping tile_231350_251650.laz, out-of-bounds file
Skipping tile_231350_251700.laz, out-of-bounds file
Skipping tile_231350_251750.laz, out-of-bounds file
Skipping tile_231350_251800.laz, out-of-bounds file
Skipping tile_231350_251850.laz, out-of-bounds file
Skipping tile_231350_251950.laz, out-of-bounds file
Skipping tile_231350_252000.laz, out-of-bounds file
Skipping tile_231350_252050.laz, out-of-bounds file
Skipping tile_231350_252100.laz, out-of-bounds file
Skipping tile_231350_252200.laz, out-of-bounds file
Skipping tile_231350_252500.laz, out-of-bounds file
Skipping tile_231350_252550.laz, out-of-bounds file
Skipping tile_231350_252700.laz, out-of-bounds file
Skipping tile_231400_250700.laz, out-of-bounds file
Skipping tile_231400_250750.laz, out-of-bounds file
Skipping tile_231400_250850.laz, out-of-bounds file
Skipping tile_231400_250900.laz, out-of-bounds file
Skipping tile_231400_250950.laz, out-of-bounds file
Skipping tile_231400_251000.laz, out-of-bounds file
Skipping tile_231400_251050.laz, out-of-bounds file
Skipping tile_231400_251100.laz, out-of-bounds file
Skipping tile_231400_251150.laz, out-of-bounds file
Skipping tile_231400_251200.laz, out-of-bounds file
Skipping tile_231400_251250.laz, out-of-bounds file
Skipping tile_231400_251300.laz, out-of-bounds file
Skipping tile_231400_251350.laz, out-of-bounds file
Skipping tile_231400_251400.laz, out-of-bounds file
Skipping tile_231400_251450.laz, out-of-bounds file
Skipping tile_231400_251500.laz, out-of-bounds file
Skipping tile_231400_251550.laz, out-of-bounds file
Skipping tile_231400_251650.laz, out-of-bounds file
Skipping tile_231400_251700.laz, out-of-bounds file
Skipping tile_231400_251750.laz, out-of-bounds file
Skipping tile_231400_251800.laz, out-of-bounds file
Skipping tile_231400_251850.laz, out-of-bounds file
Skipping tile_231400_251900.laz, out-of-bounds file
Skipping tile_231400_251950.laz, out-of-bounds file
Skipping tile_231400_252000.laz, out-of-bounds file
Skipping tile_231400_252050.laz, out-of-bounds file
Skipping tile_231400_252100.laz, out-of-bounds file
Skipping tile_231400_252150.laz, out-of-bounds file
Skipping tile_231400_252400.laz, out-of-bounds file
Skipping tile_231400_252500.laz, out-of-bounds file
Skipping tile_231400_252550.laz, out-of-bounds file
Skipping tile_231400_252600.laz, out-of-bounds file
Skipping tile_231400_252700.laz, out-of-bounds file
Skipping tile_231450_250700.laz, out-of-bounds file
Skipping tile_231450_250750.laz, out-of-bounds file
Skipping tile_231450_250800.laz, out-of-bounds file
Skipping tile_231450_250850.laz, out-of-bounds file
Skipping tile_231450_250900.laz, out-of-bounds file
Skipping tile_231450_250950.laz, out-of-bounds file
Skipping tile_231450_251000.laz, out-of-bounds file
Skipping tile_231450_251050.laz, out-of-bounds file
Skipping tile_231450_251100.laz, out-of-bounds file
Skipping tile_231450_251150.laz, out-of-bounds file
Skipping tile_231450_251200.laz, out-of-bounds file
Skipping tile_231450_251250.laz, out-of-bounds file
Skipping tile_231450_251300.laz, out-of-bounds file
Skipping tile_231450_251350.laz, out-of-bounds file
Skipping tile_231450_251450.laz, out-of-bounds file
Skipping tile_231450_251500.laz, out-of-bounds file
Skipping tile_231450_251550.laz, out-of-bounds file
Skipping tile_231450_251600.laz, out-of-bounds file
Skipping tile_231450_251650.laz, out-of-bounds file
Skipping tile_231450_251700.laz, out-of-bounds file
Skipping tile_231450_251750.laz, out-of-bounds file
Skipping tile_231450_251800.laz, out-of-bounds file
Skipping tile_231450_251850.laz, out-of-bounds file
Skipping tile_231450_251900.laz, out-of-bounds file
Skipping tile_231450_251950.laz, out-of-bounds file
Skipping tile_231450_252000.laz, out-of-bounds file
Skipping tile_231450_252050.laz, out-of-bounds file
Skipping tile_231450_252100.laz, out-of-bounds file
Skipping tile_231450_252150.laz, out-of-bounds file
Skipping tile_231450_252250.laz, out-of-bounds file
Skipping tile_231450_252300.laz, out-of-bounds file
Skipping tile_231450_252350.laz, out-of-bounds file
Skipping tile_231450_252400.laz, out-of-bounds file
Skipping tile_231450_252450.laz, out-of-bounds file
Skipping tile_231450_252500.laz, out-of-bounds file
Skipping tile_231450_252550.laz, out-of-bounds file
Skipping tile_231450_252600.laz, out-of-bounds file
Skipping tile_231450_252650.laz, out-of-bounds file
Skipping tile_231450_252700.laz, out-of-bounds file
Skipping tile_231450_252750.laz, out-of-bounds file
Skipping tile_231500_250700.laz, out-of-bounds file
Skipping tile_231500_250750.laz, out-of-bounds file
Skipping tile_231500_250800.laz, out-of-bounds file
Skipping tile_231500_250850.laz, out-of-bounds file
Skipping tile_231500_250950.laz, out-of-bounds file
Skipping tile_231500_251000.laz, out-of-bounds file
Skipping tile_231500_251050.laz, out-of-bounds file
Skipping tile_231500_251100.laz, out-of-bounds file
Skipping tile_231500_251150.laz, out-of-bounds file
Skipping tile_231500_251200.laz, out-of-bounds file
Skipping tile_231500_251250.laz, out-of-bounds file
Skipping tile_231500_251300.laz, out-of-bounds file
Skipping tile_231500_251350.laz, out-of-bounds file
Skipping tile_231500_251400.laz, out-of-bounds file
Skipping tile_231500_251450.laz, out-of-bounds file
Skipping tile_231500_251500.laz, out-of-bounds file
Skipping tile_231500_251550.laz, out-of-bounds file
Skipping tile_231500_251600.laz, out-of-bounds file
Skipping tile_231500_251650.laz, out-of-bounds file
Skipping tile_231500_251750.laz, out-of-bounds file
Skipping tile_231500_251800.laz, out-of-bounds file
Skipping tile_231500_251850.laz, out-of-bounds file
Skipping tile_231500_251900.laz, out-of-bounds file
Skipping tile_231500_251950.laz, out-of-bounds file
Skipping tile_231500_252000.laz, out-of-bounds file
Skipping tile_231500_252050.laz, out-of-bounds file
Skipping tile_231500_252100.laz, out-of-bounds file
Skipping tile_231500_252150.laz, out-of-bounds file
Skipping tile_231500_252200.laz, out-of-bounds file
Skipping tile_231500_252250.laz, out-of-bounds file
Skipping tile_231500_252300.laz, out-of-bounds file
Skipping tile_231500_252350.laz, out-of-bounds file
Skipping tile_231500_252400.laz, out-of-bounds file
Skipping tile_231500_252450.laz, out-of-bounds file
Skipping tile_231500_252550.laz, out-of-bounds file
Skipping tile_231500_252600.laz, out-of-bounds file
Skipping tile_231500_252650.laz, out-of-bounds file
Skipping tile_231500_252700.laz, out-of-bounds file
Skipping tile_231550_250600.laz, out-of-bounds file
Skipping tile_231550_250700.laz, out-of-bounds file
Skipping tile_231550_250750.laz, out-of-bounds file
Skipping tile_231550_250800.laz, out-of-bounds file
Skipping tile_231550_250850.laz, out-of-bounds file
Skipping tile_231550_250900.laz, out-of-bounds file
Skipping tile_231550_250950.laz, out-of-bounds file
Skipping tile_231550_251000.laz, out-of-bounds file
Skipping tile_231550_251050.laz, out-of-bounds file
Skipping tile_231550_251100.laz, out-of-bounds file
Skipping tile_231550_251150.laz, out-of-bounds file
Skipping tile_231550_251250.laz, out-of-bounds file
Skipping tile_231550_251300.laz, out-of-bounds file
Skipping tile_231550_251350.laz, out-of-bounds file
Skipping tile_231550_251400.laz, out-of-bounds file
Skipping tile_231550_251450.laz, out-of-bounds file
Skipping tile_231550_251500.laz, out-of-bounds file
Skipping tile_231550_251550.laz, out-of-bounds file
Skipping tile_231550_251600.laz, out-of-bounds file
Skipping tile_231550_251650.laz, out-of-bounds file
Skipping tile_231550_251700.laz, out-of-bounds file
Skipping tile_231550_251750.laz, out-of-bounds file
Skipping tile_231550_251800.laz, out-of-bounds file
Skipping tile_231550_251850.laz, out-of-bounds file
Skipping tile_231550_251900.laz, out-of-bounds file
Skipping tile_231550_251950.laz, out-of-bounds file
Skipping tile_231550_252050.laz, out-of-bounds file
Skipping tile_231550_252100.laz, out-of-bounds file
Skipping tile_231550_252150.laz, out-of-bounds file
Skipping tile_231550_252200.laz, out-of-bounds file
Skipping tile_231550_252250.laz, out-of-bounds file
Skipping tile_231550_252300.laz, out-of-bounds file
Skipping tile_231550_252350.laz, out-of-bounds file
Skipping tile_231550_252400.laz, out-of-bounds file
Skipping tile_231550_252450.laz, out-of-bounds file
Skipping tile_231550_252500.laz, out-of-bounds file
Skipping tile_231550_252550.laz, out-of-bounds file
Skipping tile_231550_252600.laz, out-of-bounds file
Skipping tile_231550_252650.laz, out-of-bounds file
Skipping tile_231550_252700.laz, out-of-bounds file
Skipping tile_231600_250550.laz, out-of-bounds file
Skipping tile_231600_250700.laz, out-of-bounds file
Skipping tile_231600_250750.laz, out-of-bounds file
Skipping tile_231600_250800.laz, out-of-bounds file
Skipping tile_231600_250850.laz, out-of-bounds file
Skipping tile_231600_250900.laz, out-of-bounds file
Skipping tile_231600_250950.laz, out-of-bounds file







    




---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-45-0de4aa946d5c> in <module>()
      4         pbar = tqdm_notebook(total=len(jobs))
      5         completed = []
----> 6         for output_fn in workers.imap(process_las_file, jobs):
      7             print output_fn
      8             completed.append(output_fn)

/home/liux13/.conda/envs/box/lib/python2.7/multiprocessing/pool.pyc in next(self, timeout)
    666         if success:
    667             return value
--> 668         raise value
    669 
    670     __next__ = next                    # XXX

AttributeError: 'module' object has no attribute 'from_epsg'

In [48]:

    

!gdalbuildvrt stack.vrt tile_*.tif


    

    




0...10...20...30...40...50...60...70...80...90...100 - done.
ERROR 4: tile_*.tif: No such file or directory
Warning 1: Can't open tile_*.tif. Skipping it

In [49]:

    

output_ds = rasterio.open('try2/stack.vrt')


    

    




---------------------------------------------------------------------------
RasterioIOError                           Traceback (most recent call last)
<ipython-input-49-c3d949a462d3> in <module>()
----> 1 output_ds = rasterio.open('try2/stack.vrt')

/home/liux13/.conda/envs/box/lib/python2.7/site-packages/rasterio/__init__.pyc in open(path, mode, driver, width, height, count, crs, transform, dtype, nodata, **kwargs)
    191         raise ValueError(
    192             "mode string must be one of 'r', 'r+', or 'w', not %s" % mode)
--> 193     s.start()
    194     return s
    195 

rasterio/_base.pyx in rasterio._base.DatasetReader.start (rasterio/_base.c:2980)()

RasterioIOError: try2/stack.vrt: No such file or directory

In [ ]:

    

vector_file = fiona.open('./boxes11/boxes_section11.json')
print vector_file.meta


    

In [ ]:

    

hotspots = [shapely.geometry.shape(f['geometry']).centroid for f in vector_file if f['geometry'] is not None]


    

In [ ]:

    

hs = hotspots[randint(len(hotspots))]
ctr = np.concatenate(hs.xy)


    

In [ ]:

    

print ctr
cr = np.asarray(~output_ds.affine * ctr).astype(int)
c, r = cr
print cr

data = output_ds.read(window=((r-32, r+32), (c-32, c+32)), boundless=True)

imshow((data.transpose(1,2,0)/20.).clip(0, 1))


    

In [ ]:

    

rgb = rasterio.open('./rgb/rgb.vrt')


    

In [ ]:

    

import pyproj


    

In [ ]:

    

EPSG2223 = pyproj.Proj(init="epsg:2223", preserve_units=True)
EPSG26949 = pyproj.Proj(init="epsg:26949", preserve_units=True)
print EPSG2223.srs
print EPSG26949.srs


    

In [ ]:

    

rgb_ctr = pyproj.transform(EPSG26949, EPSG2223, ctr[0], ctr[1])
print ctr, rgb_ctr


    

In [ ]:

    

rgb_cr = np.asarray(~rgb.affine * rgb_ctr).astype(int)
print rgb_cr
rgb_data = rgb.read(window=((rgb_cr[1]-8, rgb_cr[1]+8), (rgb_cr[0]-8, rgb_cr[0]+8)), boundless=True)
imshow(rgb_data.transpose(1,2,0))


    

In [ ]:

    

%%file data_provider.py

import pyproj
import rasterio
import numpy as np
import scipy.ndimage
from scipy.ndimage import zoom
from math import cos, sin, radians, atan2, degrees
from collections import namedtuple

Predictions = namedtuple('Predictions', ['x', 'y', 'label', 'dx', 'dy', 'angle'])
Sample = namedtuple('Sample', ['data', 'predicted'])

EPSG2223 = pyproj.Proj(init="epsg:2223", preserve_units=True)
EPSG26949 = pyproj.Proj(init="epsg:26949", preserve_units=True)
        
def normalize_angle(angle):
    angle = radians(angle)
    angle = atan2(sin(angle), cos(angle))
    angle = degrees(angle)
    return angle
    
class DataProvider(object):
    def __init__(self, radius_in_pixels=32, jitter=0.25):        
        self.densities_path = '/home/shared/srp/try2/stack.vrt'
        self.colors_path = '/home/shared/srp/rgb/rgb.vrt'
        self.sample_path = '/home/shared/srp/sample_locations_epsg26949.npz'
        self.radius_in_pixels = radius_in_pixels
        self.jitter = jitter
        self._open_datasets()
      
    def _open_datasets(self):
        self.densities = rasterio.open(self.densities_path)
        self.colors = rasterio.open(self.colors_path)
        samples= np.load(self.sample_path)
        self.pos_xy = samples['pos_xy']
        self.neg_xy = samples['neg_xy']
        self.pos_angles = samples['pos_angles']
    
    def random_positive(self, jitter=None, rotate=True):
        jitter = self.jitter if jitter is None else jitter
        index = np.random.randint(len(self.pos_xy))
        dx, dy = jitter*np.random.randn(2)
        x, y = self.pos_xy[index] 
        x -= dx
        y -= dy
        label=1
        if rotate:
            rotation = np.random.randint(0, 360)
            patch = self.get_patch_xyr(x, y, rotation)        
            angle = normalize_angle(rotation + self.pos_angles[index])
        else:
            patch = self.get_patch_xy(x, y)
            angle = self.pos_angles[index]
            
        return Sample(patch, Predictions(x, y, label, dx, dy, angle))
    
    def random_negative(self, jitter=None, rotate=True):
        jitter = self.jitter if jitter is None else jitter
        index = np.random.randint(len(self.neg_xy))
        dx, dy = jitter*np.random.randn(2)
        x, y = self.neg_xy[index] 
        x -= dx
        y -= dy
        label = 0
        
        if rotate:
            rotation = np.random.randint(0, 360)
            patch = self.get_patch_xyr(x, y, rotation)        
            angle = normalize_angle(rotation)
        else:
            patch = self.get_patch_xy(x, y)
            angle = 0
            
        return Sample(patch, Predictions(x, y, label, dx, dy, angle))
    
    def random_sample(self, jitter=None, rotate=True, prob_positive=0.5):
        jitter = self.jitter if jitter is None else jitter
        if np.random.rand() > prob_positive:
            return self.random_negative(jitter, rotate)
        else:
            return self.random_positive(jitter, rotate)
        
    def get_patch_xyr(self, x, y, angle, radius_in_pixels=None):
        radius_in_pixels = radius_in_pixels or self.radius_in_pixels
        
        source_patch = self.get_patch_xy(x, y, radius_in_pixels*2)
        width = height = 2*radius_in_pixels
                
        radians = np.radians(angle)
        c, s = cos(radians), sin(radians)
        R = np.matrix([[c, -s], 
                       [s, c]])
        X = np.asarray([width, height])
        X = np.asarray(X-R.dot(X)).flatten()
        
        rotated_patch = np.empty_like(source_patch)
        for i in range(len(rotated_patch)):
            scipy.ndimage.affine_transform(source_patch[i],
                                           matrix=R, offset=X, 
                                           output_shape = rotated_patch[i].shape,
                                           output=rotated_patch[i])
        
        x, y = int((source_patch.shape[2]-width)/2), int( (source_patch.shape[1]-height)/2)
        cropped_patch = rotated_patch[:, y:y+height, x:x+width].copy()
        return cropped_patch    
        
    
    def get_patch_xy(self, x, y, radius_in_pixels=None):
        radius_in_pixels = radius_in_pixels or self.radius_in_pixels
        R = radius_in_pixels
        x_2223, y_2223 = pyproj.transform(EPSG26949, EPSG2223, x, y)
        c_2223, r_2223 = np.asarray(~self.colors.affine * (x_2223, y_2223)).astype(int)
        c, r = np.asarray(~self.densities.affine*(x, y)).astype(int)
                
        colors = self.colors.read(window=((r_2223-R, r_2223+R), (c_2223-R, c_2223+R)), 
                                  boundless=True).astype(np.float32)/255.
        
        densities = self.densities.read(window=((r-R, r+R), (c-R, c+R)),
                                        boundless=True).astype(np.float32)
        
        combined = np.concatenate([colors, densities])
        
        return combined


    

In [46]:

    

import data_provider
reload(data_provider)
DataProvider = data_provider.DataProvider


    

    




---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-46-14f1eeaf372f> in <module>()
----> 1 import data_provider
      2 reload(data_provider)
      3 DataProvider = data_provider.DataProvider

ImportError: No module named data_provider

In [ ]:

    

from pandas import DataFrame


    

In [ ]:

    

gtp = DataProvider(radius_in_pixels=32, jitter = 0.25)


    

In [ ]:

    

patch, labels = gtp.random_positive()
print patch.shape

subplot(1,2,1)
imshow(patch[:3,:,:].transpose(1,2,0).clip(0,1))
title('colors')

subplot(1,2,2)
imshow((patch[3:6,:,:]/40).transpose(1,2,0).clip(0,1))
title('densities');

display(DataFrame([labels], columns=data_provider.Predictions._fields))


    

In [ ]:

    

patch, labels = gtp.random_sample()

subplot(1,2,1)
imshow(patch[:3,:,:].transpose(1,2,0).clip(0,1))
title('colors')

subplot(1,2,2)
imshow((patch[3:6,:,:]/20).transpose(1,2,0).clip(0,1))
title('densities')

display(DataFrame([labels], columns=data_provider.Predictions._fields))


    

In [ ]:

    

def make_datum(patch, labels):
     #Create a Datum
    datum = Datum()
    datum.channels, datum.height, datum.width = patch.shape
    datum.data = patch.tobytes()
    datum.label = labels.label    
    return datum


    

In [ ]:

    

from tqdm import tnrange
import lmdb
from caffe.proto.caffe_pb2 import Datum

def make_lmdb(dbname, nsamples, data_generator):
    train_data_env = lmdb.open(dbname)
    train_data_env.set_mapsize(1e12)

    for i in tnrange(nsamples):
        sample = data_generator.random_sample()
        datum = make_datum(sample.data, sample.predicted)
        key = '{:06}{:06}'.format(int(100*sample.predicted.x), int(100*sample.predicted.y))
        with train_data_env.begin(write=True) as txn:
            print '\r', key,
            txn.put(key, datum.SerializeToString())


    

In [ ]:

    

make_lmdb('/home/shared/srp/ds1.lmdb', 1000, gtp)


    

In [ ]:

    

import caffe, caffe.net_spec


    

In [ ]:

    

from caffe import layers as L
from caffe import params as P


    

In [ ]:

    

batch_size=64
n = caffe.NetSpec()
n.data, n.label = L.DummyData(shape=[dict(dim=[batch_size, 1, 28, 28]),
                                     dict(dim=[batch_size, 1, 1, 1])],
                              transform_param=dict(scale=1./255), ntop=2)

n.conv1 = L.Convolution(n.data,  kernel_size=5, num_output=20, weight_filler=dict(type='xavier'))
n.pool1 = L.Pooling(n.conv1,  kernel_size=2, stride=2, pool=P.Pooling.MAX) 
n.conv2 = L.Convolution(n.pool1, kernel_size=5, num_output=50, weight_filler=dict(type='xavier'))
n.pool2 = L.Pooling(n.conv2,  kernel_size=2, stride=2, pool=P.Pooling.MAX)
n.ip1 = L.InnerProduct(n.pool2,  num_output=500, weight_filler=dict(type='xavier'))
n.relu1 = L.ReLU(n.ip1,  in_place=True)
n.ip2 = L.InnerProduct(n.relu1,  num_output=10, weight_filler=dict(type='xavier'))
n.loss = L.SoftmaxWithLoss(n.ip2, n.label, )


    

In [47]:

    

from caffe.draw import draw_net


    

In [ ]:

    

import pydot

g = caffe.draw.get_pydot_graph(n.to_proto(), 'TB')

from IPython.display import Image, display
Image(g.create_png())


    

In [ ]:

    

from IPython.display import Image, display
Image(g.create_png())