Author: TRW (Page 3 of 4)

If you know my #WarTraces post, you’ll see that I use satellite imagery from the European Space Agency. Today I’ll show how to download and manipulate these images from mission Sentinel 2, as well as Landsat, being very similar to do.

I’ll do it in QGIS, so you need to download this software with the GRASS GIS module included.

Sign up in platforms.

You need to create an account to access both data sources, just filling a form and confirming via email.

Copernicus Open Access Hub

https://scihub.copernicus.eu/dhus/#/self-registration

NASA Earthdata Search

https://urs.earthdata.nasa.gov/users/new

Download data.

Let's make a difference according to the access platform.

Copernicus Open Access Hub

Access the viewer and log in with the new credentials.

https://scihub.copernicus.eu/dhus/

Now, draw a delimiting polygon of the area of interest. Every corner is set with right click, and the polygon is closed with double right click.

Next, click the filter icon in the search bar and tick the box for mission Sentinel-2 (you can also set dates, but it will be sorted from recent to old by default).

Finally, all available products will appear, so you can preview or download in a compressed .zip file.

NASA Earthdata Search

Like before, log in on the viewer.

https://search.earthdata.nasa.gov/search

Define an area of interest using the geometry tools in the lateral bar.

Filter the products, for example selecting "Imagery". Then select the HLS Landsat dataset in the catalogue.

In this case, we need to download each band separately, which we'll see next.

TCI image (true colour)

Simplifying (a lot), satellites capture images in different wave longitude ranges, called bands.

Also, digital images or pictures that we use to represent reality seen by the naked eye use to be defined in three layers: one for color Red, another for color Green, and another one for Blue, defining in each of these layers one value for every pixel with the amount of this color. The combination of all pixels is what simulates an image similar to reality.

TIP! Do the test and zoom in a lot in any picture and see the mosaic of pixels that make it up.

Lets use the bands captured by satellites in the visible wave range to generate an image as if we saw it from space.

Sentinel 2

In the case of Sentinel 2, we can show a true color image directly downloading the product L1C and using the layer _TCI.jp2 located in /GRANULE/DATASET_CODE_NAME/IMG_DATA/

The TCI image will show like this:

Landsat

TIP! https://gisgeography.com/landsat-8-bands-combinations/

For Landsat, we need to combine the different bands captured individually by the satellite in the visible wave range. These are:

RED	B04
GREEN	B03
BLUE	B02

Now in QGIS, we'll use function i.group from GRASS module which combines the selected layers in a multiband Red, Green, Blue layer.

To make it easier, I rename the layers like this, to make sure the script reads them properly:

Looking like this:

In this case it's nothing similar to the previous image or reality, so we'll adjust the properties to extract more color from the data:

• First, let's make sure we use all the data reading real max/min values rather than estimated.

• The image already looks better, but missing some brightness and color.

• So let's adjust Color Rendering parameters:

And it now looks more realistic:

SWIR (short wave infra-red)

Using bands in infra-red wave length, we can see further more than with the naked eye, or even across certain objects like some clouds or smoke.

During summer wildfires, I saw how Copernicus used these images to trace the extent of fires, and I thought they could be used to locate bomb impacts under the smoke they produce. That's how this ended in #WarTraces.

Sentinel 2

Using the process above, now we combine the following bands:

RED	B12
GREEN	B8A
BLUE	B04

Landsat

Make a new combination of these bands:

RED	B07
GREEN	B06
BLUE	B04

Conclusion

As a final note, it's important to highlight the pixel size, as it will affect image detail (resolution):

• Sentinel 2: 10 meters.
• Landsat: 30 meters.

On the other hand, Landsat provides images since 2013, while Sentinel only provides data from 2015.

Also, you can preview all these layers and combinations in online viewers like EO Browser by Sentinel Hub.

If you get stuck or want to comment, let me know on 🐦 Twitter!

🐦 @RoamingWorkshop

If you’ve used desktop GIS software, you will be surprised that many geospatial tools are also available to use in your web maps using a free and open-source javascript library.

Today I show you TurfJS. Press the button 🌤 in the map below and see one of many possibilities: inverse distance weighting or IDW.

