Meteorology 101: How To Read Data Elevation Model On Basemap
Well, Hello Friends!
Welcome to my fifth post explaining Meteorology 101! Meteorology 101 is a series that explains basic tutorial ‘how to’ code about meteorology. Previously, I wrote about How To Plot Windrose using two different ways of visualization the heatmap. You can check it out if you missed. And, if you missed another topics about Meteorology 101, you can read it by visiting my Medium home profile.
First of all, this topic is a next story from the last topic — which still about my final project when i was in college. When I took my fourth final project-defense, my examiner asked me how the topography affected the wind speed at the site. So my answer is in the Figure 1!!
But, at that time I was confused about how to make an elevation map in python because I usually use ArcGIS/QGIS. So, I searched some tutorial on Google about how to make an elevation map in python. The result I found is that you can use GDAL or Geospatial Data Abstraction Library. But after I tried using GDAL, there were some problems which are location coordinates not appearing on the map, and also I had to overlay with the wind direction variable so I can make an analysis of how the topography affects the wind speed.
I will show you the problem and the solution that I found.
But before that, let me give a litle introduction about what the elevation data is. Elevation data is a digital representation of the Earth’s terrain. Data can be represented as points (x,y,z) or in a gridded format (raster). Digital Elevation Models (DEMs) are gridded raster where each pixel of the image has a value, the elevation above a model of the Earth (WGS84 Ellipsoid, for example). The size of each rectangular pixel (e.g. 2 meters x 2 meters) is referred to as the spatial resolution.
DOWNLOAD THE DATA
First, we need to download the elevation data. The data that we will use is from SRTM or Shuttle Radar Topography Mission. SRTM is elevation data on a near-global scale, acquired from the Space Shuttle, to generate the most complete high-resolution digital topographic database of Earth. You can read overview of SRTM data here.
There are several websites that provide the SRTM data. First, you can download it officially from earthexplorer.usgs.gov. There you will get all products from remote sensing SRTM. Second, if you want to get data in simple way, you can download from dwtkns.com/srtm30m/. Lastly, for Indonesian region, you can download the data from www.indonesia-geospasial.com/2020/01/download-dem-srtm-30-meter-se-indonesia.html, which has been divided into data per province.
PLOT THE ELEVATION DATA
After downloading the data, let’s open the jupyter notebook!
First, i will show you how do I plot with GDAL library. If you haven’t installed the library yet, I suggest you to read the tutorial here or watch this youtube tutorial, because installing this library is quite difficult. And, if you have finished installing the library, let’s into the code!
This is the library you have to declared.
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import xarray as xr
from mpl_toolkits.basemap import BasemapThen, for the GDAL library you have to declare this.
from osgeo import gdalNext, let’s open the data that we have downloaded.
gdal.UseExceptions()
filename = "pulausumba.tif"
gdal_data = gdal.Open(filename)
gdal_band = gdal_data.GetRasterBand(1)
nodataval = gdal_band.GetNoDataValue()By default the GDAL/OGR Python bindings do not raise exceptions when errors occur. Instead they return an error value such as None and write an error message to sys.stdout. You can enable exceptions by calling the UseExceptions() function (GDAL/OGR General). Next, define filename with name of the data. In this case, my data file name is “pulausumba” where the data that I used is the elevation data of Sumba Island. Then, define band by taking band 1 and define no data value which represents the NoData value of the raster band.
Oh by the way, if you want to know the informations about the data, you can simply using gdal.Info.
info = gdal.Info("pulausumba.tif")
print(info)
Next, convert the data to a numpy array and replace missing values if necessary.
# convert to a numpy array
data_array = gdal_data.ReadAsArray().astype(np.float)
data_array
# replace missing values if necessary
if np.any(data_array == nodataval):
data_array[data_array == nodataval] = np.nanAnd after that, we can visualize the data by using imshow or contourf and selecting ‘terrain’ as the colormap. See the difference of visualization between imshow and contourf here.
fig = plt.figure(figsize = (12, 8))
plt.imshow(data_array, cmap = "terrain", origin='lower')
#levels = list(range(0, 1100, 100)))
cbar = plt.colorbar()
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
As you can see, the result of the elevation map is inverted. It’s important by the way to know the region of your data, so if it turns out like it’s upside down or something you can fix it to the right position. Since I knew something was wrong on the result, I inverted it using flipud to the data array.
data_array_flip = np.flipud(data_array)
data_array_flip
And here it is, the elevation map now is on the right position.
As I mentioned earlier about the problems that came up with this library, the location coordinates don’t appear on the map so I can’t overlay it with another variable which is wind direction.
So, what is the solution?
The solution is very simple it turns out! You just need to open and read the elevation data with xarray!!! There are many dataset methods that provide by xarray library — one of them is open_rasterio. This function work with any file that rasterio can open (most often: geoTIFF). The x and y coordinates are generated automatically from the file’s geoinformation, shifted to the center of each pixel (see “PixelIsArea” Raster Space for more information). You can read here for full documentation.
So, let’s get started by open up the data with open_rasterio function.
data = xr.open_rasterio('pulausumba.tif')
data
This is the advantage of using rasterio function, the coordinates of the data will be immediately identified and we can use it while we plot on the basemap — like previous tutorials.
Then, define latitutde and longitude by taking the ‘x’ and ‘y’ columns, and raster band by selecting column 0.
lon = data_sumba_elevasi['x']
lat = data_sumba_elevasi['y']
data_band = data[0]And this is the final result of the elevation map.
fig= plt.figure(figsize=(20,8))
bm = Basemap(projection='cyl', llcrnrlon=118.75, llcrnrlat=-10.5, urcrnrlon=121, urcrnrlat=-9, resolution='i')
bm.drawcoastlines(1)
bm.drawcountries()
parallels = np.arange(-10.5, -8.9, 0.25)
bm.drawparallels(parallels, labels=[1,0,0,0], linewidth=1)
meridians = np.arange(118.75, 121.25, 0.25)
bm.drawmeridians(meridians, labels=[0,0,0,1], linewidth=1)
cf = plt.contourf(lon,lat,data_band, levels = list(range(0, 1100, 100)),cmap='terrain')
cb = plt.colorbar(cf)
cb.set_label('Elevation (m)', fontsize=15)
plt.show()
I use contourf because it can carry information of the coordinates where to plot the data, whereas imshow cannot (see the documentation of contourf here and imshow here). And of course, you can overlay with another variable after this is done because the coordinate of the data is already the same — like in the Figure 1 for example. If you are confused about how to overlay the data, you can check my topic on this one “Meteorology 101: How To Plot Wind Map”.
