<jupyter_start><jupyter_text>This script is used to conver the point level data into the street-segmentUsing the street-segment rather than the point level data would be much better for the data visualization.Copyright(C) <NAME>, MIT Senseable City Lab, Nov 21, 2018<jupyter_code>import os, os.path
from PIL import Image
import numpy as np
from matplotlib import pyplot as plt


root = '/home/xiaojiang/xiaojiang/tokyo-proj/shibuya-proj'
segImgfiles = os.path.join(root, 'gsv-images/tif-seg-panos')
segImgfile = os.path.join(segImgfiles, '9 - TwUGzndW9KW4tdupf0c6Hg - 139.709491267 - 35.662483952 - 2013-06 - 53.015842437_sky.png')

segImg = np.asarray(Image.open(segImgfile))

plt.imshow(segImg)
plt.show()

np.unique(segImg)<jupyter_output><empty_output><jupyter_text>Generate the hemispherical images and rotate the hemispherical images to make sure the north direction is up<jupyter_code>import sys
sys.path.append("libraries")
import SunExpoLib as sunexpo


segHemiImgRoot = os.path.join(root, 'gsv-images/seg-hemi')

if not os.path.exists(segHemiImgRoot):
    os.mkdir(segHemiImgRoot)
    print('The segHemi folder is not existing, we are going to make a folder')
    

i = 0
for imgfile in os.listdir(segImgfiles):
    basename = imgfile.split('_sky.')[0]
    yaw = basename.split(' - ')[-1]
    
    # the output seg hemifile name
    segHemiImgName = basename + '.png'
    segHemiImgFullName = os.path.join(segHemiImgRoot, segHemiImgName)
    
    segPanoImg = np.asarray(Image.open(os.path.join(segImgfiles, imgfile)))
    segHemiImg = sunexpo.cylinder2fisheyeImage(segPanoImg, yaw, segHemiImgFullName)
    
    if os.path.exists(segHemiImgFullName): continue
    
    i = i + 1
    if i % 10000 == 0: print('You have convereted %s pano into hemispherical images' %(i))
    
#     plt.imshow(segPanoImg)
#     plt.show()
    
#     plt.imshow(segHemiImg)
#     plt.show()
#     break<jupyter_output><empty_output><jupyter_text>Example to check whether the direction is right<jupyter_code>panoId = '97NIppgAnsJXwAv0XFu0Yw'
cylinderPanoImg = sunexpo.GSVpanoramaDowloader_GoogleMaps(panoId, '97NIppgAnsJXwAv0XFu0Yw.png')

fisheyeImg = sunexpo.cylinder2fisheyeImage(cylinderPanoImg, 78.730598449, '97NIppgAnsJXwAv0XFu0Yw_hemi.png')

plt.imshow(cylinderPanoImg)
plt.show()

plt.imshow(fisheyeImg)
plt.show()<jupyter_output>The panoid is:------ 97NIppgAnsJXwAv0XFu0Yw
The size of the GSV panorama is: (832, 1664, 3)<jupyter_text>Calculate the spatial temporal distribution of the sunlight exposure<jupyter_code>import sys
sys.path.append("libraries")
import SunposLib as sunpos


year = 2018
month = 7
day = 15
hour = 9
minute = 0
second = 0
zone = -9 # https://www.esrl.noaa.gov/gmd/grad/solcalc/azel.html


latitude = 35.668263
longitude = -139.697001 # watch out the longitude for the west hemisphere is positive and east is negative


for segImgfile in os.listdir(segHemiImgRoot):
    print (segImgfile)
    segImg = np.asarray(Image.open(os.path.join(segHemiImgRoot, segImgfile)))
    plt.imshow(segImg)
    plt.show()
    
    # plot the result on the hemispherical images
    sunexpo.plot_SunPathOnFisheyeimage_noaa(segImg, 'hemiImg_sunpath.png', latitude, longitude, zone, year, month, day)
    
    plotImg = np.asarray(Image.open('hemiImg_sunpath.png'))
    plt.imshow(plotImg)
    plt.show()
    
    
    # Plot the sun pos on the hemispherical image
    [azimuth, sunele] = sunpos.calcSun(latitude, longitude, zone, year, month, day, hour, minute, second)
    # print ('The hour, sunele and azimuth are:=======', hh, sunele, azimuth)
    [hemiImg, px, py] = sunexpo.plot_SunPosOnFisheyeimage_noaa(segImg, 'hemiImg_sunpos.png', azimuth, sunele)
    
    hemiSunposImg = np.asarray(Image.open("hemiImg_sunpos.png"))
    plt.imshow(hemiSunposImg)
    plt.show()
    
    
    # Judge whether the sunlight is blocked or not, 0 is not shaded, 1 is shaded. 
    shade = sunexpo.Shaded_judgement_noaa(segImg, 0, 4, azimuth, sunele)
    print('This site shade is:', shade)
    
    
    # delete the temperary image file
    if os.path.exists('hemiImg_sunpath.png'): os.remove('hemiImg_sunpath.png')
    if os.path.exists('hemiImg_sunpos.png'): os.remove('hemiImg_sunpos.png')
    
    
    break
    
    