Installation

TurfJS hosts a practical website explaining all available functions, as well as the initial setup:

https://turfjs.org/getting-started

To use Turf you only need to link the library in your map using the CDN link:

<script src='https://unpkg.com/@turf/turf@6/turf.min.js'></script>

Base map

You can create a basemap, for example, with Leaflet. TurfJS uses JSON sintaxis, so it's highly compatible with the geoJSON capabilities in Leaflet.

TIP! Leaflet has been recently updated to version 1.9.1. This post hasn't been updated and still uses version 1.7.1.

I will recycle this example from one of my first posts with a few changes:

• Add the Turf library link in the <head>.
• Define the markers in a data matrix or array, similar to how they are obtained from any data source in JSON format:

https://www.w3schools.com/js/js_arrays.asp

var markers=[
{"name":"Estacion1",
"lat":25.77,
"lon":-80.13,
"dato":26.4},
{"name":"Estacion2",
"lat":25.795,
"lon":-80.21,
"dato":29.8},
{"name":"Estacion3",
"lat":25.70,
"lon":-80.275,
"dato":31.35}
];

TIP! Numeric values won't go in quotation marks, otherwise Turf will take it as a String (something that doesn't happen in Leaflet).

• Add the markers to the map in a loop, using forEach(), widely used in these cases:

https://www.w3schools.com/jsref/jsref_foreach.asp

markers.forEach(function(m){
    var marker=L.marker([m.lat, m.lon]).addTo(map);
    marker.bindTooltip("<b>"+m.name+"</b><br>"+m.dato).openTooltip();
});

//WE ADDED A SINGLE MARKER LIKE THIS ON LEAFLET
//var marker = L.marker([25.77, -80.13]).addTo(map);
//marker.bindTooltip("<b>Estacion1</b><br>m.dato").openTooltip();

Looking like this:

and this one being its code:

<head>
<link rel="stylesheet" href="https://unpkg.com/[email protected]/dist/leaflet.css"
  integrity="sha512-xodZBNTC5n17Xt2atTPuE1HxjVMSvLVW9ocqUKLsCC5CXdbqCmblAshOMAS6/keqq/sMZMZ19scR4PsZChSR7A=="
  crossorigin=""/>
<script src="https://unpkg.com/[email protected]/dist/leaflet.js"
  integrity="sha512-XQoYMqMTK8LvdxXYG3nZ448hOEQiglfqkJs1NOQV44cWnUrBc8PkAOcXy20w0vlaXaVUearIOBhiXZ5V3ynxwA=="
  crossorigin=""></script>
<script src='https://unpkg.com/@turf/turf@6/turf.min.js'></script>
</head>

<style>
#map{border-radius: 30px;}
</style>

<body>
<div id="map" style="height:500px;"></div>
</body>

<script>
var map = L.map('map').setView([25.75, -80.2], 12);
L.tileLayer("http://a.tile.openstreetmap.org/{z}/{x}/{y}.png", {/*no properties*/
}).addTo(map);

var markers=[
{"name":"Estacion1",
"lat":25.77,
"lon":-80.13,
"dato":26.4},
{"name":"Estacion2",
"lat":25.795,
"lon":-80.21,
"dato":29.8},
{"name":"Estacion3",
"lat":25.70,
"lon":-80.275,
"dato":31.35}
];

markers.forEach(function(m){
    var marker=L.marker([m.lat, m.lon]).addTo(map);
    marker.bindTooltip("<b>"+m.name+"</b><br>"+m.dato).openTooltip();
});
</script>

IDW interpolation with TurfJS

In the form above, data is now easily used by TurfJS.

One of the most interesting tools is the interpolation of data in a regular grid.

https://turfjs.org/docs/#interpolate

This is done using a very usual method: inverse distance weighting (IDW).

Following the example given by TurfJS:

1. Convert data into a Turf point layer, a featureCollection:

https://turfjs.org/docs/#featureCollection

//empty points array
var points=[];

//loop to fill the array with data
markers.forEach(function(m){
    points.push(turf.point([m.lon, m.lat],{name: m.name, dato: m.dato}));
});

