Land Surface Temperature Data

Report date: May 29, 2024
BACK TO SURFACE HEAT PROJECT DESCRIPTION.

Project summary

Get the right data from Google Earth Engine and understand it.
Bring it into QGIS to look at it; make sure the temperature patterns it illustrates make sense.
Attach the proper color scheme.
Output as MB Tiles for inclusion in the maplibre display map.

1. Get and inspect the data.

Here is the code (written in Javascript) to acquire the data from Google Earth Engine.

// Define the geography of interest
var bexarCounty = ee.FeatureCollection('TIGER/2016/Counties')
                    .filter(ee.Filter.eq('NAME', 'Bexar'));

// Function to load and preprocess data for a given year
function processYear(year) {
  var startDate = ee.Date.fromYMD(year, 5, 1);
  var endDate = ee.Date.fromYMD(year, 9, 30);

  var dataset = ee.ImageCollection('MODIS/061/MOD11A1')
                    .filterDate(startDate, endDate)
                    .filterBounds(bexarCounty)
                    .select('LST_Day_1km');

  var meanLST = dataset.mean().multiply(0.02);
  var LSTCelsius = meanLST.subtract(273.15);
  var LSTFahrenheit = LSTCelsius.multiply(1.8).add(32);

  var reductionParams = {
    reducer: ee.Reducer.minMax(),
    geometry: bexarCounty,
    scale: 1000,
    maxPixels: 1e9
  };
  var minMax = LSTFahrenheit.reduceRegion(reductionParams);
  minMax.evaluate(function(result) {
    console.log('Year ' + year + ' - Min: ' + result['LST_Day_1km_min'] + ', Max: ' + result['LST_Day_1km_max']);
  });

  Export.image.toDrive({
    image: LSTFahrenheit,
    description: 'Average_Daytime_LST_' + year,
    folder: 'MODIS_LST_Fahrenheit_Day',
    fileNamePrefix: 'Bexar_County_LST_' + year,
    region: bexarCounty.geometry(),
    scale: 1000,
    fileFormat: 'GeoTIFF',
    maxPixels: 1e9
  });

  // Reproject the LST image to UTM Zone 14N (adjust based on your study area)
  var reprojectedLST = LSTFahrenheit.reproject({
    crs: 'EPSG:32614',
    scale: 1000
  });

  var visParams = {min: 84, max: 107, palette: ['green', 'yellow', 'red']};
  Map.addLayer(reprojectedLST, visParams, 'Average Daytime LST ' + year);
}

// Loop over each year from 2013 to 2023
for (var year = 2013; year <= 2023; year++) {
  processYear(year);
}

Once aquired, send it to Google Drive and then download it onto the local machine as a TIFF,
Use QGIS to assess. In this case, that means get the low and high temperature points.
We need the temperatures so we can create a color scheme for the data viz. (Claude.ai helped create the colors I used).

2. Clean and process the data

I created a Jupyter notebook for this part. Here’s the code where I give it a color scheme that makes sense:

import rasterio
from rasterio.plot import show
from matplotlib import pyplot as plt
from matplotlib.colors import LinearSegmentedColormap

# Define the minimum and maximum values from your CSV
min_value = 69.61
max_value = 111.82

# Create a custom color map
colors = [
    '#006837', '#1a9850', '#66bd63', '#a6d96a', '#d9ef8b',
    '#ffffbf', '#fee08b', '#fdae61', '#f46d43', '#d73027'
]  # Define your desired colors
custom_cmap = LinearSegmentedColormap.from_list('custom', colors)

# Your raster file paths
raster_files = ['./MODIS_A1_daytimeLST/TIFF/Bexar_County_LST_2013.tif',
                './MODIS_A1_daytimeLST/TIFF/Bexar_County_LST_2014.tif',
                './MODIS_A1_daytimeLST/TIFF/Bexar_County_LST_2015.tif',
                './MODIS_A1_daytimeLST/TIFF/Bexar_County_LST_2016.tif',
                './MODIS_A1_daytimeLST/TIFF/Bexar_County_LST_2017.tif',
                './MODIS_A1_daytimeLST/TIFF/Bexar_County_LST_2018.tif',
                './MODIS_A1_daytimeLST/TIFF/Bexar_County_LST_2019.tif',
                './MODIS_A1_daytimeLST/TIFF/Bexar_County_LST_2020.tif',
                './MODIS_A1_daytimeLST/TIFF/Bexar_County_LST_2021.tif',
                './MODIS_A1_daytimeLST/TIFF/Bexar_County_LST_2022.tif',
                './MODIS_A1_daytimeLST/TIFF/Bexar_County_LST_2023.tif']

