# Read in the .env file
from dotenv import dotenv_values, load_dotenv

# load env variables
config = dotenv_values(".env")
load_dotenv()

# set env variables
GCS_PROJECT = config['GCS_PROJECT']
BUCKET = config['BUCKET']
HOME_PATH = config['HOME_PATH']
LOCAL_PATH = config['LOCAL_PATH']
DATA_PATH = config['DATA_PATH']
FOLDER = config['FOLDER']
GCLOUD_PATH = config['GCLOUD_PATH']
CUDA_VISIBLE_DEVICES = config['CUDA_VISIBLE_DEVICES']

print(config)

OrderedDict([('GCS_PROJECT', 'master-thesis-ilg'), ('BUCKET', 'csci5922-proj'), ('HOME_PATH', '/home/isly9493/'), ('LOCAL_PATH', '/local/isly9493/'), ('DATA_PATH', 'csci5922/proj/data/'), ('FOLDER', 'viirs-africa-many/'), ('GCLOUD_PATH', '/home/isly9493/google-cloud-sdk/bin/'), ('CUDA_VISIBLE_DEVICES', '1')])

# Import, authenticate and initialize the Earth Engine library.
import ee
ee.Authenticate()
ee.Initialize()

# Tensorflow setup.
import tensorflow as tf
print(tf.__version__)

# Folium setup.
import folium
print(folium.__version__)

# supress warnings
import warnings
warnings.filterwarnings('ignore')

2.15.0
0.15.1

# variables for image names
TRAINING_BASE = 'training_patches'
EVAL_BASE = 'eval_patches'

# Specify inputs (Landsat bands) to the model and the response variable.
opticalBands = ['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7']
thermalBands = ['B10', 'B11']
BANDS = opticalBands + thermalBands
RESPONSE = 'avg_rad' ##### 
FEATURES = BANDS + [RESPONSE]

# Specify the size and shape of patches expected by the model.
KERNEL_SIZE = 256
KERNEL_SHAPE = [KERNEL_SIZE, KERNEL_SIZE]
COLUMNS = [
  tf.io.FixedLenFeature(shape=KERNEL_SHAPE, dtype=tf.float32) for k in FEATURES
]
FEATURES_DICT = dict(zip(FEATURES, COLUMNS))

# Sizes of the training and evaluation datasets.
TRAIN_SIZE = 4800
EVAL_SIZE = 200

# Specify model training parameters.
BATCH_SIZE = 48 
EPOCHS = 1
BUFFER_SIZE = 2000
OPTIMIZER = 'SGD'
LOSS = 'MeanSquaredError'
METRICS = ['RootMeanSquaredError']

# Use Landsat 8 surface reflectance data.
l8sr = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')

# Cloud masking function.
def maskL8sr(image):
  cloudShadowBitMask = ee.Number(2).pow(3).int()
  cloudsBitMask = ee.Number(2).pow(5).int()
  qa = image.select('pixel_qa')
  mask1 = qa.bitwiseAnd(cloudShadowBitMask).eq(0).And(
    qa.bitwiseAnd(cloudsBitMask).eq(0))
  mask2 = image.mask().reduce('min')
  mask3 = image.select(opticalBands).gt(0).And(
          image.select(opticalBands).lt(10000)).reduce('min')
  mask = mask1.And(mask2).And(mask3)
  return image.select(opticalBands).divide(10000).addBands(
          image.select(thermalBands).divide(10).clamp(273.15, 373.15)
            .subtract(273.15).divide(100)).updateMask(mask)

# The image input data is a cloud-masked median composite.
image = l8sr.filterDate('2020-01-01', '2022-12-31').map(maskL8sr).median()

# printn out the bands
print(image.bandNames().getInfo())