//convert to Turf featureCollection
points=turf.featureCollection(points);

TIP! Note that Turf taked coordinates as longitude-latitude, instead of the usual latitude-longitude!

2. Define interpolation options, having the following:
   • gridType: or cell shape
     • points
     • square
     • hex
     • triangle
   • property: field that contains the interpolation values.
   • units: spatial units used in the interpolation:
     • miles
     • kilometers
     • radians
     • degrees
   • weight: power used in the inverse distance weighting, usually 2, although higher values will make a smoother result.

//interpolation options: rectangular grid using 'dato' variable in kilometers using the power of 2

var options = {gridType: 'square', property: 'dato', units: 'kilometers', weight: 2};

3. Finally, generate the interpolated grid and add it to the map, being the numeric value the cell size of the grid (0.5 Km):

//create Turf grid
var malla=turf.interpolate(points, 0.5, options);

//add to Leaflet
L.geoJSON(malla).addTo(map);

Generating this result (nothing conclusive):

Format grid and show values

Seems not, but the interpolation is done, it's just shown incorrectly.

I'll add a color ramp to represent data, and also show values when hovering the cells. Let's modify the Leaflet geoJSON object properties for this:

• style: function to modify object style, according to its value.
• onEachFeature: function added to each object to interact with it.

https://leafletjs.com/reference.html#geojson

https://leafletjs.com/examples/choropleth/

L.geoJSON(malla,{
  style:function(feature){
  var val=feature.properties.dato;
  if(val>30){
    return {fillColor: "orangered", fillOpacity: 0.5, weight:0};
  }else if(val>28 && val < 30){
    return {fillColor: "yellowgreen", fillOpacity: 0.5, weight:0};
  }else if(val>26 && val < 28){
    return {fillColor: "darkgreen", fillOpacity: 0.5, weight:0};
  }
},
onEachFeature:function(feature,layer){
  layer.on({
  mouseover:function(e){
    val=e.target.feature.properties.dato.toFixed(2);
    e.target.bindTooltip(val).openTooltip();
  }
  });
}
}).addTo(map);

Ending like this:

And here is the full code:

<head>
<link rel="stylesheet" href="https://unpkg.com/[email protected]/dist/leaflet.css"
  integrity="sha512-xodZBNTC5n17Xt2atTPuE1HxjVMSvLVW9ocqUKLsCC5CXdbqCmblAshOMAS6/keqq/sMZMZ19scR4PsZChSR7A=="
  crossorigin=""/>
<script src="https://unpkg.com/[email protected]/dist/leaflet.js"
  integrity="sha512-XQoYMqMTK8LvdxXYG3nZ448hOEQiglfqkJs1NOQV44cWnUrBc8PkAOcXy20w0vlaXaVUearIOBhiXZ5V3ynxwA=="
  crossorigin=""></script>
<script src='https://unpkg.com/@turf/turf@6/turf.min.js'></script>
</head>

<style>
#map{border-radius: 30px;}
</style>

<body>
<div id="map" style="height:500px;"></div>
</body>

<script>
var map = L.map('map').setView([25.75, -80.2], 12);
L.tileLayer("http://a.tile.openstreetmap.org/{z}/{x}/{y}.png", {/*no properties*/
}).addTo(map);

var markers=[
{"name":"Estacion1",
"lat":25.77,
"lon":-80.13,
"dato":26.4},
{"name":"Estacion2",
"lat":25.795,
"lon":-80.21,
"dato":29.8},
{"name":"Estacion3",
"lat":25.70,
"lon":-80.275,
"dato":31.35}
];

markers.forEach(function(m){
    var marker=L.marker([m.lat, m.lon]).addTo(map);
    marker.bindTooltip("<b>"+m.name+"</b><br>"+m.dato).openTooltip();
});

//empty points array
var points=[];

//loop to fill the array with data
markers.forEach(function(m){
    points.push(turf.point([m.lon, m.lat],{name: m.name, dato: m.dato}));
});

//convert to Turf featureCollection
points=turf.featureCollection(points);

