#################################################################################
## https://www.kaggle.com/c/dstl-satellite-imagery-feature-detection
## Dstl Satellite Imagery Feature Detection
## Santiago Mota
## santiago_mota@yahoo.es
## https://es.linkedin.com/in/santiagomota/en

# Can you train an eye in the sky?

# The proliferation of satellite imagery has given us a radically improved 
# understanding of our planet. It has enabled us to better achieve everything 
# from mobilizing resources during disasters to monitoring effects of global 
# warming. What is often taken for granted is that advancements such as these 
# have relied on labeling features of significance like building footprints and 
# roadways fully by hand or through imperfect semi-automated methods.

# As these large, complex datasets continue to increase exponentially in number, 
# the Defence Science and Technology Laboratory (Dstl) is seeking novel solutions 
# to alleviate the burden on their image analysts. In this competition, Kagglers 
# are challenged to accurately classify features in overhead imagery. Automating 
# feature labeling will not only help Dstl make smart decisions more quickly 
# around the defense and security of the UK, but also bring innovation to 
# computer vision methodologies applied to satellite imagery.

# Started: 12:00 am, Thursday 15 December 2016 UTC
# Ends: 11:59 pm, Tuesday 7 March 2017 UTC (82 total days)
# Points: this competition awards standard ranking points
# Tiers: this competition counts towards tiers

# Data Files
# File Name 	      Available Formats
# train_wkt.csv 	      .zip (9.82 mb)
# sample_submission.csv .zip (14.89 kb)
# train_geojson 	      .zip (10.78 mb)
# grid_sizes.csv 	      .zip (2.17 kb)
# sixteen_band 	      .zip (7.30 gb)
# three_band 	      .zip (12.87 gb)

# In this competition, Dstl provides you with 1km x 1km satellite images in both 
# 3-band and 16-band formats. Your goal is to detect and classify the types of 
# objects found in these regions. 

# 3- and 16-bands images

# There are two types of imagery spectral content provided in this competition. 
# The 3-band images are the traditional RGB natural color images. The 16-band 
# images contain spectral information by capturing wider wavelength channels. 
# This multi-band imagery is taken from the multispectral (400 – 1040nm) and 
# short-wave infrared (SWIR) (1195-2365nm) range. All images are in GeoTiff 
# format and might require GeoTiff viewers (such as QGIS) to view. Please refer 
# to our tutorial on how to programmatically view the images.

# Ortho Ready Standard Imagery © 2016 DigitalGlobe, Inc.

# All imagery credit to: Satellite Imagery © DigitalGlobe, Inc.

# Imagery details

# Sensor : WorldView 3
# Wavebands :
#  - Panchromatic: 450-800 nm
#  - 8 Multispectral: (red, red edge, coastal, blue, green, yellow, near-IR1 and 
#    near-IR2) 400 nm - 1040 nm
#  - 8 SWIR: 1195 nm - 2365 nm
# Sensor Resolution (GSD) at Nadir :
#  - Panchromatic: 0.31m 
#  - Multispectral: 1.24 m
#  - SWIR: Delivered at 7.5m
# Dynamic Range
#  - Panchromatic and multispectral : 11-bits per pixel
#  - SWIR : 14-bits per pixel

# Object types

# In a satellite image, you will find lots of different objects like roads, 
# buildings, vehicles, farms, trees, water ways, etc. Dstl has labeled 10 
# different classes:
      
#  1. Buildings - large building, residential, non-residential, fuel storage 
#                 facility, fortified building
#  2. Misc. Manmade structures 
#  3. Road 
#  4. Track - poor/dirt/cart track, footpath/trail
#  5. Trees - woodland, hedgerows, groups of trees, standalone trees
#  6. Crops - contour ploughing/cropland, grain (wheat) crops, row (potatoes, 
#             turnips) crops
#  7. Waterway 
#  8. Standing water
#  9. Vehicle Large - large vehicle (e.g. lorry, truck,bus), logistics vehicle
# 10. Vehicle Small - small vehicle (car, van), motorbike

# Every object class is described in the form of Polygons and MultiPolygons, 
# which are simply a list of polygons. We provide two different formats for 
# these shapes: GeoJson and WKT. These are both open source formats for 
# geo-spatial shapes. 