# for hour in range(5, 20):
#     [azimuth, sunele] = sunpos.calcSun(latitude, longitude, zone, year, month, day, hour, minute, second)
#     print ('The sun azimuth and sunele are:', hour, minute, azimuth, sunele)<jupyter_output>9473 - wmb9ClKHbgFaowPhzbd_5Q - 139.70188071 - 35.6675999707 - 2012-04 - 183.45747375.png<jupyter_text>Loop all segmented images and calculate the sunlight exposure<jupyter_code>root = '/home/xiaojiang/xiaojiang/tokyo-proj/shibuya-proj'

segImgroot = os.path.join(root, 'gsv-images/tif-seg-panos')

# for segImgfile in os.listdir(segImgroot):
#     segImgfilename = os.path.join(segImgroot, segImgfile)
#     segImg = np.array(Image.open(segImgfilename))
    

from PIL import Image
import numpy as np
import os, os.path
import math
import ogr,osr
import sys
sys.path.append("libraries")
import SunposLib as sunpos
import SunExpoLib as sunexpo


year = 2018
month = 7
day = 15
hour = 9
minute = 0
second = 0
zone = -9 # https://www.esrl.noaa.gov/gmd/grad/solcalc/azel.html


# watch out the longitude for the west hemisphere is positive and east is negative
latitude = 35.668263
longitude = -139.697001

segHemiImgRoot = os.path.join(root, 'gsv-images/seg-hemi')

# the output sunexpo shapefile
shpfile = r'/home/xiaojiang/xiaojiang/tokyo-proj/shibuya-proj/spatial-data/sunexpo-July-15-9am.shp'

# create a shpafile to save the sun duration
driver = ogr.GetDriverByName("ESRI Shapefile")

if os.path.exists(shpfile):
    driver.DeleteDataSource(shpfile)
data_source = driver.CreateDataSource(shpfile)

targetSpatialRef = osr.SpatialReference()
targetSpatialRef.ImportFromEPSG(4326)

outLayer = data_source.CreateLayer("Sunexpo", targetSpatialRef, ogr.wkbPoint)
panoId = ogr.FieldDefn('panoid', ogr.OFTString)
outLayer.CreateField(panoId)

# create fields for different time stamps
hours = [9, 12, 14, 17]
for hour in hours:
    fieldname = 'expo%s'%(hour)
    exposureField = ogr.FieldDefn(fieldname, ogr.OFTInteger) # 0 means not shaded or exposed to sunlight, 1 mean shaded or not exposed to sunlight
    outLayer.CreateField(exposureField)


# loop all the sky image in the folder
i = 0
for skyimgfile in os.listdir(segHemiImgRoot):
    i = i + 1
    if i % 1000 == 0: print('You have processed %s'%(i))
    
    
    skyImgFileFullname = os.path.join(segHemiImgRoot, skyimgfile)
    skyImg = np.asarray(Image.open(skyImgFileFullname))
    
    basename = skyimgfile.split('_sky.')[0]
    fields = basename.split(' - ')
    
    panoID = fields[1]
    lon = float(fields[2]) # watch out in the NOAA sunpos model, the eastern hemisphere has longitude negative, west has longitude positive
    lat = float(fields[3])
    
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(lon, lat)
    
    # create feature and set the attribute values
    featureDefn = outLayer.GetLayerDefn()
    outFeature = ogr.Feature(featureDefn)
    
    # set the geometrical feature
    outFeature.SetGeometry(point)
    outFeature.SetField('panoid', panoID)
    
    # calculate the several hours sun exposure
    for hour in hours:
        # calculate the sun position
        [azimuth, sunele] = sunpos.calcSun(latitude, longitude, zone, year, month, day, hour, minute, second)
        
        # Judge whether the sunlight is blocked or not, 0 is not shaded, 1 is shaded. 
        shade = sunexpo.Shaded_judgement_noaa(skyImg, 0, 4, azimuth, sunele)
        exposure = 1 - shade # exposed to sunlight 1, not exposed to sunlight 0
        fieldname = 'expo%s'%(hour)
        outFeature.SetField(fieldname, exposure)
        
    
    # push the feature to the layer
    outLayer.CreateFeature(outFeature)
    outFeature.Destroy()
    
    