# Use folium to visualize the imagery.
mapid = image.getMapId({'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.3})
map = folium.Map(location=[-1, 37], zoom_start=6)
folium.TileLayer(
    tiles=mapid['tile_fetcher'].url_format,
    attr='Map Data &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
    overlay=True,
    name='median composite',
  ).add_to(map)

mapid = image.getMapId({'bands': ['B10'], 'min': 0, 'max': 0.5})
folium.TileLayer(
    tiles=mapid['tile_fetcher'].url_format,
    attr='Map Data &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
    overlay=True,
    name='thermal',
  ).add_to(map)
map.add_child(folium.LayerControl())
map

['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'B11']

viirs = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG').select('avg_rad').filterDate('2020-01-01', '2022-12-31').median()
viirs = viirs.divide(60).float() # normalize to 0-1 range

mapid = viirs.getMapId({'min': 0, 'max': 0.02}) 
map = folium.Map(location=[-1, 37], zoom_start=6)
folium.TileLayer(
    tiles=mapid['tile_fetcher'].url_format,
    attr='Map Data &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
    overlay=True,
    name='viirs dnb',
  ).add_to(map)
map.add_child(folium.LayerControl())
map

featureStack = ee.Image.cat([
  image.select(BANDS),
  viirs.select(RESPONSE)
]).float()

list = ee.List.repeat(1, KERNEL_SIZE)
lists = ee.List.repeat(list, KERNEL_SIZE)
kernel = ee.Kernel.fixed(KERNEL_SIZE, KERNEL_SIZE, lists)

arrays = featureStack.neighborhoodToArray(kernel)

trainingPolys = ee.FeatureCollection('projects/'+GCS_PROJECT+'/assets/trainingPolysAfricaBig')
evalPolys = ee.FeatureCollection('projects/'+GCS_PROJECT+'/assets/evalPolysAfrica')

polyImage = ee.Image(0).byte().paint(trainingPolys, 1).paint(evalPolys, 2)
polyImage = polyImage.updateMask(polyImage)

mapid = polyImage.getMapId({'min': 1, 'max': 2, 'palette': ['gold', 'blue']})
map = folium.Map(location=[0, 12], zoom_start=5)
folium.TileLayer(
    tiles=mapid['tile_fetcher'].url_format,
    attr='Map Data &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
    overlay=True,
    name='training polygons',
  ).add_to(map)
map.add_child(folium.LayerControl())
map

# Convert the feature collections to lists for iteration.
trainingPolysList = trainingPolys.toList(trainingPolys.size())
evalPolysList = evalPolys.toList(evalPolys.size())

# These numbers determined experimentally.
n = 20 # Number of shards in each polygon. 
N = 200 # Total sample size in each polygon.

# Export all the training data (in many pieces), with one task
# per geometry.
for g in range(trainingPolys.size().getInfo()):
  geomSample = ee.FeatureCollection([])
  for i in range(n):
    sample = arrays.sample(
      region = ee.Feature(trainingPolysList.get(g)).geometry(),
      scale = 30,
      numPixels = N / n, # Size of the shard.
      seed = i,
      tileScale = 8
    )
    geomSample = geomSample.merge(sample)
  
  desc = TRAINING_BASE + '_g' + str(g)
  print(g, FOLDER, BUCKET, desc)
  
  task = ee.batch.Export.table.toCloudStorage(
    collection = geomSample,
    description = desc,
    bucket = BUCKET,
    fileNamePrefix = FOLDER + desc,
    fileFormat = 'TFRecord',
    selectors = BANDS + [RESPONSE]
  )
  task.start()
print('Training imagery tasks submitted. Check on the Tasks in GEE for their status (https://code.earthengine.google.com/tasks).')

# Export all the evaluation data.
for g in range(evalPolys.size().getInfo()):
  geomSample = ee.FeatureCollection([])
  for i in range(n):
    sample = arrays.sample(
      region = ee.Feature(evalPolysList.get(g)).geometry(),
      scale = 30,
      numPixels = N / n,
      seed = i,
      tileScale = 8
    )
    geomSample = geomSample.merge(sample)

  desc = EVAL_BASE + '_g' + str(g)
  print(g, FOLDER, BUCKET, desc)
  
  task = ee.batch.Export.table.toCloudStorage(
    collection = geomSample,
    description = desc,
    bucket = BUCKET,
    fileNamePrefix = FOLDER + desc,
    fileFormat = 'TFRecord',
    selectors = BANDS + [RESPONSE]
  )
  task.start()
print('Eval imagery tasks submitted. Check on the Tasks in GEE for their status (https://code.earthengine.google.com/tasks).')

# Copy from your bucket to local path (note -r is for recursive call)
!{GCLOUD_PATH}gsutil -m cp -r 'gs://'{BUCKET}'/'{FOLDER}{TRAINING_BASE}'*' {LOCAL_PATH}{DATA_PATH}{FOLDER}
!{GCLOUD_PATH}gsutil -m cp -r 'gs://'{BUCKET}'/'{FOLDER}{EVAL_BASE}'*' {LOCAL_PATH}{DATA_PATH}{FOLDER}

# see what devices are available
print("Num GPUs Available: ", len(tf.config.list_physical_devices('GPU')))
print("GPUs Available: ", tf.config.list_physical_devices('GPU'))
!nvidia-smi

Num GPUs Available:  1
GPUs Available:  [PhysicalDevice(name='/physical_device:GPU:0', device_type='GPU')]

2023-12-12 16:01:15.741987: I external/local_xla/xla/stream_executor/cuda/cuda_executor.cc:901] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero. See more at https://github.com/torvalds/linux/blob/v6.0/Documentation/ABI/testing/sysfs-bus-pci#L344-L355
2023-12-12 16:01:15.773299: I external/local_xla/xla/stream_executor/cuda/cuda_executor.cc:901] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero. See more at https://github.com/torvalds/linux/blob/v6.0/Documentation/ABI/testing/sysfs-bus-pci#L344-L355
2023-12-12 16:01:15.773568: I external/local_xla/xla/stream_executor/cuda/cuda_executor.cc:901] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero. See more at https://github.com/torvalds/linux/blob/v6.0/Documentation/ABI/testing/sysfs-bus-pci#L344-L355

