In everyday use, images are typically stored in 3-band 8-bit formats; for example a red pixel in a .png or .jpg file is denoted as (255, 0 , 0). Unsigned 16-bit imagery has a greater range of pixel values (from 0 to 65,535) and stores more information at the expense of greater storage space and reduced compatibility with analysis programs. While we encourage SpaceNet Dataset users to utilize the native 16-bit imagery, we find that visualizing results is often far easier with 8-bit imagery. Conversion is possible within packages such as NumPy, but in order to retain geographic data we utilize the gdal library. The code snippet below permits an end user to specify the percentile range. Specifically, this allows the removal of some outliers from 16-bit imagery, and preserves more relevant data when rescaling to the 8-bit range. The code can be found in apls_tools.py on the CosmiQ github repository.

def convert_to_8Bit(inputRaster, outputRaster,

outputPixType='Byte',

outputFormat='GTiff',

rescale_type='rescale',

percentiles=[2, 98]):

'''

Convert 16bit image to 8bit

rescale_type = [clip, rescale]

if clip, scaling is done strictly between 0 65535

if rescale, each band is rescaled to a min and max

set by percentiles

''' srcRaster = gdal.Open(inputRaster)

cmd = ['gdal_translate', '-ot', outputPixType, '-of',

outputFormat]



# iterate through bands

for bandId in range(srcRaster.RasterCount):

bandId = bandId+1

band = srcRaster.GetRasterBand(bandId)

if rescale_type == 'rescale':

bmin = band.GetMinimum()

bmax = band.GetMaximum()

# if not exist minimum and maximum values

if bmin is None or bmax is None:

(bmin, bmax) = band.ComputeRasterMinMax(1)

# else, rescale

band_arr_tmp = band.ReadAsArray()

bmin = np.percentile(band_arr_tmp.flatten(),

percentiles[0])

bmax= np.percentile(band_arr_tmp.flatten(),

percentiles[1]) else:

bmin, bmax = 0, 65535 cmd.append('-scale_{}'.format(bandId))

cmd.append('{}'.format(bmin))

cmd.append('{}'.format(bmax))

cmd.append('{}'.format(0))

cmd.append('{}'.format(255)) cmd.append(inputRaster)

cmd.append(outputRaster)

print "Conversin command:", cmd

subprocess.call(cmd)

3. Road Masks

The ultimate goal of the SpaceNet Roads Challenge is to extract a graph structure of road networks. Since the desired output is not the segmentation mask typically utilized for scoring (Post 1 provides a rationale for why we refrain from pixel-based metrics), any number of algorithmic approaches are possible. Nevertheless, one of the more obvious approaches to the challenge is to infer road masks from the 400 meter (1300 x 1300 pixel) imagery cutouts, and subsequently refine those inferred masks into a graph structure. In support of the segmentation approach, we create ground truth road masks for algorithm training.

Since the goal is to identify road centerlines, we refrain from attempting to mask the entire road width. Rather, we simply create a buffer (we use a buffer of 2 meters, yielding a total lane width of 4m) about the road centerline for use as a ground truth mask. Creating training masks is a two-step process.

Using GeoPandas we ingest the SpaceNet GeoJSON labels into a GeoDataFrame. Utilizing the geometry values of the data, we then create a buffer about the road centerline. The code snippet below is from apls_tools.py and allows the user to specify the desired buffer width (bufferDistanceMeters).

def create_buffer_geopandas(geoJsonFileName,

bufferDistanceMeters=2,

bufferRoundness=1, projectToUTM=True):

'''

Create a buffer around the lines of the geojson.

Return a geodataframe.

'''



inGDF = gpd.read_file(geoJsonFileName)



# set a few columns that we will need later

inGDF['type'] = inGDF['road_type'].values

inGDF['class'] = 'highway'

inGDF['highway'] = 'highway'



if len(inGDF) == 0:

return [], [] # Transform gdf Roadlines into UTM so that Buffer makes sense

if projectToUTM:

tmpGDF = ox.project_gdf(inGDF)

else:

tmpGDF = inGDF gdf_utm_buffer = tmpGDF # perform Buffer to produce polygons from Line Segments

gdf_utm_buffer['geometry'] =

tmpGDF.buffer(bufferDistanceMeters,

bufferRoundness) gdf_utm_dissolve = gdf_utm_buffer.dissolve(by='class')

gdf_utm_dissolve.crs = gdf_utm_buffer.crs if projectToUTM:

gdf_buffer = gdf_utm_dissolve.to_crs(inGDF.crs)

else:

gdf_buffer = gdf_utm_dissolve return gdf_buffer

2. We leverage the gdal library to convert the GeoDataFrame to a NumPy array, saving the array as an image. This process also ensures that the newly created mask image is snapped to the corresponding RGB image, and that the analogous pixel locations, rows, columns, and resolutions are identical. See apls_tools.py for code.

def gdf_to_array(gdf, im_file, output_raster, burnValue=150):



'''

Turn geodataframe to array, save as image file with non-null

pixels set to burnValue

''' NoData_value = 0 gdata = gdal.Open(im_file)



# set target info

target_ds =

gdal.GetDriverByName('GTiff').Create(output_raster,

gdata.RasterXSize, gdata.RasterYSize,

1, gdal.GDT_Byte)

target_ds.SetGeoTransform(gdata.GetGeoTransform())



# set raster info

raster_srs = osr.SpatialReference()

raster_srs.ImportFromWkt(gdata.GetProjectionRef())

target_ds.SetProjection(raster_srs.ExportToWkt())



band = target_ds.GetRasterBand(1)

band.SetNoDataValue(NoData_value)



outdriver=ogr.GetDriverByName('MEMORY')

outDataSource=outdriver.CreateDataSource('memData')

tmp=outdriver.Open('memData',1)

outLayer = outDataSource.CreateLayer("states_extent",

raster_srs, geom_type=ogr.wkbMultiPolygon)

# burn

burnField = "burn"

idField = ogr.FieldDefn(burnField, ogr.OFTInteger)

outLayer.CreateField(idField)

featureDefn = outLayer.GetLayerDefn()

for geomShape in gdf['geometry'].values:

outFeature = ogr.Feature(featureDefn)

outFeature.SetGeometry(ogr.CreateGeometryFromWkt(

geomShape.wkt))

outFeature.SetField(burnField, burnValue)

outLayer.CreateFeature(outFeature)

outFeature = 0



gdal.RasterizeLayer(target_ds, [1], outLayer,

burn_values=[burnValue])

The create_spacenet_masks.py script executes all the code included above, yielding training masks and illustrations of the workflow.

github.com/CosmiQ/apls/blob/master/src/create_spacenet_masks.py ... # iterate through images, convert to 8-bit, and create masks

im_files = os.listdir(path_images_raw)

for im_file in im_files:



name_root = im_file.split('_')[-1].split('.')[0]



# create 8-bit image

im_file_raw = os.path.join(path_images_raw, im_file)

im_file_out = os.path.join(path_images_8bit, im_file)

# convert to 8bit

apls_tools.convert_to_8Bit(im_file_raw, im_file_out,

outputPixType='Byte',

outputFormat='GTiff',

rescale_type='rescale',

percentiles=[2,98]) # determine output files

label_file = os.path.join(path_labels,

'spacenetroads_AOI_2_Vegas_' \

+ name_root + '.geojson')

label_file_tot = os.path.join(path_labels, label_file)

output_raster = os.path.join(path_masks, 'mask_' \

+ name_root + '.png')

plot_file = os.path.join(path_masks_plot, 'mask_' \

+ name_root + '.png') # create masks

mask, gdf_buffer = apls_tools.get_road_buffer(

label_file_tot, im_file_out,

output_raster,

buffer_meters=args.buffer_meters,

burnValue=args.burnValue,

bufferRoundness=6,

plot_file=plot_file,

figsize= (6,6),

fontsize=8, dpi=200,

show_plot=False, verbose=False)

The images below are examples of the output of the create_spacenet_masks.py script. We utilize the converted 8-bit images from Section 2 for visualization purposes.