data_source.Destroy()<jupyter_output>You have processed 1000
You have processed 2000
You have processed 3000
You have processed 4000
You have processed 5000
You have processed 6000
You have processed 7000
You have processed 8000
You have processed 9000
You have processed 10000
You have processed 11000
You have processed 12000
You have processed 13000
You have processed 14000
You have processed 15000
You have processed 16000
You have processed 17000
You have processed 18000
You have processed 19000
You have processed 20000
You have processed 21000
You have processed 22000
You have processed 23000
You have processed 24000
You have processed 25000
You have processed 26000
You have processed 27000
You have processed 28000
You have processed 29000
You have processed 30000
You have processed 31000
You have processed 32000
You have processed 33000
You have processed 34000
You have processed 35000
You have processed 36000
You have processed 37000
You have processed 38000
You have processed 39000
You have processed 40000
You have [...]<jupyter_text>This script is used to aggregate the point level features to linear, which would look better for visualization and analysis<jupyter_code>import os, os.path
import sys
import fiona
from shapely.geometry import LineString, mapping, shape, Point
from shapely.ops import transform, unary_union
import rtree

root = r'/home/xiaojiang/xiaojiang/tokyo-proj/shibuya-proj/spatial-data/'

# the output sunexpo shapefile
shpfile = os.path.join(root, 'sunexpo-July-15-9am-WGS84.shp')
lineshpFile = os.path.join(root, 'Japan_highway_shibuya.shp')

datasource = fiona.open(shpfile, 'r')

crs = datasource.crs
datasource.schema<jupyter_output><empty_output><jupyter_text>Save the point feature class information into list<jupyter_code>expo_records = []
fieldlist = ['panoid', 'expo9', 'expo12', 'expo14', 'expo17']

with fiona.open(shpfile, 'r') as expo:
    for feat in expo:
        item = {}
        for field in fieldlist:
            attr = feat['properties'][field]
            item.update({field: attr})
            
#         lon = feat['properties']['longitude']
#         lat = feat['properties']['latitude']
#         point = Point(float(lon), float(lat))
#         print ('The created geom object is:', point)
        
        # add the point geom to the dictionary
        pnt_geom = shape(feat['geometry'])
        item.update({'point': pnt_geom})
        
        expo_records.append(item)


len(expo_records)<jupyter_output><empty_output><jupyter_text>Split the roads into shorter linesThe following code recursively splits the line in half until it is shorter than the maximum desired length.In order make the result more reasonable, we need to split the long street segment into shorter street segments, cite from https://www.azavea.com/blog/2016/10/05/philippines-road-safety-using-shapely-fiona-locate-high-risk-traffic-areas/<jupyter_code># using recursive statement to split longer roads into shorter segments
# Reference: https://www.azavea.com/blog/2016/10/05/philippines-road-safety-using-shapely-fiona-locate-high-risk-traffic-areas/