def parse_tfrecord(example_proto):
  return tf.io.parse_single_example(example_proto, FEATURES_DICT)


def to_tuple(inputs):
  inputsList = [inputs.get(key) for key in FEATURES]
  stacked = tf.stack(inputsList, axis=0)
  # Convert from CHW to HWC
  stacked = tf.transpose(stacked, [1, 2, 0])
  return stacked[:,:,:len(BANDS)], stacked[:,:,len(BANDS):]


def get_dataset(path, prefix, n_files):
  glob = []
  for i in range(0, n_files):
    glob = path + prefix + '_g' + str(i) + '.tfrecord.gz'
  dataset = tf.data.TFRecordDataset(glob, compression_type='GZIP')
  dataset = dataset.map(parse_tfrecord, num_parallel_calls=5)
  dataset = dataset.map(to_tuple, num_parallel_calls=5)
  return dataset

training = get_dataset(LOCAL_PATH + DATA_PATH + FOLDER, TRAINING_BASE, 24)
training = training.shuffle(BUFFER_SIZE).batch(BATCH_SIZE).repeat()
# print(iter(training.take(1)).next())

evaluation = get_dataset(LOCAL_PATH + DATA_PATH + FOLDER, EVAL_BASE, 2)
evaluation = evaluation.batch(1).repeat()
# print(iter(evaluation.take(1)).next()[0][0,:,:,1:3])

# visualize of the of the training / evaluation image patch pairs
import numpy as np
import matplotlib.pyplot as plt

