In which we use raster data with Leaflet, and undertake spatial analysis to understand green space availability in Kampala, Uganda.

Part 1

Raster Data

Before we do anything, make a new web page using the usual template, name it something like week12.html

This week we will be working in the city of Kampala, Uganda - set the map() constructor to the following:

map = L.map('map').setView([0.348384, 32.568841], 14);

We want some sattelite imagery rather than a map this week, so Change the L.tileLayer() constructor to the following:

L.tileLayer('https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}', {
					attribution: 'Tiles © Esri — Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community'
}).addTo(map);

As you have seen in the course so far, web GIS is dominated by vector data, and it is very unusual to see raster data being used for anything more than adding as an image overlay to a map. There are a number of reasons for this, chief amongst are the difficulty in getting raster data into a format that works well in the browser; and the complexity of dealing with raster data compared to vector. However, there is a great deal of potential for the use of raster data in web GIS, and so to facilitate this I have created a piece of software that allows you to convert any raster dataset into a JSON object for use in JavaScript. Even better than that, the resulting object includes lots of functions and properties that make working with raster data quite simple - effectively removing both of the major barriers to using raster data in web GIS! In doing so we will gain some valuable experience both in working with rasters and in using new libraries.

Have a read over the README file in my library to understand how the object works.

Before we do anything else we therefore need to get our data. To save us some time I have already converted our data into a javaScript object for you, so all that you need to do is add it and proj4.js to your page:

Click here to access the dataset and press Ctrl+S to save the resulting page as kampala.js in the same directory as your week12.html file

Add script import tags for https://cdnjs.cloudflare.com/ajax/libs/proj4js/2.5.0/proj4.js and kampala.js - make sure that it is in that order, as kampala.js will not work unless proj4.js has already been imported

That is all that you need to do in order to load your raster data - not bad eh! To make sure that it has loaded okay, let’s test the properties - here is a reminder of what they are:

To make sure that it has loaded properly and that all is as it should be, let’s print out some of the properties of the layer:

Add the below snippet to an appropriate location in your script to see the propeties of the dataset

console.log(kampala.bands);
console.log(kampala.bl);
console.log(kampala.tl);
console.log(kampala.br);
console.log(kampala.tr);
console.log(kampala.bounds);
console.log(kampala.resolution);
console.log(kampala.width);
console.log(kampala.height);
console.log(kampala.proj);
console.log(kampala.transformer);

Did it work?

If so, you should be able to see all of the information about the dataset, including that it has a resolution of 10m and that we have two bands available:

• band 3 (red)

• band 4 (near infra red)

  We have left out all of the other (unused) bands in order to save space.

The Normalised Difference Vegetation Index (NDVI)

The reason that we are interested in the red and near infra red bands in particular is that we are going to do some remote sensing analysis to understand verying levels of availability of green space across the city of Kampala. In order to be able to acheive this, we need to be able to define ‘green space’, and we will do that using the NDVI, which is a simple equation:

This returns a value of between -1 (definitely not green, probably water) and 1 (definitely green). This will be the basis of our entire analysis!

Convert the above equation into a function called ndvi() that takes in two arguments: red and nir

Now we have done this, we can calculate NDVI. We have a few choices for how we want to do this - we could calculate it for a little bit of the dataset each time that we want to know some values (real-time calculation: e.g. when a user clicks on the map), or we could calculate it all in advance and store the values in a global variable to be retrieved later (pre-calculation i.e. when the page loads). Making decisions like this is always a trade off - but as a rule of thumb you would normally use real-time calculation for very small calculations (which won’t impact on the user experience) or very large calculations (which are too big to pre-calculate without making the user wait a long time). In this case, we are somewhere in the middle, so we will pre-calculate the NDVI values and store them for retrieval later on:

Add a global variable called ndviLayer

Complete the below function and add it to your code - it should loop through the whole dataset and calculate the NDVI for each individual cell, storing the result in a new array that is the same dimensions as the dataset:

(The trick is that you need to think about how you can access the number of columns and the number of rows in the dataset - have a look at the properties above to find the answer…)