def split_line(line, max_line_units):
    '''
    The input:
        line: the LineString object
        max_line_units: the split distance, be careful of the units
    output:
        a list of LineString segment
    '''
    
    if line.length <= max_line_units:
        return [line]
    
    half_length = line.length / 2
    coords = list(line.coords)
    
    for idx, point in enumerate(coords):
        proj_dist = line.project(Point(point))
        if proj_dist == half_length:
            return [LineString(coords[:idx + 1]), LineString(coords[idx:])]
        
        if proj_dist > half_length:
            mid_point = line.interpolate(half_length)
            head_line = LineString(coords[:idx] + [(mid_point.x, mid_point.y)])
            tail_line = LineString([(mid_point.x, mid_point.y)] + coords[idx:])
            return split_line(head_line, max_line_units) + split_line(tail_line, max_line_units)<jupyter_output><empty_output><jupyter_text>Split lines into shorter street segments, keep the attribute of streetThis session also output the shorter street segments into a shapefile, although this is not necessary. The returned results will be saved into a list split_lines, which stores all LineString items.There is one more option, you can go to the SamplingLib.py to find the function of "createStreetSegment"<jupyter_code>def SplitStreets2Segments(streetMap, segStreetMap, distance):
    '''
    This function is used to split the streets into shorter street segments for visualization
    and spatial analyses and modeling purposes
    
    parameters:
        streetMap: the original street feature class
        segStreetMap: the filename of the output street segment
        distance: the distance of the spliting the distance
    
    example:
        root = r'/home/xiaojiang/xiaojiang/tokyo-proj/shibuya-proj/spatial-data/'
        streetMap = os.path.join(root, 'Japan_highway_shibuya.shp')
        segStreetMap = os.path.join(root, 'seg_Japan_highway_shibuya.shp')
        
        SplitStreets2Segments(streetMap, segStreetMap, 40)
    
    
    Copyright(c) <NAME>, MIT Senseable City Lab 
    Last modified Dec 2ed, 2018
    '''
    
    import pyproj
    from shapely.ops import transform
    from functools import partial
    
    split_lines = []
    
    with fiona.open(streetMap, 'r') as streets:
        crs = streets.crs
        schema = streets.schema.copy()
        
        with fiona.open(segStreetMap, 'w', driver="ESRI Shapefile", crs=crs, schema=schema) as output_segment:
            for street in streets:            
                street_attr = street['properties']            
                street_geom = street['geometry']
                featureType = street_geom['type']
                
                # check if the linear feature is LineString or Multi-Line String
                if featureType == 'MultiLineString':
                    print('This is a MultiLineString')
                    multi_street_geom_degree = shape(street_geom)
                    
                    for singleLine_degree in multi_street_geom_degree:
                        length = singleLine_degree.length
                        # partition each single line in the multiline
                        project = partial(pyproj.transform,pyproj.Proj(init='EPSG:4326'),pyproj.Proj(init='EPSG:3857')) #3857 is psudo WGS84 the unit is meter
                        singleLine_meter = transform(project, singleLine_degree)
                        
                        street_segments = split_line(singleLine_meter, distance)
                        split_lines.extend(street_segments)
                        
                        # save the splitted segment to a shapefile
                        for street_segment in street_segments:                
                            project2 = partial(pyproj.transform,pyproj.Proj(init='EPSG:3857'),pyproj.Proj(init='EPSG:4326'))
                            street_segment_deg = transform(project2, street_segment)

                            output_segment.write({
                                'geometry': mapping(street_segment_deg),
                                'properties': dict(street_attr)
                            })
                
                elif featureType == 'LineString': 
                    # RouteName = street['properties']['FULLNAME']
                    street_geom_degree = shape(street_geom)
                    
                    # convert degree to meter, in order to split by distance in meter
                    project = partial(pyproj.transform,pyproj.Proj(init='EPSG:4326'),pyproj.Proj(init='EPSG:3857')) #3857 is psudo WGS84 the unit is meter
                    street_geom_meter = transform(project, street_geom_degree)
                    
                    street_segments = split_line(street_geom_meter, distance)
                    split_lines.extend(street_segments)
                    
                    # save the splitted segment to a shapefile
                    for street_segment in street_segments:                
                        project2 = partial(pyproj.transform,pyproj.Proj(init='EPSG:3857'),pyproj.Proj(init='EPSG:4326'))
                        street_segment_deg = transform(project2, street_segment)
                        
                        output_segment.write({
                            'geometry': mapping(street_segment_deg),
                            'properties': dict(street_attr) # this is important statement to save all the previous properties
                        })
                        
                        
    return split_lines
    
root = r'/home/xiaojiang/xiaojiang/tokyo-proj/shibuya-proj/spatial-data/'
streetMap = os.path.join(root, 'Japan_highway_shibuya.shp')
segStreetMap = os.path.join(root, 'seg_Japan_highway_shibuya.shp')

streetMap = os.path.join(root, 'Shibuya_streetNetwork_noMotorway.shp')
segStreetMap = os.path.join(root, 'seg_Shibuya_streetNetwork_noMotorway.shp')

split_lines = SplitStreets2Segments(streetMap, segStreetMap, 100)

# keys = street_attr.keys()
# att = dict(street_attr)
# print('The attribute keys are:', keys)
# print('The dictionary created is:', att)

file = fiona.open(segStreetMap)
file.schema
print('the crs of the output shapefile is:', file.crs)<jupyter_output>the crs of the output shapefile is: {'init': 'epsg:4326'}<jupyter_text>Create linear segments with buffer distance of 50 meter, because the split_line records the geo in projection of EPSG:3857, the unit is meter<jupyter_code>segments_rtree_index = rtree.index.Index()

i = 0
for idx, element in enumerate(split_lines):
    i = i + 1
    segments_rtree_index.insert(idx, element.buffer(50).bounds)
    