def display_batch_of_images(databatch):
    imgs, labels = databatch
    # landsat in false color
    plt.figure(figsize=(10,10))
    label = labels.numpy().astype(np.float32)
    for i in range(25):
        plt.subplot(5,5,i+1)
        plt.imshow(imgs[i,:,:,0:3]*8) # select landsat bands here
        plt.axis('off')
    plt.show()
    # show viirs in grayscale
    plt.figure(figsize=(10,10))
    for i in range(25):
        plt.subplot(5,5,i+1)
        plt.imshow(labels[i,:,:,0], cmap='gray')
        plt.axis('off')
    plt.show()

display_batch_of_images(iter(training.take(1)).next())

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

from tensorflow import keras
from keras import layers
from keras import losses
from keras import models
from keras import metrics
from keras import optimizers

def conv_block(input_tensor, num_filters):
	encoder = layers.Conv2D(num_filters, (3, 3), padding='same')(input_tensor)
	encoder = layers.BatchNormalization()(encoder)
	encoder = layers.Activation('relu')(encoder)
	encoder = layers.Conv2D(num_filters, (3, 3), padding='same')(encoder)
	encoder = layers.BatchNormalization()(encoder)
	encoder = layers.Activation('relu')(encoder)
	return encoder

def encoder_block(input_tensor, num_filters):
	encoder = conv_block(input_tensor, num_filters)
	encoder_pool = layers.MaxPooling2D((2, 2), strides=(2, 2))(encoder)
	return encoder_pool, encoder


def decoder_block(input_tensor, concat_tensor, num_filters):
	decoder = layers.Conv2DTranspose(num_filters, (2, 2), strides=(2, 2), padding='same')(input_tensor)
	decoder = layers.concatenate([concat_tensor, decoder], axis=-1)
	decoder = layers.BatchNormalization()(decoder)
	decoder = layers.Activation('relu')(decoder)
	decoder = layers.Conv2D(num_filters, (3, 3), padding='same')(decoder)
	decoder = layers.BatchNormalization()(decoder)
	decoder = layers.Activation('relu')(decoder)
	decoder = layers.Conv2D(num_filters, (3, 3), padding='same')(decoder)
	decoder = layers.BatchNormalization()(decoder)
	decoder = layers.Activation('relu')(decoder)
	return decoder

def get_model():
	inputs = layers.Input(shape=[None, None, len(BANDS)]) # 256
	encoder0_pool, encoder0 = encoder_block(inputs, 32) # 128
	encoder1_pool, encoder1 = encoder_block(encoder0_pool, 64) # 64
	encoder2_pool, encoder2 = encoder_block(encoder1_pool, 128) # 32
	encoder3_pool, encoder3 = encoder_block(encoder2_pool, 256) # 16
	encoder4_pool, encoder4 = encoder_block(encoder3_pool, 512) # 8
	center = conv_block(encoder4_pool, 1024) # center
	decoder4 = decoder_block(center, encoder4, 512) # 16
	decoder3 = decoder_block(decoder4, encoder3, 256) # 32
	decoder2 = decoder_block(decoder3, encoder2, 128) # 64
	decoder1 = decoder_block(decoder2, encoder1, 64) # 128
	decoder0 = decoder_block(decoder1, encoder0, 32) # 256

	# to change your output activation function, change this part
	outputs = layers.Conv2D(1, (1, 1), activation='sigmoid')(decoder0)

	model = models.Model(inputs=[inputs], outputs=[outputs])
	
	# set optimizer with custom learning rate
	OPTIMIZER = 'SGD'
	optimizer = optimizers.get(OPTIMIZER)
	
	# to adjust the learning rate:
	# optimizer.learning_rate = 0.1

	# compile the model
	model.compile(
		optimizer=optimizer,
		loss=losses.get(LOSS),
		metrics=[metrics.get(metric) for metric in METRICS])

	return model

# print the current time to log how long it takes to train
from datetime import datetime
print(datetime.now())

# get the model fresh and train it
m = get_model()
history = m.fit(
    x=training,
    epochs=EPOCHS,
    steps_per_epoch=int(TRAIN_SIZE / BATCH_SIZE),
    validation_data=evaluation,
    validation_steps=EVAL_SIZE
)