//interpolation options: rectangular grid using 'dato' variable in kilometers using the power of 2

var options = {gridType: 'square', property: 'dato', units: 'kilometers', weight: 2};

//create Turf grid
var malla=turf.interpolate(points, 0.5, options);

malla=L.geoJSON(malla,{
  style:function(feature){
  var val=feature.properties.dato;
  if(val>30){
    return {fillColor: "orangered", fillOpacity: 0.5, weight:0};
  }else if(val>28 && val < 30){
    return {fillColor: "yellowgreen", fillOpacity: 0.5, weight:0};
  }else if(val>26 && val < 28){
    return {fillColor: "darkgreen", fillOpacity: 0.5, weight:0};
  }
},
onEachFeature:function(feature,layer){
  layer.on({
  mouseover:function(e){
    val=e.target.feature.properties.dato.toFixed(2);
    e.target.bindTooltip(val).openTooltip();
  }
  });
}
}).addTo(map);

</script>

Applications

Applying this to real data, I modified the leafMET library to add the interpolation using isobands.

It's the map you saw at the beginning which you can see fullscreen:

https://theroamingworkshop.cloud/leafMET/leafMET+Turf.html

The code is also available as a new branch of leafMET in github:

https://github.com/TheRoam/leafMET/tree/leafMET+Turf

And that's all! Hope you find it useful!

Any doubts or comments are awaited on 🐦 Twitter!

🐦 @RoamingWorkshop

In this post I explained how to directly add AEMET stations with certain data, using my leafMET plugin for Leaflet (available on github):

https://github.com/TheRoam/leafMET

Let’s now pay attention to the detail of the API and how to access all the available data.

API access

To use the API we need an API key. Go to the AEMET OpenData webpage and click on "Solicitar" to request it.

Next we'll be asked for an email address where the key will be sent.

You'll receive a first confirmation email, and then a second email with the access key. Let's copy it and keep moving.

Documentation and examples

Get back to AEMET OpenData and into the developer's section "Acceso Desarrolladores".

We're interested in the dynamic documentation, which shows all the data requests available that we can try live.

Click on each topic to display the syntax for every data request.

Click on observación-convencional and then /api/observacion/convencional/todas.

Now click on the red sign on the right, where you'll paste and authorize the API key.

Now press the Try it out! button and you'll get:

• a curl request example
• the request URL
• the body, code and headers of the response

If we open the URL inside the "datos" field we'll see the full list of data for all the stations. A JSON lot like this:

Loading...

The main fields are the following:

• idema: station ID or code.
• lat and lon: latitude and longitude coordinates.
• alt: station altitude.
• fint: data interval date.
• prec: precipitation in this period.
• vv: wind speed in the station.
• ta: ambience temperature in the station.

Data request using javascript

HTML structure

I'll develop an example to access the data and make use of them.

I'll simply create a .html document with a request button, a text line and a request function:

Request

Click the button to request data

<html>
<body>
<button id="solicitud" onclick="solicitar()">Request</button>
<p id="texto">Click the button to request data</p>
</body>
<script>
function solicitar(){
document.getElementById("texto").innerHTML="I didn't program that yet..";
}
</script>
</html>

TIP! Copy and paste the code above in a text document with .html extension and open it in your browser.

HTTP request

We'll do the data request in javascript using the XMLHttpRequest object. The example from w3schools is really simple:

https://www.w3schools.com/js/js_ajax_http.asp

As we saw in the dynamic documentation, after doing the request, we got a link to the data, which is another request, no we need to nest two requests like this:

<script>
function solicitar(){
// Define loading message
document.getElementById("texto").innerHTML="Cargando datos...";
// Define API Key
var AK=tUApiKey;
// Define request URLs
var URL1="https://opendata.aemet.es/opendata/api/observacion/convencional/todas";
var URL2="";

// Create XMLHttpRequest objects
const xhttp1 = new XMLHttpRequest();
const xhttp2 = new XMLHttpRequest();

// Define response function for the first request:
// We want to access te "datos" URL
xhttp1.onload = function() {
  // 1º Parse the response to JSON:
  URL2=JSON.parse(this.responseText);
  // 2º Get the "datos" field:
  URL2=URL2.datos;  
  // 3º Make the new request:
  xhttp2.open("GET", URL2);
  xhttp2.send();
}

// Define response function for the second request:
// We'll modify the text line with some information in the data
xhttp2.onload = function() {
  // 1º Parse the response to JSON (much better to work with):
  var datos=JSON.parse(this.responseText);

  // 2º Get the length of the JSON data
  // This is equivalent to the individual hourly entries per stations
  var registros=datos.length;

  // 3º Get first entry date.
  // We'll format it as it's in ISO standard
  var fechaini=datos[0].fint;
  fechaini=new Date(fechaini+"Z");
  fechaini=fechaini.toLocaleDateString()+" "+fechaini.toLocaleTimeString();

  // 4º Get last entry date.
  // We'll format in the same way
  var fechafin=datos[(datos.length-1)].fint;
  fechafin=new Date(fechaini+"Z");
  fechafin=fechafin.toLocaleDateString()+" "+fechafin.toLocaleTimeString();

  // 5º Merge it all and display in the text line
  document.getElementById("texto").innerHTML="Se han descargado <b>"+registros+"</b> registros desde el <b>"+fechaini+"</b> hasta el <b>"+fechafin+"</b>. <a href='"+URL2+"' target='_blank'>Ver todos los datos.</a>";
}

// Send the first request
xhttp1.open("GET", URL1+"/?api_key="+AK);
xhttp1.send();

}
</script>

Let's see the result below:

Request

Click the button to request data

And done! This way you can obtain coordinates and values for different metheorological stations from AEMET and drop them in a table or web map like we've seen for Leaflet here.

Leave any comments or queries on Twitter!🐦

🐦 @RoamingWorkshop

In July, Nature published a high definition Digital Surface Model of the new orography built by La Palma’s Tajogaite volcano eruption in 2021. A spectacular job by Istituto Nazionale di Geofisica e Vulcanologia together with Instituto Vulcanológico de Canarias.

The usual way to see and work with these files is using a GIS software like QGIS or ArcGIS, but how can we make this friendly to standard users? Wouldn’t it be better seen directly in 3D and visit the place virtually?

At this point is where I found CesiumJS: a 3D map viewer styling [G] Earth but open-source, free and highly customizable to make your own projects.

Create a Cesium ion account

To use CesiumJS you need to register an account in order to obtain an Access Token.

Registering, we can also upload our own files to customize our 3D map creations. In my case you can see that I uploaded the Digital Surface Model and imagery of the viewer. I drop the links in case you want to use the same files:

• Digital Surface Model of La Palma volcano in January 2022.
  Format geoTIFF with 0'2 metre grid size.
  https://portal.opentopography.org/dataspace/dataset?opentopoID=OTDS.032022.4083.1#panel_head_files
• 50cm satellite imagery from January 2022 by IDE Canarias
  https://idecan1.grafcan.es/ServicioWMS/ErupcionLP2021?SERVICE=WMS&request=GetMap&version=1.1.1&layers=ORTO_14012022_SAT50&format=image%2Fjpeg&transparent=true&width=1200&height=800&srs=EPSG:4326&bbox=-17.93,28.58,-17.83,28.65

TIP! You can use the WMS service in your GIS software to extract the imagery, or download the one I use in this example (270 MB with 1m resolution)

Include CesiumJS in your web.

To use CesiumJS we configure a .html file as usual.

In this case, I'll copy the code directly from their quickstart example as it is only about 25 lines:

<!DOCTYPE html>
<html lang="en">
<head>
  <meta charset="utf-8">
  <!-- Include the CesiumJS JavaScript and CSS files -->
  <script src="https://cesium.com/downloads/cesiumjs/releases/1.97/Build/Cesium/Cesium.js"></script>
  <link href="https://cesium.com/downloads/cesiumjs/releases/1.97/Build/Cesium/Widgets/widgets.css" rel="stylesheet">