# segments_rtree_index

print('The number of element is:', i)

element.buffer(1).bounds
# print('The boundary of the expo site is:', element.buffer(50))
element.buffer(50)


expo_records[0]['point'].bounds
# expo_records_deg.bounds<jupyter_output><empty_output><jupyter_text>Using r-tree to do the spatial operations<jupyter_code>from functools import partial
import pyproj

# model to convert the unit of WGS84 from degree to meter
wgs_deg2meter_model = partial(pyproj.transform,pyproj.Proj(init='EPSG:4326'),pyproj.Proj(init='EPSG:3857'))

# use a dictionary to store the intersected point features
segments_with_records = {}

for record in expo_records:
    record_pnt = record['point']
    record_pnt_meter = transform(wgs_deg2meter_model, record_pnt)
    
    # use rtree index to find the intersected roads
    inter_roads = segments_rtree_index.intersection(record_pnt_meter.bounds)
    road_id_inter = [(segment_id, split_lines[segment_id].distance(record_pnt_meter)) for segment_id in inter_roads]
    
    
    # find the neareast pnts
    if len(road_id_inter):
        nearest = min(road_id_inter, key=lambda tup: tup[1])
        segment_id = nearest[0]
        
        # create an empty list to store the sunglare records for all street segments
        if segment_id not in segments_with_records:
            segments_with_records[segment_id] = []
        
        segments_with_records[segment_id].append(record)
        

len(segments_with_records)
# segments_with_records[22674]<jupyter_output><empty_output><jupyter_text>Create output shapefile based on the intersection result using R-tree<jupyter_code>import numpy as np

outputShp = os.path.join(root, 'SunExpo_shibuya.shp')
schema = {
    'geometry': 'LineString',
    'properties': {
        'expo9': 'int',
        'expo12': 'int',
        'expo14': 'int',
        'expo17': 'int'
    }
}


with fiona.open(outputShp, 'w', driver="ESRI Shapefile", schema=schema, crs=crs) as output:
    for idx, segment in enumerate(split_lines):
        segment_geom = segment
        
        records = segments_with_records.get(idx)
        
        if records is None: #6 No GSV point on street segment, then assign them to -1
            sunexpo9_ave = -1
            sunexpo12_ave = -1
            sunexpo14_ave = -1
            sunexpo17_ave = -1
        else: 
            # create empty list of the sunrise and sunset glare
            expo9_list = []
            expo12_list = []
            expo14_list = []
            expo17_list = []

            # loop the records to calculate the mean sunglare duration
            for record in records:
                expo9_list.append(record['expo9'])
                expo12_list.append(record['expo12'])
                expo14_list.append(record['expo14'])
                expo17_list.append(record['expo17'])

            sunexpo9_ave = int(sum(expo9_list)/float(len(expo9_list) + 0.0001) + 0.5)
            sunexpo12_ave = int(sum(expo12_list)/float(len(expo12_list) + 0.0001) + 0.5)
            sunexpo14_ave = int(sum(expo14_list)/float(len(expo14_list) + 0.0001) + 0.5)
            sunexpo17_ave = int(sum(expo17_list)/float(len(expo17_list) + 0.0001) + 0.5)
    #         print('The sunlight exposure is:', sunexpo9_ave)
        
        # the data of the record
        data = {
            'expo9': sunexpo9_ave,
            'expo12': sunexpo12_ave,
            'expo14': sunexpo14_ave,
            'expo17': sunexpo17_ave
        }
        
        
        wgs_meter2deg_model = partial(pyproj.transform,pyproj.Proj(init='EPSG:3857'),pyproj.Proj(init='EPSG:4326'))
        segment_geom_deg = transform(wgs_meter2deg_model, segment_geom)
        
        
        # write to the output shapefile
        output.write({
            'geometry': mapping(segment_geom_deg),
            'properties': data
        })
    
    
    print('You saved the shapefile successfully!')
segment_geom_deg.bounds<jupyter_output><empty_output><jupyter_text>Do the routing algorithm based on the spatio-temporal database<jupyter_code>outputShp<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09161426538761647, 'fraction_numerical': 0.031073575985679534, 'mean_word_length': 3.0288461538461537, 'pattern_counts': {'":': 0, '<': 52, '<?xml version=': 0, '>': 52, 'https://': 4, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '372535', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 7622, 'n_tokens_neox': 6862, 'n_words': 1753, 'path': 'libraries/Pnt2Lines.ipynb', 'repo_name': 'xiaojianggis/shadefinder', 'size': '21799'}