# print the end time to log how long it takes to train
print(datetime.now())

2023-12-12 16:18:50.358788
100/100 [==============================] - 42s 361ms/step - loss: 0.2118 - root_mean_squared_error: 0.4602 - val_loss: 0.1753 - val_root_mean_squared_error: 0.4187
2023-12-12 16:19:33.359108

# set some parameters for saving the model
MODEL_NAME = 'viirs_africa_many_e1.keras'
HISTORY_NAME = 'viirs_africa_many_e1_history.pickle'
ASSET_SUFFIX = 'viirs_africa_many_e1'

# save the model to a local file 
m.save(LOCAL_PATH + DATA_PATH + FOLDER + MODEL_NAME)

# # save the history variable to a local file with pickle
import pickle
with open(LOCAL_PATH + DATA_PATH + FOLDER + HISTORY_NAME, 'wb') as f:
  pickle.dump(history, f)

# import a saved model if needed
# m = tf.keras.models.load_model(LOCAL_PATH + DATA_PATH + FOLDER + MODEL_NAME)
# import pickle
# with open(HISTORY_NAME, 'rb') as f:
# with open(LOCAL_PATH + DATA_PATH + FOLDER + HISTORY_NAME, 'rb') as f:
#   history = pickle.load(f)

# plot the training and validation loss
import matplotlib.pyplot as plt

plt.plot(history.history['root_mean_squared_error'])
plt.plot(history.history['val_root_mean_squared_error'])
plt.title('Model RMSE for ' + ASSET_SUFFIX + ' for ' + str(EPOCHS) + ' epochs')
plt.ylabel('RMSE')
plt.xlabel('Epoch')
plt.legend(['Training RMSE', 'Validation RMSE'], loc='upper left')
plt.show()

## evaluate the model on the evaluation dataset to get the loss and RMSE
val_loss, val_rmse = m.evaluate(evaluation, steps=EVAL_SIZE)
print('Model validation loss: {:0.4f}. RMSE: {:0.4f}'.format(val_loss, val_rmse))

200/200 [==============================] - 1s 7ms/step - loss: 0.1753 - root_mean_squared_error: 0.4187
Model validation loss: 0.1753. RMSE: 0.4187

## evaluate the model on the evaluation dataset to get the loss and RMSE
train_loss, train_rmse = m.evaluate(training, steps=10)
print('Model train loss: {:0.4f}. RMSE: {:0.4f}'.format(train_loss, train_rmse))

10/10 [==============================] - 3s 212ms/step - loss: 0.2333 - root_mean_squared_error: 0.4830
Model train loss: 0.2333. RMSE: 0.4830

ee_user_folder = 'projects/'+GCS_PROJECT+'/assets'

pred_image_base = 'FCNN_demo_kinshasa_384_'
# Half this will extend on the sides of each patch.
pred_kernel_buffer = [128, 128]
# Kinshasa/Brazzaville
pred_region = ee.Geometry.Polygon(
  [[[15.02, -4.04],
    [15.02, -4.54],
    [15.52, -4.54],
    [15.52, -4.04]]], 
  None, 
  False
)

def doExport(out_image_base, kernel_buffer, region):
  task = ee.batch.Export.image.toCloudStorage(
    image = image.select(BANDS),
    description = out_image_base,
    bucket = BUCKET,
    fileNamePrefix = FOLDER + out_image_base,
    region = region.getInfo()['coordinates'],
    scale = 30,
    fileFormat = 'TFRecord',
    maxPixels = 1e10,
    formatOptions = {
      'patchDimensions': KERNEL_SHAPE,
      'kernelSize': kernel_buffer,
      'compressed': True,
      'maxFileSize': 104857600
    }
  )
  task.start()
  print('Prediction image export task submitted. Check on the task in GEE for its status (https://code.earthengine.google.com/tasks).')

# Run the export.
doExport(pred_image_base, pred_kernel_buffer, pred_region)