</head>
<body>
  <div id="cesiumContainer"></div>
  <script>
    // Your access token can be found at: https://cesium.com/ion/tokens.
    // This is the default access token from your ion account
    Cesium.Ion.defaultAccessToken = 'eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJqdGkiOiI1MTYwNGM0Yi1iYzlkLTRkMTUtOGQyOS04YjQxNmUwZDQ0YjkiLCJpZCI6MTA0Mjc3LCJpYXQiOjE2NjAxMDk4MzR9.syNWeKPLWA2eMrEh4K9hvkyp1oGdbLMaL0Ozk1kaksk';
    // Initialize the Cesium Viewer in the HTML element with the `cesiumContainer` ID.
    const viewer = new Cesium.Viewer('cesiumContainer', {
      terrainProvider: Cesium.createWorldTerrain()
    });    
    // Add Cesium OSM Buildings, a global 3D buildings layer.
    const buildingTileset = viewer.scene.primitives.add(Cesium.createOsmBuildings());   
    // Fly the camera to San Francisco at the given longitude, latitude, and height.
    viewer.camera.flyTo({
      destination : Cesium.Cartesian3.fromDegrees(-122.4175, 37.655, 400),
      orientation : {
        heading : Cesium.Math.toRadians(0.0),
        pitch : Cesium.Math.toRadians(-15.0),
      }
    });
  </script>
 </div>
</body>
</html>

<body style="margin:0;width:100vw;height:100vh;">
<div style="height:100%;" id="cesiumContainer"></div>

viewer.timeline.container.style.visibility = "hidden";
viewer.animation.container.style.visibility = "hidden";

viewer.camera.flyTo({
    destination : Cesium.Cartesian3.fromDegrees(-17.8195, 28.6052, 3000),
    orientation : {
    heading : Cesium.Math.toRadians(-85.0),
    pitch : Cesium.Math.toRadians(-30.0),
    }
});

<html lang="en">
<head>
<meta charset="utf-8">

<script src="https://cesium.com/downloads/cesiumjs/releases/1.96/Build/Cesium/Cesium.js"></script>
<link href="https://cesium.com/downloads/cesiumjs/releases/1.96/Build/Cesium/Widgets/widgets.css" rel="stylesheet">
</head>
<body style="margin:0;width:100vw;height:100vh;">
<div style="height:100%;" id="cesiumContainer"></div>
<script>
    // Your access token can be found at: https://cesium.com/ion/tokens.
    // Replace `your_access_token` with your Cesium ion access token.

    Cesium.Ion.defaultAccessToken = 'eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJqdGkiOiI1MTYwNGM0Yi1iYzlkLTRkMTUtOGQyOS04YjQxNmUwZDQ0YjkiLCJpZCI6MTA0Mjc3LCJpYXQiOjE2NjAxMDk4MzR9.syNWeKPLWA2eMrEh4K9hvkyp1oGdbLMaL0Ozk1kaksk';

    // Initialize the Cesium Viewer in the HTML element with the `cesiumContainer` ID.
    const viewer = new Cesium.Viewer('cesiumContainer', {
      //terrainProvider: Cesium.createWorldTerrain()//terreno original
		terrainProvider: new Cesium.CesiumTerrainProvider({//terreno modificado
			url: Cesium.IonResource.fromAssetId(1255858),
		}),	  
    });    
	
	// Esconder widgest inferiores
	viewer.timeline.container.style.visibility = "hidden";
	viewer.animation.container.style.visibility = "hidden";
	
	// Add Sentinel2 imagery
	const layer = viewer.imageryLayers.addImageryProvider(
	  new Cesium.IonImageryProvider({ assetId: 1256129 })
	);
    // Add Cesium OSM Buildings, a global 3D buildings layer.
    const buildingTileset = viewer.scene.primitives.add(Cesium.createOsmBuildings());   
    // Fly the camera to San Francisco at the given longitude, latitude, and height.
    viewer.camera.flyTo({
      destination : Cesium.Cartesian3.fromDegrees(-17.8195, 28.6052, 3000),
      orientation : {
        heading : Cesium.Math.toRadians(-85.0),
        pitch : Cesium.Math.toRadians(-30.0),
      }
    });
  </script>
</div>
</html>