for raster_file in raster_files:
    # Open the raster file
    with rasterio.open(raster_file) as src:
        raster_data = src.read(1)

        # Normalize the raster data based on the min and max values
        normalized_data = (raster_data - min_value) / (max_value - min_value)

        # Apply the custom color map to the normalized data
        colored_data = custom_cmap(normalized_data, bytes=True)

        # Remove the alpha channel from the colored data
        colored_data = colored_data[:, :, :3]

        # Transpose the colored data to (bands, height, width)
        colored_data = colored_data.transpose((2, 0, 1))

        # Save the colored raster as a new file
        colored_raster_file = 'colored_' + raster_file.split('/')[-1]
        with rasterio.open(
            colored_raster_file,
            'w',
            driver='GTiff',
            width=src.width,
            height=src.height,
            count=3,  # RGB bands
            dtype=rasterio.uint8,
            crs=src.crs,
            transform=src.transform,
        ) as dst:
            dst.write(colored_data)

And here’s the code where I convert the GEOTiff files to MBTiles:

import os
import subprocess
import rasterio
from rasterio.crs import CRS
from rasterio.warp import calculate_default_transform, reproject, Resampling

# Your colored raster file paths
colored_raster_files = ['./MODIS_A1_daytimeLST/TIFF/colored/colored_Bexar_County_LST_2013.tif',
                        './MODIS_A1_daytimeLST/TIFF/colored/colored_Bexar_County_LST_2014.tif',
                        './MODIS_A1_daytimeLST/TIFF/colored/colored_Bexar_County_LST_2015.tif',
                        './MODIS_A1_daytimeLST/TIFF/colored/colored_Bexar_County_LST_2016.tif',
                        './MODIS_A1_daytimeLST/TIFF/colored/colored_Bexar_County_LST_2017.tif',
                        './MODIS_A1_daytimeLST/TIFF/colored/colored_Bexar_County_LST_2018.tif',
                        './MODIS_A1_daytimeLST/TIFF/colored/colored_Bexar_County_LST_2019.tif',
                        './MODIS_A1_daytimeLST/TIFF/colored/colored_Bexar_County_LST_2020.tif',
                        './MODIS_A1_daytimeLST/TIFF/colored/colored_Bexar_County_LST_2021.tif',
                        './MODIS_A1_daytimeLST/TIFF/colored/colored_Bexar_County_LST_2022.tif',
                        './MODIS_A1_daytimeLST/TIFF/colored/colored_Bexar_County_LST_2023.tif']

for colored_raster_file in colored_raster_files:
    # Open the colored raster file
    with rasterio.open(colored_raster_file) as src:
        # Define the target CRS for MBTiles (Web Mercator)
        dst_crs = CRS.from_epsg(3857)

        # Calculate the target transform and dimensions
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds)

        # Create the reprojected raster file path
        reprojected_raster_file = colored_raster_file.replace('.tif', '_reprojected.tif')

        # Reproject the colored raster to the target CRS
        with rasterio.open(
            reprojected_raster_file,
            'w',
            driver='GTiff',
            width=width,
            height=height,
            count=3,
            dtype=rasterio.uint8,
            crs=dst_crs,
            transform=transform
        ) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    resampling=Resampling.nearest,
                    num_threads=4
                )

        # Create the MBTiles file path
        mbtiles_file = colored_raster_file.replace('.tif', '.mbtiles')

        # Use the gdal_translate command to create the MBTiles file
        subprocess.run([
            'gdal_translate',
            '-of', 'MBTiles',
            '-co', 'TILE_FORMAT=JPEG',
            '-co', 'QUALITY=90',
            '-co', 'ZOOM_LEVELS=0-5',
            reprojected_raster_file,
            mbtiles_file
        ])

        # Remove the reprojected raster file
        os.remove(reprojected_raster_file)

Get the MBTiles onto Wheeljack and then create the URL:

Use this code in a Jupyter notebook:

baseURL = "https://maps.hdndevhub.com/"
extension = ".png"

def getTileURL(baseURL, folderPath, fileName, extension):
    return baseURL + folderPath + "tileserver.php?/" + fileName + ".json?/" + fileName + "/{z}/{x}/{y}" + extension

# Generate URLs for each year from 2013 to 2023
for year in range(2013, 2024):
    folderPath = f"devhub/tx-heat/day_{year}/"
    fileName = f"colored_Bexar_County_LST_{year}"
    
    tileURL = getTileURL(baseURL, folderPath, fileName, extension)
    print(tileURL)