/**
 * Calculate NDVI for entire dataset
 */
function calculateNdviLayer(){

  // initialise global array for ndvi results
  ndviLayer = [];

  // loop through every single column
  for(let x = 0; ){

    // add a new column
    ndviLayer.push([]);

    // loop through every single pixel in the column
    for(let y = 0; ){

      // calculate ndvi for each pixel
      ndviLayer[x][y] = ndvi(kampala.band3[x][y], kampala.band4[x][y]);
    }
  }
}

Once you have done this, you should have finished part 1 - all that you need to do is verify that it has worked.

Paste the below snippet into your code at an appropriate location in order to check that it has worked:

// check the width of the dataset and the ndvi layer - both numbers should be the same
console.log(kampala.width, ndviLayer.length);

// check the height of the dataset and the ndvi layer - both numbers should be the same
console.log(kampala.height, ndviLayer[0].length);

// check a random point in the dataset
console.log(ndvilayer[100][100]);	

If all has gone to plan, you should see this:

1766 1766
1711 1711
0.4264202600958248

If so, then well done!! Now it’s time to go back to the lecture…

Part 2

Our main goal today is to produce a website that can quantitatively assess the availability of green space within the city of Kampala, Uganda. This is important, and the Kampala City Authorities are currently trying very hard to identify areas of low green space availability so that they can address it in order to improve the health of their citizens. A simple (if rather arbitrary) measure for green space availability which is common in the literature is the amount of green space within 500m of a given location, expressed as a green ratio. We can use NDVI to calculate this by setting a simple threshold for what we consider to be ‘Green’ (say, 0.5 NDVI) and representing this a ratio of the number of cells of ‘Green’ against the number of cells within 500m of a given location:

To do this, we simply need to set a click listener that:

1. identifies all of the cells within 500m of the click location
2. assesses them all for ‘greenness’
3. calculates the Green Space Ratio and reports it to the screen

That doesn’t sound so difficult, right? With a little help from my raster-js library, turf.js and a new library called chart.js, we should be able to do this quite easily, but first we need to make sure that we understand a little more about rasters…

Working with raster data

As we discussed in the lecture, there is a bit of maths to do in order to be able to relate WGS84 (longitude, latitude) coordinates to the corresponding location in a 2D array of raster data. It works like this:

Spherical Coordinate Space ↔ Projected Coordinate Space ↔ Image Space

So to go from coordinates in spherical coordinate space (latitude longitude) to an array position in Image Space (and vice-versa), you need to go via projected coordinate space. Just like in week 2, we can use proj4js in order to transform between spherical coordinate space data (used by leaflet) and projected coordinate space (in which the dataset is stored). There is no library to convert between projected coordinate space (in which the dataset is stored) and image space (which relates to the array position needed to access data in JavaScript), but this can be achieved manually using some relatively simple maths.

Fortunately, my library contains all of the code that you need in order to be able to do this, but it is still vitally important that you understand it, otherwise you won’t be able to use it effectively, which can be very frustrating!

Have a close look at the functions in the object and make sure that you can see the above processes working

Now that we have access to functions like this, we can start to interact with the dataset in a more meaningful way:

Add a click listener to the map that logs NDVI value for that location to the console - use the snippet below to help you:

const imageCoords = geo2image([e.latlng.lng, e.latlng.lat]);
console.log(ndviLayer[imageCoords[0]][imageCoords[1]]);

Does it go up and down as expected if you click on green and non-green areas?

Calculating the Green Ratio

In order to calcuate the green ratio for a 500m area surrounding our click location, we are going to need a bit of spatial analysis using turf.js.

