Categories
Reanalysis

Downloading MERRA2 data

I tried to download some MERRA2 data, and realised that I’d deleted all my old scripts! So, here I go again…

Preliminary steps:

  1. Register / Log in to the NASA EARTHDATA site: https://urs.earthdata.nasa.gov/home
  2. Follow the steps to link GES DISC with your account: https://disc.gsfc.nasa.gov/earthdata-login
  3. Verify by downloading this example data file URL

Go to the NASA Goddard Earth Sciences (GES) Data and Information Services Center (DISC) page: https://disc.gsfc.nasa.gov/

  • Source: Models/Analyses MERRA-2
  • Temporal Resolution: 1 hour

Choose “Subset / Get data” for the instantaneous ‘Single-Level Diagnostics’ file.

Under “Download Method” there are two relevant options, as we want to download a spatial domain:

  • Subset using OPeNDAP
  • Subset using GES DISC. This option seems to also allow regridding.

I selected OPeNDAP, as I don’t need to regrid the data.

I didn’t change the dates under ‘Refine Date Range’: 1980-01-01 to 2021-08-31

I used ‘Refine Region’ to select a region around Ireland: -11,51,-5,56 and I kept the ‘Use ‘Refine Region’ box checked.

Under ‘Variables’ I selected surface wind speed, and I kept output format as netCDF.

Click ‘Get Data’ and download the links list

Follow the steps to set up wget or curl here: https://disc.gsfc.nasa.gov/data-access#mac_linux_wget

I downloaded the data to the mccdata server using curl.

Categories
Uncategorized

New European Wind Atlas

The new European wind atlas is an great tool to download and visualise high-resolution and publicly available data of wind resources specifically suitable for wind energy development, site prospecting and wind farm design.

The mesoscale modelling covers the entire EU plus Turkey and 100 km offshore as well as the complete North and Baltic Seas. The WRF model was used in a NEWA consortium developed configuration, with a 3 km grid spacing and simulation period covering 30 years (1989-2018). The temporal resolution is 30 minutes and there are eight heights available from 10m, 50m, 75m up to 500m. Parameters available include wind speed, wind direction, air temperature, air density, shortwave radiation, among many others.

Currently there are data from 2009-2018 available to download from the website. Unfortunately the project ended before they were able to convert all of the data to the format used for distribution on the website. They are currently working on finding some resources to finish this conversion and they hope to have the remaining years available around the end of this year.

Categories
Uncategorized

Taylor Diagrams and high resolution plots of land use and elevation.

This is an overview for producing these following plots.

Taylor Diagram

The Taylor diagram graphically illustrates how well a model represents the variability of the observations it is attempting to predict, in terms of correlation, standard deviation and RMSE.

This code was adapted from the python file in this GitHub repository to suit our needs, credit goes to them.

This involves the writing of two python functions, one for generating circles on a graph and one that pulls all the data together, draws the appropriate circles and produces the plot. These are shown in the code blocks below.

def drawCircle(radius, xcenter=0, ycenter=0, maxradius=np.inf,style="--", color="k", lw=1, label=""):
    angles = np.linspace(-np.pi / 2, np.pi / 2, 361)
    x = np.sin(angles) * radius + xcenter
    y = np.cos(angles) * radius + ycenter

  # Only keep points within the circle
    I = np.where(x ** 2 + y ** 2 < maxradius ** 2)[0]
    if(len(I) == 0):
        return
    x = x[I]
    y = y[I]
    plt.plot(x, y, style, color=color, lw=lw, zorder=-100, label=label)
    plt.plot(x, -y, style, color=color, lw=lw, zorder=-100)
    
    
def makeTaylor(corrs,stds,stdobs,circ_LIM=1.0,loc='location'):
    '''This function takes a list of correlations and standard deviations 
       corresponding to 6 different model outputs in this case (should alter 
       the code for your application), along with the standard deviation
       of the observations. circ_LIM determines the radius of the outer
       circle in the plot and loc is a name identifier (like station name)
       for plot titles and filenames for saving plots
    '''


    plot_colors = ['blue','purple','brown','orange','green','red']
    ang = np.arccos(corrs)
    x = stds * np.cos(ang)
    y = stds * np.sin(ang)
    
    fig, axs = plt.subplots(1, 1,figsize=(8,8))

    plt.scatter(x[:3], y[:3], marker='o',s=120, color=plot_colors[:3],edgecolor='k', label=plot_cols)
    plt.scatter(x[3:], y[3:], marker='X',s=180, color=plot_colors[3:],edgecolor='k', label=plot_cols)

    plt.ylim(0,circ_LIM+0.05)
    plt.xlim(0,circ_LIM+0.05)
    xticks = np.arange(0,circ_LIM,0.1)

    drawCircle(stdobs, xcenter=0, ycenter=0, style="--", color="cyan", lw=3)

    plt.xticks(xticks[xticks >= 0])
    plt.text(np.sin(np.pi / 4) * circ_LIM, np.cos(np.pi / 4) * circ_LIM,
        "Correlation", rotation=-45,horizontalalignment="center", verticalalignment="bottom")
    plt.gca().yaxis.set_visible(False)
    plt.gca().xaxis.set_ticks_position('bottom')


    plt.plot(stdobs, 0, 's-', color='cyan', label="Observation", mew=2, clip_on=False)


    corrs = [0, 0.5, 0.8, 0.9, 0.95,
                0.99]
    for i in range(0, len(corrs)):
        angl = np.arccos(corrs[i]) 
        xl = np.cos(angl) * circ_LIM
        yl = np.sin(angl) * circ_LIM
        plt.plot([0, xl], [0, yl], 'k--')
        plt.text(xl, yl, str(corrs[i]), verticalalignment="bottom")

    xticks = plt.xticks()[0]

    Rs = np.linspace(0.1, 1.0, 10)
    for R in Rs:
        if(R > 0):
            drawCircle(R, xcenter=stdobs, ycenter=0, maxradius=circ_LIM, style="-", color="gray", lw=3)
            xc = np.sin(-np.pi / 4) * R + stdobs
            yc = np.cos(np.pi / 4) * R
            if(xc >0):
                plt.text(xc, yc, str(R)[:3], horizontalalignment="right",
                     verticalalignment="bottom",
                     color="gray")

    for X in plt.xticks()[0]:
        if(X <= circ_LIM):
            drawCircle(X, style=":")
    drawCircle(circ_LIM, style="-", lw=3)

    plt.text(0.55,0.85,loc,fontsize=48)

    plt.tight_layout()

    plt.savefig('plots/'+loc.replace(' ','_')+'_TaylorWind',dpi=300)
    plt.close(fig)