# Your submission will be in the WKT format. 

# Geo Coordinates

# In this dataset that we provide, we create a set of geo-coordinates that are 
# in the range of x = [0,1] and y = [-1,0]. These coordinates are transformed 
# such that we obscure the location of where the satellite images are taken 
# from. The images are from the same region on Earth.

# To utilize these images, we provide the grid coordinates of each image so you 
# know how to scale them and align them with the images in pixels. You need the 
# Xmax and Ymin for each image to do the scaling (provided in our 
# grid_sizes.csv) Please refer to our tutorial on how to programmatically view 
# the images.

# File descriptions

# train_wkt.csv - the WKT format of all the training labels
#  - ImageId - ID of the image
#  - ClassType - type of objects (1-10)
#  - MultipolygonWKT - the labeled area, which is multipolygon geometry 
#                      represented in WKT format 
# three_band.zip - the complete dataset of 3-band satellite images. The three 
#                  bands are in the images with file name = {ImageId}.tif. MD5 
#                  = 7cf7bf17ba3fa3198a401ef67f4ef9b4 
# sixteen_band.zip - the complete dataset of 16-band satellite images. The 16 
#                    bands are distributed in the images with file name = 
#                    {ImageId}_{A/M/P}.tif. MD5 = e2949f19a0d1102827fce35117c5f08a
# grid_sizes.csv - the sizes of grids for all the images
#  - ImageId - ID of the image
#  - Xmax - maximum X coordinate for the image
#  - Ymin - minimum Y coordinate for the image
# sample_submission.csv - a sample submission file in the correct format
# ImageId - ID of the image
# ClassType - type of objects (1-10)
# MultipolygonWKT - the labeled area, which is multipolygon geometry represented 
#                   in WKT format
# train_geojson.zip - the geojson format of all the training labels (essentially 
#                     these are the same information as train_wkt.csv) 

# filename_to_classType = {
#       '001_MM_L2_LARGE_BUILDING':1,
#       '001_MM_L3_RESIDENTIAL_BUILDING':1,
#       '001_MM_L3_NON_RESIDENTIAL_BUILDING':1,
#       '001_MM_L5_MISC_SMALL_STRUCTURE':2,
#       '002_TR_L3_GOOD_ROADS':3,
#       '002_TR_L4_POOR_DIRT_CART_TRACK':4,
#       '002_TR_L6_FOOTPATH_TRAIL':4,
#       '006_VEG_L2_WOODLAND':5,
#       '006_VEG_L3_HEDGEROWS':5,
#       '006_VEG_L5_GROUP_TREES':5,
#       '006_VEG_L5_STANDALONE_TREES':5,
#       '007_AGR_L2_CONTOUR_PLOUGHING_CROPLAND':6,
#       '007_AGR_L6_ROW_CROP':6, 
#       '008_WTR_L3_WATERWAY':7,
#       '008_WTR_L2_STANDING_WATER':8,
#       '003_VH_L4_LARGE_VEHICLE':9,
#       '003_VH_L5_SMALL_VEHICLE':10,
#       '003_VH_L6_MOTORBIKE':10}

################################################################################
# Some kernel information
dir('../input')
# list.dirs('../input/train_geojson')
system("df")
system("lsblk -a")
cat(system("lsblk", intern = TRUE), sep = '\n')

# Load some packages
# devtools::install_github("ropensci/geojsonio")
library("geojsonio")

# install.packages("rgdal", type = "source")
# install.packages("rgeos", type = "source")
library('rgdal')
library('rgeos')
library('ggplot2')
library('raster')
library('readr')
library('plyr')

##### Grid sizes #####

grid_sizes <- read_csv('../input/grid_sizes.csv')
summary(grid_sizes)
class(grid_sizes)

# Number of images 450

head(grid_sizes)

grid_sizes <- as.data.frame(grid_sizes)

table(grid_sizes$Xmax)

table(grid_sizes$Ymin)

##### Train wkt #####
# New wkt file
train_wkt <- read_csv('../input/train_wkt_v4.csv')

summary(train_wkt)
class(train_wkt)

head(train_wkt)

##### Three_band #####

raster_6010_1_2 <- raster("../input/three_band/6010_1_2.tif")
plot(raster_6010_1_2, main='Three Band 6010_1_2')