Import turf.js (https://npmcdn.com/@turf/turf/turf.min.js)

Make a function called calculateGreenRatio() that takes in a single argument called lngLat

Add the following function call to your click listener:

// extract the coordinates of the click and pass to the function
const greenRatio = parseFloat((calculateGreenRatio([e.latlng.lng, e.latlng.lat])* 100).toFixed(2));

Note how we have used *100 to convert the result to a percentage, .toFixed(2) to convert the number to a String with only two decimal places, and then parseFloat() to convert it back to a number again. This is a cheeky workaround that gives us the effect of turning messy numbers like 0.83635268293836252723893 that we might get back from our function to much nicer looking numbers like 83 (meaning 83%). Nice touches like this are important for when you are reporting values back to users (espacially in Assessments…).

In order to do this, we will construct a 500m circular buffer around our click location, and then draw it to the map.

Add a global variable called circleLyr then add the below code to calculateGreenRatio()

// get a polygon representing a 500m (0.5km) radius around the clicked location
const circle = turf.circle(turf.point(lngLat), 0.5, {
  units: 'kilometers',
  steps: 64,	// how many points make up the circle (more = smoother circle)
});

// remove previous circle layer if there is one 
if (circleLyr) map.removeLayer(circleLyr);

Now add the circle to the map as a GeoJSON layer using the L.geoJson() constructor and store the result in the circleLyr global variable

Does it work? If so you should now get a circle every time you click.

Now for the tricky part - we need to get the values out of all of the cells in the raster dataset that fall within our circle - this requires some spatial analysis, which turf can handle for us - but unfortunately our circle is in spherical WGS84 coordinates, our dataset is in projected UTM (Universal Transverse Mercator) coordinates and is only accessible to JavaScript via in the form of a 2D array in image space - this is where our chain comes in handy:

Spherical Coordinate Space ↔ Projected Coordinate Space ↔ Image Space

Obviously, to be able to do our calculations we must have both datasets in the same coordinate system (either spherical, projected or image); it doesn’t really matter which way around we do this, as long as we pick one.

Point in Polygon analysis - pre-calculation with Bounding Box

Because we are going to work out which datasets fall within the circle, it is also highly beneficial if we can reduce the number of points that we test, because testing if a point is within a polygon is a computationally expensive test and if we test the entire dataset it will result in quite a slow answer.

We will therefore use turf to calculate the bounding box of the circle. We can then convert the corners of the bounding box to image space and use those coordinates to extract only the points from the array that are in the immediate vicinity of the circle. This means that we don’t need to test anything that is further away. The process is quite simple

Add the below statement to calculateGreenRatio() calculate the image space coordinates of the bottom left corner, and then add another similar statement to calculate the image space coordinates of the top right corner and store in a variable called tl

// get bounds of the circle and extract the corners in image space
const bl = kampala.geo2image([circleLyr.getBounds().getWest(), circleLyr.getBounds().getSouth()]);

Now we have these, we can loop through the sections of the array that we are interested in, convert each cell into spherical (WGS84) coordinates so that they are comparable with the circle and store as a turf.js point. While we’re at it, we can add a property to each point recording it’s NDVI value.

Getting the bounds in this way will give us a square that perfectly fits the circle (e.g. the red square and the blue circle below) - this is really useful, as we can easly extract a square from our dataset, wehere is would be very difficult to extract a circle.

Add the below to calculateGreenRatio()

// initialise an array to store the resulting points
let points = [];

// loop through each column (x) and row (y)
for (let x = bl[0]; x <= tr[0]; x++) {
  for (let y = tr[1]; y <= bl[1]; y++){

    // add point with property for ndvi value
    points.push(turf.point(kampala.image2geo([x, y]),{ ndvi: ndviLayer[x][y] }));
  }
}

We now have an array called points that contains a small square subset of locations from our raster dataset that immediately surrounds our circle. To illustrate, here are the resulting points drawn as markers onto the map (we aren’t going to do this as it takes ages to draw them all!).

Now we only have to assess those few points for whether or not they fall within the circle, rather than the entire dataset!! To do this, we convert our points array into a turf.featureCollection() and extract the ones that are within the circle using the turf.pointsWithinPolygon() function.

// get the coordinates from the array that are in the aoi
const pointsWithin = turf.pointsWithinPolygon(turf.featureCollection(points), circle);

Now we have only the points from the dataset that are actually inside the circle:

Now we have this, we just need to loop through each feature in the pointsWithin featureCollection, and count how many of them are > 0.5, and use this to calculate and return the green ratio.

Add the below snippet to your code and complete it (you will need to edit the conditional statement to work out if the current cell in the loop is >= 0.5).

// count the number of green cells (NDVI >= 0.5 threshold)
let green = 0;

// loop through each cell
for (let p = 0; p < pointsWithin.features.length; p++) {

  // is the current cell greater than or equal to 0.5?
  if ( ) {
    
    // if so, then increment the green counter
    green++;
  }
}

// calculate green ratio and return
return green / pointsWithin.features.length;

Edit the console.log() statement to your click listener so that it logs the green ratio rather than the ndvi.

Using Chart.js for Real-time Reporting

This is great! The only remaining question is how we can display our results in an effective way. In order to do this, we are going to write the result onto the page and display it visually using an animated chart.

Replace the <div id='map'></div> element in your HTML with the below snippet:

<!-- container to control the layout -->
<div id="container">

  <!-- the map -->
  <div id='map'></div>

  <!-- container for the pie chart -->
  <div id="canvas-holder">

    <!-- chart -->
    <canvas id="chart"></canvas>

    <!-- report text -->
    <div id='green-ratio'>
      <p>Click the map to calculate Green Ratio</p>
    </div>
  </div>
</div>

And update the CSS accordingly:

/* style the map container*/
#map {
  width: 400px;
  height: 500px;
  float: left;
}