Prediction image export task submitted. Check on the task in GEE for its status (https://code.earthengine.google.com/tasks).

# copy the prediction file from cloud storage to local
!{GCLOUD_PATH}gsutil -m cp -r 'gs://'{BUCKET}'/'{FOLDER}{pred_image_base}'*' {LOCAL_PATH}{DATA_PATH}{FOLDER}

def doPrediction(out_image_base, kernel_buffer):
  print('Looking for TFRecord files...')

  # # Get a list of all the files in the output bucket.
  filesList = !ls {LOCAL_PATH}{DATA_PATH}{FOLDER}

  # Get only the files generated by the image export.
  exportFilesList = [s for s in filesList if out_image_base in s]

  # Get the list of image files and the JSON mixer file.
  imageFilesList = []
  jsonFile = None
  for f in exportFilesList:
    if f.endswith('.tfrecord.gz'):
      imageFilesList.append(f)
    elif f.endswith('.json'):
      jsonFile = f
  imageFilesList.sort()

  from pprint import pprint
  pprint(imageFilesList)
  print(jsonFile)

  import json
  jsonText = !cat {LOCAL_PATH}{DATA_PATH}{FOLDER}{jsonFile}
  print(jsonText)
  mixer = json.loads(jsonText.nlstr)
  pprint(mixer)
  patches = mixer['totalPatches']

  buffered_shape = [
      KERNEL_SHAPE[0] + kernel_buffer[0],
      KERNEL_SHAPE[1] + kernel_buffer[1]]

  imageColumns = [
    tf.io.FixedLenFeature(shape=buffered_shape, dtype=tf.float32)
      for k in BANDS
  ]

  imageFeaturesDict = dict(zip(BANDS, imageColumns))

  def parse_image(example_proto):
    return tf.io.parse_single_example(example_proto, imageFeaturesDict)

  def toTupleImage(inputs):
    inputsList = [inputs.get(key) for key in BANDS]
    stacked = tf.stack(inputsList, axis=0)
    stacked = tf.transpose(stacked, [1, 2, 0])
    return stacked

  # append absolute path to imageFilesList
  imageFilesList = [LOCAL_PATH + DATA_PATH + FOLDER + f for f in imageFilesList]
  print(imageFilesList)
  
  # Create a dataset from the TFRecord file(s) in Cloud Storage.
  imageDataset = tf.data.TFRecordDataset(imageFilesList, compression_type='GZIP')
  imageDataset = imageDataset.map(parse_image, num_parallel_calls=5)
  imageDataset = imageDataset.map(toTupleImage).batch(1)

  # Perform inference.
  print('Running predictions...')
  predictions = m.predict(imageDataset, steps=patches, verbose=1)
  # print(predictions[0])
  return predictions

def writePredictions(predictions, out_image_base, kernel_buffer):
  print('Writing predictions...')
  out_image_file = 'gs://' + BUCKET + '/' + FOLDER + out_image_base + '.TFRecord'

  # Get set up for prediction.
  x_buffer = int(kernel_buffer[0] / 2)
  y_buffer = int(kernel_buffer[1] / 2)
  
  
  writer = tf.io.TFRecordWriter(out_image_file)
  patches = 0
  for predictionPatch in predictions:
    # print('Writing patch ' + str(patches) + '...')
    predictionPatch = predictionPatch[
        x_buffer:x_buffer+KERNEL_SIZE, y_buffer:y_buffer+KERNEL_SIZE]

    # Create an example.
    example = tf.train.Example(
      features=tf.train.Features(
        feature={
          'impervious': tf.train.Feature(
            float_list=tf.train.FloatList(
              value=predictionPatch.flatten()
            )
          )
        }
      )
    )
    # Write the example.
    writer.write(example.SerializeToString())
    patches += 1

  print('Done writing predictions.')
  writer.close()

def uploadPredictionsToEE(out_image_base, ee_user_folder):
  # Start the upload.
  print('Uploading to Google Earth Egnine...')
  
  # define file paths
  out_image_file = 'gs://' + BUCKET + '/' + FOLDER + out_image_base + '.TFRecord'
  out_image_asset = ee_user_folder + '/' + out_image_base + ASSET_SUFFIX
  jsonFile = 'gs://' + BUCKET + '/' + FOLDER + out_image_base + 'mixer.json'
  
  print('out_image_asset: ', out_image_asset)
  print('out_image_file: ', out_image_file)
  print('jsonFile: ', jsonFile)
  
  # upload file to google earth engine from gcs
  !earthengine upload image --asset_id={out_image_asset} {out_image_file} {jsonFile}

  print('Prediction image upload task submitted. Check on the task in GEE for its status (https://code.earthengine.google.com/tasks).')

# Run the prediction.
predictions = doPrediction(pred_image_base, pred_kernel_buffer)

# Write the predictions.
writePredictions(predictions, pred_image_base, pred_kernel_buffer)

# upload to earth engine
uploadPredictionsToEE(pred_image_base, ee_user_folder)

map = folium.Map(location=[-4.29, 15.27], zoom_start=11)

# add landsat layer clipped to prediction region
landsat = image
landsat = landsat.clip(pred_region)
mapid = landsat.getMapId({'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.3})
folium.TileLayer(
    tiles=mapid['tile_fetcher'].url_format,
    attr='Map Data &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
    overlay=True,
    name='landsat',
  ).add_to(map)

# actual viirs
viirs = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG').select('avg_rad').filterDate('2021-01-01', '2022-12-31').median().divide(30).float()
viirs = viirs.clip(pred_region)
# actual_image = viirs
mapid = viirs.getMapId({'min': 0, 'max': 1}) # normally max 1.0
folium.TileLayer(
    tiles=mapid['tile_fetcher'].url_format,
    attr='Map Data &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
    overlay=True,
    name='viirs dnb',
  ).add_to(map)

# prediction
out_image = ee.Image(ee_user_folder + '/' + pred_image_base + ASSET_SUFFIX)
mapid = out_image.getMapId({'min': 0, 'max': 1})
folium.TileLayer(
    tiles=mapid['tile_fetcher'].url_format,
    attr='Map Data &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
    overlay=True,
    name='predicted nightlights',
  ).add_to(map)


# difference
diff = out_image.subtract(viirs)
mapid = diff.getMapId({'min': -0.5, 'max': 0.5, 'palette': ['FF0000', '000000', '00FF00']})
folium.TileLayer(
    tiles=mapid['tile_fetcher'].url_format,
    attr='Map Data &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
    overlay=True,
    name='difference',
  ).add_to(map)

# visualize the map
map.add_child(folium.LayerControl())
map

# export to HTML for webpage
import os

# os.system('jupyter nbconvert --to html mod1.ipynb')
os.system('jupyter nbconvert --to html unet_regression_virrs.ipynb --HTMLExporter.theme=dark')

Earth at Night in HD: Predicting High-Resolutions Nighttime Radiance from Daytime Satellite Imagery¶

Introduction¶

Methods¶

Code for Implementing the Model¶

Configure the Environment¶

Import Key Libraries and Authenticate to Google Earth Engine¶

Set other Global Variables¶

Data Preparation¶

Define Training and Evaluation Areas¶

Sampling¶

Export Training & Evaluation Data from Earth Engine to Google Cloud Storage¶

Copy Files from Google Cloud Storage to Local Storage¶

Configure GPU Usage¶

Load Training and Evaluation Data in Memory for Training.¶

Evaluation data¶

Data Analysis & Deep Learning Model¶

Training the model¶

Evaluation¶

Prediction¶

Specify Region and Outputs for Prediction¶

Export Sample Imagery for Prediction¶

Define Code for Running Model on Prediction Data¶

Run Predictions¶

Display the Prediction¶

Results¶

Conclusions¶

References¶

Appendix A: Export to HTML¶

Appendix B: To Do List¶