# Copied from https://www.kaggle.com/jeffhebert/dstl-satellite-imagery-feature-detection/notebook1b5997d25c/discussion
# Thanks @JeffH
img_6010_1_2 <- stack("../input/three_band/6010_1_2.tif")
plotRGB(img_6010_1_2, stretch = "lin", main='Three Band 6010_1_2')

gdal_6010_1_2 <- readGDAL(paste0("../input/three_band/", '6010_1_2', ".tif"))
plot(gdal_6010_1_2, main='Three Band 6010_1_2')

##### Submission #####
submission <- read_csv('../input/sample_submission.csv')

# Number of train (21) and test (429) 429 + 21 = 450 
test_id  <- unique(submission$ImageId)
train_id <- unique(train_wkt$ImageId)

print(train_id)
print(test_id)

# print('Train geojson')

##### Train geojson #####
# One train geojson
grid_6010 <- geojson_read("../input/train_geojson_v3/6010_1_2/Grid_6010.geojson", method = local, what= 'sp')
                           
class(grid_6010)

# library('maps')
# plot(grid_6010)

# ggplot(grid_6010, aes(long, lat, group = group)) +
#         geom_polygon()

# Load geojson
grid_6010_1_2 <- geojson_read("../input/train_geojson_v3/6010_1_2/Grid_6010.geojson",
                               method = local, what= 'sp')
class(grid_6010_1_2)

# plot(grid_6010_1_2)

# ggplot(grid_6010_1_2, aes(long, lat, group = group)) +
#         geom_polygon()

##### File system #####
print('File system')
dir('../input')
# [1] "grid_sizes.csv"  "sample_submission.csv" "sixteen_band" "three_band"           
# [5] "train_geojson"   "train_wkt.csv" 

directories         <- list.dirs('../input')                  # 27
print('Directories')
print(directories)

directories_geojson <- list.dirs('../input/train_geojson_v3') # 21
print('Directories geojson')
print(directories_geojson)

# theree_band_files   <- dir('../input/three_band')    # 450
# theree_band_files   <- list.dirs('../input/three_band')    # 450
# print(three_band_files)

# sixteen_band_files  <- dir('../input/sixteen_band')  # 1350 (450 * 3)
# print(sixteen_band_files)

geojson_files       <- lapply(directories_geojson, dir)

##### Locate train files #####
train_id[!(train_id %in% geojson_files[[1]])]
# [1] "6010_4_2" "6140_3_1" "6110_4_0" "6150_2_3"

geojson_files[[1]][!(geojson_files[[1]] %in% train_id)]
# [1] "6010_1_2" "6040_4_4" "6100_2_2"

manual <- c('6100_1_3', '6010_4_2', '6010_4_4', '6140_3_1', '6170_2_4',
            '6040_1_3', '6040_2_2', '6170_4_1', '6110_4_0', '6120_2_2',
            '6100_2_3', '6120_2_0', '6150_2_3', '6110_1_2', '6170_0_4',
            '6160_2_1', '6090_2_0', '6140_1_2', '6060_2_3', '6110_3_1',
            '6040_1_0')

train_id[!(train_id %in% manual)]
# character(0)

manual[!(manual %in% train_id)]
# character(0)

# Check out some problem here
geojson_files[[1]][!(geojson_files[[1]] %in% manual)]
# [1] "6010_1_2" "6040_4_4" "6100_2_2"

##### Plot some geojson #####

# First directory
base <- '../input/train_geojson_v3/6010_1_2/'
# st_files <- paste0(base, st_names)
# st_names <- Filter(function(x) grepl("\\.geojson", x), sapply(content(res), "[[", "name"))
# st_names <- Filter(function(x) grepl("\\.geojson", x), dir(base))
st_names <- dir(base)
length(st_names)
# [1] 8
st_files <- paste0(base, geojson_files[[2]])
st_files

st_use <- st_files
geo_6010_1_2 <- lapply(st_use, geojson_read, method = "local", what = "sp")
df_6010_1_2  <- ldply(setNames(lapply(geo_6010_1_2, fortify), 
                               gsub("\\.geojson", "", st_names)))
ggplot(df_6010_1_2, aes(long, lat, group = group)) +
      geom_polygon() +
      facet_wrap(~.id, scales = "free")