/* style the chart container */
#canvas-holder {
  width: 400px;
  float: right;
}

/* style the container */
#container {
  width: 800px;
}

/* style the green ration report text */
#green-ratio {
  font-family: sans-serif;
  text-align: center;
}

Now, we are going to load our chart into chart and our text result into green-ratio.

Add your result in text to greenRatio so that it says n% Green Space (where n is the number)

Now to make our chart. We are going to use an excellent Open Source library called Chart.js to do this. You can use this library to make almost any kind of chart quite easily - it is well worth checking out ahead of Assessment 2… You can see many examples of the types of charts that yoiu can make here, and the instructions for use are here.

Firstly, import the Chart.js library from https://www.chartjs.org/dist/2.9.3/Chart.min.js

Chart.js works like this:

1. Get a reference (known as a context) to the HTML element (<canvas>) where you are going to draw the chart
2. Create anobject containing all of the options that define what the chart will look like
3. Use the Chart() constructor to build a chart in your desired context (ctx), using your defined options (chartConfig)

To build a chart in your chart canvas, and with no data in it, add the below to your initMap() function:

// get the reference to the chart div
const ctx = document.getElementById('chart').getContext('2d');

// initialise chart data
let chartConfig = {
  type: 'pie',
  data: {
    datasets:[{
      data: [],	// initialise with empty array as no data yet
      backgroundColor: ["#999999", "#27b827"],
    }],
    labels: ['Not Green', 'Green'],
  },
};

// initialise the chart
let pieChart = new Chart(ctx, chartConfig);

This is a pie chart with no data in it - which obviously can’t be drawn, so it will look something like this:

Note that there is no chart drawn, and that both of the classes in the legend are crossed out. This is exactly what we are expecting, safe in the knowledge that the chart will spring into life once we give it some data. To illustrate this, click the button below to add some data to the above chart.

Click Me!

Now all that we have to do is to update the data in chartConfig to represent the percentage of green and the non-green within the circle, and then tell the chart to update itself.

Add the below to your click listener

// update the chart to reflect the ratio
chartConfig.data.datasets[0].data = [100 - greenRatio, greenRatio];
pieChart.update();

Now, when you click on the map, it will update the chart and the reporting text! If all has gone to plan, it should look like the below:

If so, great work!!

That’s all folks!

Thankyou very much for a great module - I have really enjoyed teaching this class and it is great to see how much progress you have all made. Just think that most of you had never even programmed before only 3 months ago!

Happy coding!!! (and Merry Christmas)

Jonny.