This will produce a plot like this:

Corine Land use data

The Corine land use dataset can be downloaded from the EPA website. This data comes in the form of a .shp file which I had not encountered before. A quick google suggested that the library geopandas would allow me to analyse the data. The data are comprised of 18882 rows × 8 columns. The important column is “Class_Desc” which gives the names of the land use classifiers.

import geopandas as gpd
import matplotlib.pyplot as plt

cr = gpd.read_file('/CLC18_IE_ITM/CLC18_IE_ITM.shp')

cr.plot(column='Class_Desc',cmap='gist_earth_r')
plt.savefig('corine_ALL_Ireland',dpi=600)

This produced a map of the ROI shown below.

The colouring represents each of the 35 land use categories in the dataset.

Looking into the data further, it becomes apparent that the data is mapped in the Irish Transverse Mercator projection on an easting and northing grid. It is possible to zoom in on certain areas by indexing the data appropriately, and you can get lat-lon coordinates converted using this handy tool from OSI. An example below zooming in on Kerry and dropping pins (black crosses) at two location in Lerrig and Dromnacurra in the north of the county.

import geopandas as gpd
import matplotlib.pyplot as plt

cr = gpd.read_file('CLC18_IE_ITM/CLC18_IE_ITM.shp')

fig, ax = plt.subplots(1, 1,figsize=(10,10))

cr.plot(column='Class_Desc',cmap='gist_earth_r',ax=ax)

#Markers, with values converted from OSI tool
plt.plot(480561.6873,624590.0764,color='k',mec='w',marker='X',markersize=12)
plt.plot(478445.6175,633327.0275,color='k',mec='w',marker='X',markersize=12)

#Arbitrary plot limits
plt.xlim(435000,525500)
plt.ylim(582000,650000)

plt.savefig('Kerry_LU',dpi=400)

Digital Elevation Map (DEM)

One option to download high resolution elevation data is outlined here in this handy tutorial. You will need to install the following libraries.

$ conda create -n gdal_test
$ conda activate gdal_test
$ conda install gdal
$ conda install elevation

These allow you to both download and plot DEM data. The follow should be entered in the command line to download your data (obviously change your bounding box to suit your application). I have chosen to focus on Kerry (again).

$ !eio clip -o Kerry-30m-DEM.tif --bounds -10.6 51. -8.9 53. 

This will download the .tif that you have named in your command above. Now you can use gdal within python to visualise the data.

from osgeo import gdal
import numpy as np
import matplotlib 
import matplotlib.pyplot as plt

filename = "Kerry-30m-DEM.tif"
gdal_data = gdal.Open(filename)
gdal_band = gdal_data.GetRasterBand(1)
nodataval = gdal_band.GetNoDataValue()

# 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.nan

#Plot out data with Matplotlib's 'contour'
fig = plt.figure(figsize = (12, 8))
ax = fig.add_subplot(111)

## I reversed the data below using [-1::-1] because the data was 
## initially plotting upside down! This fixed it.
plt.contourf(data_array[-1::-1,:], cmap = "viridis", 
            levels = list(range(0, 1050, 50)))

plt.title("Elevation Contours over Kerry")
cbar = plt.colorbar()
plt.savefig('Kerry_DEM',dpi=400)
plt.close(fig)

Kerry does not have the same terrain complexity as Mt Shasta in the tutorial, so remember to adjust your contouring levels accordingly.

These data are mapped on a “World Geodetic System 1984” projection.

Categories
Uncategorized

Using conda install on a server

conda is great, it makes installing lots of packages easy. However, if you have an account on a server which already has conda installed, you should use conda to only install packages for your own user, not system-wide. To do this, you should make use of conda environments.

For example, let’s say you want to use cdo on your server. First, create a conda environment, I’ve called it cdo-env:

conda create -y -n cdo-env

Now activate your environment, so you can work within it:

conda activate cdo-env

When you’re in the environment, the shell prompt may include the name of the environment, as a handy reminder:

(cdo-env) [sweeneyc@login ~]$

Now you can install cdo within your cdo-env environment:

conda install -c conda-forge cdo

If all goes well, this should install cdo within your cdo-env enviroment, and you can use it whenever you like!

To close/exit your conda environment, use the command:

conda deactivate

That’s it!