Forem: Jordan Bell

OSCAR 2022 sea surface velocity streamplot animation

Jordan Bell — Sat, 17 Feb 2024 23:31:03 +0000

Sea surface currents derived from OSCAR. Dec 20, 2022

We use xarray, cartopy, and matplotlib to plot sea surface current derived from sea surface velocities. We create scrollbars for these plots using ImageMagick and combine these into an mp4 file using ffmpeg.

OSCAR (Ocean Surface Current Analysis Real-time)

OSCAR third degree resolution ocean surface currents
(OSCAR_L4_OC_third-deg). PO.DAAC

OSCAR (Ocean Surface Current Analysis Real-time) contains near-surface ocean current estimates, derived using quasi-linear and steady flow momentum equations. The horizontal velocity is directly estimated from sea surface height, surface vector wind and sea surface temperature. These data were collected from the various satellites and in situ instruments. The model formulation combines geostrophic, Ekman and Stommel shear dynamics, and a complementary term from the surface buoyancy gradient. Data are on a 1/3 degree grid with a 5 day resolution. OSCAR is generated by Earth Space Research (ESR) https://www.esr.org/research/oscar/oscar-surface-currents/. This collection contains data in 5-day files. For yearly files, see https://doi.org/10.5067/OSCAR-03D1Y

We obtain the file oscar_vel2022.nc. The format of this file is NetCDF. See NetCDF Users Guide v1.1

Xarray

Now we use Xarray to interact with the NetCDF file.

import xarray as xr

# Open the dataset
ds = xr.open_dataset("local_folder/oscar_vel2022.nc")

We inspect ds:

ds.info()

xarray.Dataset {
dimensions:
    time = 72 ;
    year = 72 ;
    depth = 1 ;
    latitude = 481 ;
    longitude = 1201 ;

variables:
    datetime64[ns] time(time) ;
        time:long_name = Day since 1992-10-05 00:00:00 ;
    float32 year(year) ;
        year:units = time in years ;
        year:long_name = Time in fractional year ;
    float32 depth(depth) ;
        depth:units = meter ;
        depth:long_name = Depth ;
    float64 latitude(latitude) ;
        latitude:units = degrees-north ;
        latitude:long_name = Latitude ;
    float64 longitude(longitude) ;
        longitude:units = degrees-east ;
        longitude:long_name = Longitude ;
    float64 u(time, depth, latitude, longitude) ;
        u:units = meter/sec ;
        u:long_name = Ocean Surface Zonal Currents ;
    float64 v(time, depth, latitude, longitude) ;
        v:units = meter/sec ;
        v:long_name = Ocean Surface Meridional Currents ;
    float64 um(time, depth, latitude, longitude) ;
        um:units = meter/sec ;
        um:long_name = Ocean Surface Zonal Currents Maximum Mask ;
    float64 vm(time, depth, latitude, longitude) ;
        vm:units = meter/sec ;
        vm:long_name = Ocean Surface Meridional Currents Maximum Mask ;

// global attributes:
    :VARIABLE = Ocean Surface Currents ;
    :DATATYPE = 1/72 YEAR Interval ;
    :DATASUBTYPE = unfiltered ;
    :GEORANGE = 20 to 420 -80 to 80 ;
    :PERIOD = Jan.01,2022 to Dec.26,2022 ;
    :year = 2022 ;
    :description = OSCAR Third Degree Sea Surface Velocity ;
    :CREATION_DATE = 02:21 06-Feb-2023 ;
    :version = 2009.0 ;
    :source = Gary Lagerloef, ESR (lager@esr.org) and Kathleen Dohan, ESR (kdohan@esr.org) ;
    :contact = Kathleen Dohan (kdohan@esr.org) or John T. Gunn (gunn@esr.org) ;
    :company = Earth & Space Research, Seattle, WA ;
    :reference = Bonjean F. and G.S.E. Lagerloef, 2002 ,Diagnostic model and analysis of the surface currents in the tropical Pacific ocean, J. Phys. Oceanogr., 32, 2,938-2,954 ;
    :note1 = Maximum Mask velocity is the geostrophic component at all points + any concurrent Ekman and buoyancy components ;
    :note2 = Longitude extends from 20 E to 420 E to avoid a break in major ocean basins. Data repeats in overlap region. ;
}

Plate carrée projection with Cartopy

Cartopy

import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import os

dec = 2 # decimation

# Create a directory for the images if it doesn't exist
output_dir = "oscar_images"
os.makedirs(output_dir, exist_ok=True)

for t in range(len(ds.time)):
    plt.figure(figsize=(18, 9))
    ax = plt.axes(projection=ccrs.PlateCarree()) # plate carrée projection

    lon = ds.longitude.values[::dec]
    lon[lon > 180] = lon[lon > 180] - 360

    plt.streamplot(
        lon,
        ds.latitude.values[::dec],
        ds.u.values[t, 0, ::dec, ::dec],
        ds.v.values[t, 0, ::dec, ::dec],
        8,
        transform = ccrs.PlateCarree()
    )

    ax.coastlines()

    # Create a title using the time dimension
    plt.title(f'Sea surface currents derived from oscar_vel2022.nc, Time: {ds.time.values[t]}')

    # Save the figure with a filename based on the time step
    plt.savefig(f"{output_dir}/oscar_vel2022_t{t}.png", dpi=150)

    # Close the figure to free up memory
    plt.close()

We now use ImageMagick convert and ffmpeg to create scrollbars and combine these into an animation.

#!/bin/bash

# Define total number of time steps
N=72

# Define output directory
output_dir="ffmpeg"

# Create output directory if it doesn't exist
mkdir -p $output_dir

# Define width and height
width=2700
height=1350

# Define progress bar height
progress_height=80

# Create N images with increasing filled progress bars
for i in $(seq 1 $N); do
    progress=$((i*width/N))  # proportionally increase size
    convert -size ${width}x${progress_height} xc:grey50 -fill "rgb(0,0,0)" -draw "rectangle 0,0 $progress,${progress_height}" ${output_dir}/progress_${i}.png
done

# Overlay progress bar onto each image, create new image
for i in $(seq 1 $N); do
    # Calculate progress bar position
    y_offset=$(( height-progress_height-10 ))  # place the progress bar at the bottom, with 10 pixels padding

    # Overlay progress bar and save as new image
    convert oscar_vel2022_t${i}.png ${output_dir}/progress_${i}.png -geometry +0+${y_offset} -composite ${output_dir}/oscar_vel2022_progress_t${i}.png
done

# Create video using ffmpeg
ffmpeg -framerate 5 -i "${output_dir}/oscar_vel2022_progress_t%d.png" -c:v libx264 -r 30 -pix_fmt yuv420p ${output_dir}/oscar_vel2022.mp4

# Clean up progress bar and composite images
rm ${output_dir}/progress_*.png
rm ${output_dir}/oscar_vel2022_progress_t*.png

We have created oscar_vel2022.mp4. This is the video file hosted on YouTube at https://youtu.be/PGpqtkXvNVw. We embed the YouTube video using the following Liquid tag.

{% youtube PGpqtkXvNVw %}

We convert this to oscar_vel2022.gif,

ffmpeg -i ffmpeg/oscar_vel2022.mp4 -vf "fps=3,scale=1000:-1:flags=lanczos" -c:v gif ffmpeg/oscar_vel2022.gif

which is this

Peters equal area projection

import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import os

dec = 2 # decimation

# Create a directory for the images if it doesn't exist
output_dir = "output_Peters"
os.makedirs(output_dir, exist_ok=True)

for t in range(len(ds.time)):
    plt.figure(figsize=(18, 9))
    ax = plt.axes(projection=ccrs.EqualEarth()) # Equal Earth (similar to Peters) projection

    lon = ds.longitude.values[::dec]
    lon[lon > 180] = lon[lon > 180] - 360

    plt.streamplot(
        lon,
        ds.latitude.values[::dec],
        ds.u.values[t, 0, ::dec, ::dec],
        ds.v.values[t, 0, ::dec, ::dec],
        8,
        transform = ccrs.PlateCarree()
    )

    ax.coastlines()

    # Create a title using the time dimension
    plt.title(f'Sea surface currents derived from oscar_vel2022.nc, Time: {ds.time.values[t]}')

    # Save the figure with a filename based on the time step
    plt.savefig(f"{output_dir}/oscar_vel2022_t{t}.png", dpi=150)

    # Close the figure to free up memory
    plt.close()

We modify these images using ImageMagick convert and combine the modified images into an animation using ffmpeg.

#!/bin/bash

# Define total number of time steps
N=72

# Define output directory
output_dir="ffmpeg_Peters"

# Create output directory if it doesn't exist
mkdir -p $output_dir

# Define width and height
width=2700
height=1350

# Define progress bar height
progress_height=80

# Create N images with increasing filled progress bars
for i in $(seq 1 $N); do
    progress=$((i*width/N))  # proportionally increase size
    convert -size ${width}x${progress_height} xc:grey50 -fill "rgb(0,0,0)" -draw "rectangle 0,0 $progress,${progress_height}" ${output_dir}/progress_${i}.png
done

# Overlay progress bar onto each image, create new image
for i in $(seq 1 $N); do
    # Calculate progress bar position
    y_offset=$(( height-progress_height-10 ))  # place the progress bar at the bottom, with 10 pixels padding

    # Overlay progress bar and save as new image
    convert oscar_vel2022_t${i}.png ${output_dir}/progress_${i}.png -geometry +0+${y_offset} -composite ${output_dir}/oscar_vel2022_progress_t${i}.png
done

# Create video using ffmpeg
ffmpeg -framerate 5 -i "${output_dir}/oscar_vel2022_progress_t%d.png" -c:v libx264 -r 30 -pix_fmt yuv420p oscar_vel2022.mp4

# Clean up progress bar and composite images
rm ${output_dir}/progress_*.png
rm ${output_dir}/oscar_vel2022_progress_t*.png

This creates a file oscar_vel2022_peters.mp4. The video file is hosted on YouTube at https://youtu.be/xy0P8l4690o.

We use ffmpeg to convert oscar_vel2022_peters.mp4 to oscar_vel2022_peters.gif:

ffmpeg -i oscar_vel2022_peters.mp4 -vf "fps=3,scale=1000:-1:flags=lanczos" -c:v gif oscar_vel2022_peters.gif

which is this:

AIS vessel density maps with pyspark and h3 and animations with ffmpeg

Jordan Bell — Fri, 16 Feb 2024 22:31:40 +0000

MarineCadastre.gov AIS vessel density maps with PySpark and h3 | by Jordan Bell | Feb, 2024 | Medium

Office for Coastal Management, 2023: Nationwide Automatic Identification System 2022, https://www.fisheries.noaa.gov/inport/item/67336

See Leveraging ESG Data to Operationalize Sustainability. November 11, 2020. Antoine Amend. Databricks

We ingest AIS message data for June 2022 into a Databricks table. (See https://github.com/jordanbell2357/uscg-nais-data/blob/main/Databricks/Databricks-AIS-01-Data-Loading.ipynb and https://github.com/jordanbell2357/uscg-nais-data/blob/main/Databricks/Databricks-AIS-10-Statistical-Tests.ipynb)

We get the following:

df_ais_june = spark.read.table("ais.june")

df_ais_june.printSchema()

root
 |-- MMSI: string (nullable = true)
 |-- BaseDateTime: timestamp (nullable = true)
 |-- LAT: double (nullable = true)
 |-- LON: double (nullable = true)
 |-- SOG: double (nullable = true)
 |-- COG: double (nullable = true)
 |-- Heading: double (nullable = true)
 |-- VesselName: string (nullable = true)
 |-- IMO: string (nullable = true)
 |-- CallSign: string (nullable = true)
 |-- VesselType: integer (nullable = true)
 |-- Status: integer (nullable = true)
 |-- Length: double (nullable = true)
 |-- Width: double (nullable = true)
 |-- Draft: double (nullable = true)
 |-- Cargo: string (nullable = true)
 |-- TranscieverClass: string (nullable = true)
 |-- filename: string (nullable = true)

In https://github.com/jordanbell2357/uscg-nais-data/blob/main/Databricks/Databricks-AIS-11-Time-Filtering.ipynb we make a plot of vessel density made using pyspark and h3 on a Databricks cluster, for June 4, 2022.

We make similar plots for each day in June in the following, also from https://github.com/jordanbell2357/uscg-nais-data/blob/main/Databricks/Databricks-AIS-11-Time-Filtering.ipynb:

import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import Polygon
import h3
import h3_pyspark
from pyspark.sql import functions as F
from pyspark.sql.functions import udf
from pyspark.sql.types import StringType
import matplotlib.colors as colors

# Define the lat/lon boundaries
min_lat, max_lat, min_lon, max_lon = 10, 60, -140, -50

# Define the hexagon resolution
resolution = 4

# Iterate over the days of the month
for day_of_month in range(1, 31):
    # Read the spark dataframe
    df_ais_june = spark.read.table("ais.june")

    # Filter the dataframe for the specific day of the month
    df_ais_june = df_ais_june.filter(F.dayofmonth('BaseDateTime') == day_of_month)

    # Add a 'resolution' column
    df_ais_june = df_ais_june.withColumn('resolution', F.lit(resolution))

    # Create an h3 index for each location using the UDF
    df_ais_june = df_ais_june.withColumn('hex_id', h3_pyspark.geo_to_h3(df_ais_june['LAT'], df_ais_june['LON'], df_ais_june['resolution']))

    # Count the occurrences of each hex_id to get the density
    df_density = df_ais_june.groupBy('hex_id').count()

    # Convert the density dataframe to Pandas
    df_density_pd = df_density.toPandas()

    # Create polygons for each unique hex_id, and filter them based on the bounding box
    hex_polygons = []
    hex_ids = []
    for hex_id in df_density_pd['hex_id'].unique():
        hex_boundary = h3.h3_to_geo_boundary(hex_id)
        polygon = Polygon([(lon, lat) for lat, lon in hex_boundary])
        # create a point based on the polygon's centroid and check if it lies within the bounding box
        centroid = polygon.centroid
        if min_lat <= centroid.y <= max_lat and min_lon <= centroid.x <= max_lon:
            hex_polygons.append(polygon)
            hex_ids.append(hex_id)

    # Filter the density dataframe to only include hexagons within the bounding box
    df_density_pd = df_density_pd[df_density_pd['hex_id'].isin(hex_ids)]

    # Replace inf with max non-inf value and NaN with 1
    df_density_pd.loc[df_density_pd['count'] == np.inf, 'count'] = df_density_pd.loc[df_density_pd['count'] != np.inf, 'count'].max()
    df_density_pd['count'] = df_density_pd['count'].replace({0:1}).fillna(1)

    # Create a GeoDataFrame
    gdf = gpd.GeoDataFrame(df_density_pd, geometry=hex_polygons)

    # Save the choropleth map to a file
    output_file = f"/dbfs/tmp/ais_june_choropleth_{day_of_month}.png"
    fig, ax = plt.subplots(1, 1, figsize=(20, 20))
    gdf.plot(column='count', cmap='YlOrRd', norm=colors.LogNorm(vmin=gdf['count'].min(), vmax=gdf['count'].max()), missing_kwds={'color': 'white'}, ax=ax)
    plt.xlim(min_lon, max_lon)
    plt.ylim(min_lat, max_lat)
    plt.xlabel('Longitude')
    plt.ylabel('Latitude')
    plt.savefig(output_file)
    plt.close(fig)

    print(f"Choropleth map for day {day_of_month} saved to: {output_file}")

Choropleth map for day 1 saved to: /dbfs/tmp/ais_june_choropleth_1.png
Choropleth map for day 2 saved to: /dbfs/tmp/ais_june_choropleth_2.png
Choropleth map for day 3 saved to: /dbfs/tmp/ais_june_choropleth_3.png
Choropleth map for day 4 saved to: /dbfs/tmp/ais_june_choropleth_4.png
Choropleth map for day 5 saved to: /dbfs/tmp/ais_june_choropleth_5.png
Choropleth map for day 6 saved to: /dbfs/tmp/ais_june_choropleth_6.png
Choropleth map for day 7 saved to: /dbfs/tmp/ais_june_choropleth_7.png
Choropleth map for day 8 saved to: /dbfs/tmp/ais_june_choropleth_8.png
Choropleth map for day 9 saved to: /dbfs/tmp/ais_june_choropleth_9.png
Choropleth map for day 10 saved to: /dbfs/tmp/ais_june_choropleth_10.png
Choropleth map for day 11 saved to: /dbfs/tmp/ais_june_choropleth_11.png
Choropleth map for day 12 saved to: /dbfs/tmp/ais_june_choropleth_12.png
Choropleth map for day 13 saved to: /dbfs/tmp/ais_june_choropleth_13.png
Choropleth map for day 14 saved to: /dbfs/tmp/ais_june_choropleth_14.png
Choropleth map for day 15 saved to: /dbfs/tmp/ais_june_choropleth_15.png
Choropleth map for day 16 saved to: /dbfs/tmp/ais_june_choropleth_16.png
Choropleth map for day 17 saved to: /dbfs/tmp/ais_june_choropleth_17.png
Choropleth map for day 18 saved to: /dbfs/tmp/ais_june_choropleth_18.png
Choropleth map for day 19 saved to: /dbfs/tmp/ais_june_choropleth_19.png
Choropleth map for day 20 saved to: /dbfs/tmp/ais_june_choropleth_20.png
Choropleth map for day 21 saved to: /dbfs/tmp/ais_june_choropleth_21.png
Choropleth map for day 22 saved to: /dbfs/tmp/ais_june_choropleth_22.png
Choropleth map for day 23 saved to: /dbfs/tmp/ais_june_choropleth_23.png
Choropleth map for day 24 saved to: /dbfs/tmp/ais_june_choropleth_24.png
Choropleth map for day 25 saved to: /dbfs/tmp/ais_june_choropleth_25.png
Choropleth map for day 26 saved to: /dbfs/tmp/ais_june_choropleth_26.png
Choropleth map for day 27 saved to: /dbfs/tmp/ais_june_choropleth_27.png
Choropleth map for day 28 saved to: /dbfs/tmp/ais_june_choropleth_28.png
Choropleth map for day 29 saved to: /dbfs/tmp/ais_june_choropleth_29.png
Choropleth map for day 30 saved to: /dbfs/tmp/ais_june_choropleth_30.png

We then download these daily plot files locally. Next we use ffmpeg. ffmpeg Documentation

The following bash script uses ImageMagick convert to make daily progress bars and overlay those with the daily plots. ImageMagick convert documentation

#!/bin/bash
# ffmpeg.sh

# Create 30 images with increasing filled progress bars
for i in {1..30}; do
    progress=$((i*1000/30))  # increase size 10 times
    convert -size 1000x200 xc:grey50 -fill "rgb(100,149,237)" -draw "rectangle 0,0 $progress,200" progress_$(printf "%02d" $i).png
done

width=1440  # width of original images

# Overlay progress bar onto each image, create new image
for i in {1..30}; do
    # Calculate progress bar position
    x_offset=$(( (width-1000)/2 ))  # center the progress bar

    # Overlay progress bar and save as new image
    convert ais_june_choropleth_$i.png progress_$(printf "%02d" $i).png -geometry +$x_offset+10 -composite ais_june_choropleth_progress_$(printf "%02d" $i).png
done

# Create video using ffmpeg
ffmpeg -framerate 5 -pattern_type glob -i 'ais_june_choropleth_progress_*.png' -c:v libx264 -r 30 -pix_fmt yuv420p ais_june_choropleth.mp4

# Clean up progress bar and composite images
rm progress_*.png
rm ais_june_choropleth_progress_*.png

This creates ais_june_choropleth.mp4. This video file is hosted on YouTube:

We remark in passing that the Liquid syntax used for embedding the above YouTube video is

{% youtube Rt0RTpmlaK0 %}

Finally, we use ffmpeg to convert ais_june_choropleth.mp4 to ais_june_choropleth.gif.

ffmpeg -i ais_june_choropleth.mp4 -vf "fps=3,scale=1000:-1:flags=lanczos" -c:v gif ais_june_choropleth.gif

The file ais_june_choropleth.gif is shown below.

Configuring and using Hadoop and Spark on Ubuntu 22.04 LTS (with Canada 2021 Census data)

Jordan Bell — Wed, 14 Feb 2024 21:02:18 +0000

Hadoop and Spark on Ubuntu 22.04 LTS with Canada 2021 Census data
| by Jordan Bell | Feb, 2024 | Medium

Introduction

This is a single guide that explains every step needed to configure and use Hadoop and Spark, and uses it with a geographically granular census dataset.

This tutorial works for Ubuntu 22.04 LTS both as VM (tested using VirtualBox) and using WSL2.

Install and configure dependencies

Reference: https://hadoop.apache.org/docs/stable/hadoop-project-dist/hadoop-common/SingleCluster.html

Secondary references: https://learning.oreilly.com/library/view/hadoop-the-definitive/9781491901687/apa.html#I_sect1_d1e37797 and https://kontext.tech/article/448/install-hadoop-330-on-linux

sudo apt-get install default-jre
sudo apt-get install openjdk-11-jre-headless
sudo apt-get install openjdk-8-jre-headless
sudo apt-get install openjdk-8-jdk

sudo apt-get install openssh-client openssh-server
sudo apt-get install ssh

ssh-keygen -t rsa -P '' -f ~/.ssh/id_rsa
cat ~/.ssh/id_rsa.pub >> ~/.ssh/authorized_keys
chmod 0600 ~/.ssh/authorized_keys

Download and configure Hadoop

https://hadoop.apache.org/releases.html

cd ~
wget https://dlcdn.apache.org/hadoop/common/hadoop-3.3.6/hadoop-3.3.6.tar.gz
tar -xvzf hadoop-3.3.6.tar.gz

ls ~/hadoop-3.3.6/etc/hadoop

capacity-scheduler.xml            httpfs-env.sh               mapred-site.xml
configuration.xsl                 httpfs-log4j.properties     shellprofile.d
container-executor.cfg            httpfs-site.xml             ssl-client.xml.example
core-site.xml                     kms-acls.xml                ssl-server.xml.example
hadoop-env.cmd                    kms-env.sh                  user_ec_policies.xml.template
hadoop-env.sh                     kms-log4j.properties        workers
hadoop-metrics2.properties        kms-site.xml                yarn-env.cmd
hadoop-policy.xml                 log4j.properties            yarn-env.sh
hadoop-user-functions.sh.example  mapred-env.cmd              yarn-site.xml
hdfs-rbf-site.xml                 mapred-env.sh               yarnservice-log4j.properties
hdfs-site.xml                     mapred-queues.xml.template

core-site.xml

<configuration>
    <property>
        <name>fs.defaultFS</name>
        <value>hdfs://localhost:9000</value>
    </property>
</configuration>

hdfs-site.xml

<configuration>
    <property>
        <name>dfs.replication</name>
        <value>1</value>
    </property>
</configuration>

mapred-site.xml

<configuration>
    <property>
        <name>mapreduce.framework.name</name>
        <value>yarn</value>
    </property>
    <property>
        <name>mapreduce.application.classpath</name>
        <value>$HADOOP_MAPRED_HOME/share/hadoop/mapreduce/*:$HADOOP_MAPRED_HOME/share/hadoop/mapreduce/lib/*</value>
    </property>
</configuration>

yarn-site.xml

<configuration>
    <property>
        <name>yarn.nodemanager.aux-services</name>
        <value>mapreduce_shuffle</value>
    </property>
    <property>
        <name>yarn.nodemanager.env-whitelist</name>
        <value>JAVA_HOME,HADOOP_COMMON_HOME,HADOOP_HDFS_HOME,HADOOP_CONF_DIR,CLASSPATH_PREPEND_DISTCACHE,HADOOP_YARN_HOME,HADOOP_HOME,PATH,LANG,TZ,HADOOP_MAPRED_HOME</value>
    </property>
</configuration>

Reference: https://stackoverflow.com/questions/29260900/could-not-find-or-load-main-class-com-sun-tools-javac-main-hadoop-mapreduce

Add to environment variables; either for single user with ~/.bashrc or for system with /etc/environment

export JAVA_HOME=/usr/lib/jvm/java-1.8.0-openjdk-amd64
export HADOOP_HOME=~/hadoop-3.3.6
export PATH=$PATH:$HADOOP_HOME/bin:$HADOOP_HOME/sbin
export HADOOP_CONF_DIR=$HADOOP_HOME/etc/hadoop
export HADOOP_CLASSPATH=/usr/lib/jvm/java-8-openjdk-amd64/lib/tools.jar

hdfs namenode -format

start-all.sh

For WSL2 or VM:

ip address

1: lo: <LOOPBACK,UP,LOWER_UP> mtu 65536 qdisc noqueue state UNKNOWN group default qlen 1000
    link/loopback 00:00:00:00:00:00 brd 00:00:00:00:00:00
    inet 127.0.0.1/8 scope host lo
       valid_lft forever preferred_lft forever
    inet6 ::1/128 scope host
       valid_lft forever preferred_lft forever
2: eth0: <BROADCAST,MULTICAST,UP,LOWER_UP> mtu 1500 qdisc mq state UP group default qlen 1000
    link/ether 00:15:5d:a6:9a:40 brd ff:ff:ff:ff:ff:ff
    inet 172.17.194.43/20 brd 172.17.207.255 scope global eth0
       valid_lft forever preferred_lft forever
    inet6 fe80::215:5dff:fea6:9a40/64 scope link
       valid_lft forever preferred_lft forever

In Windows, viewing http://172.17.194.43:9870/ or http://localhost:9870/ yields

Likewise, we can view Yarn on http://172.17.194.43:8088/ or http://localhost:8088/

CSV dataset

Canada 2021 Census:

https://www12.statcan.gc.ca/census-recensement/2021/dp-pd/prof/details/download-telecharger.cfm

We use the most geographically granular data, for Aggregate Dissemination Areas (ADA).

https://www.gnu.org/software/wget/manual/wget.html

‘ --content-disposition’
This option is useful for some file-downloading CGI programs that use Content-Disposition headers to describe what the name of a downloaded file should be.

wget --content-disposition "https://www12.statcan.gc.ca/census-recensement/2021/dp-pd/prof/details/download-telecharger/comp/GetFile.cfm?Lang=E&FILETYPE=CSV&GEONO=012"
unzip 98-401-X2021012_eng_CSV.zip

Archive:  98-401-X2021012_eng_CSV.zip
  inflating: 98-401-X2021012_English_CSV_data.csv
  inflating: 98-401-X2021012_English_meta.txt
  inflating: 98-401-X2021012_Geo_starting_row.CSV
  inflating: README_meta.txt

head 98-401-X2021012_English_CSV_data.csv

CENSUS_YEAR,DGUID,ALT_GEO_CODE,GEO_LEVEL,GEO_NAME,TNR_SF,TNR_LF,DATA_QUALITY_FLAG,CHARACTERISTIC_ID,CHARACTERISTIC_NAME,CHARACTERISTIC_NOTE,C1_COUNT_TOTAL,SYMBOL,C2_COUNT_MEN+,SYMBOL,C3_COUNT_WOMEN+,SYMBOL,C10_RATE_TOTAL,SYMBOL,C11_RATE_MEN+,SYMBOL,C12_RATE_WOMEN+,SYMBOL
2021,"2021S051610010001","10010001","Aggregate dissemination area","10010001",3.1,4.6,"00000",1,"Population, 2021",1,8881,"",,"...",,"...",,"...",,"...",,"..."
2021,"2021S051610010001","10010001","Aggregate dissemination area","10010001",3.1,4.6,"00000",2,"Population, 2016",1,9334,"",,"...",,"...",,"...",,"...",,"..."
2021,"2021S051610010001","10010001","Aggregate dissemination area","10010001",3.1,4.6,"00000",3,"Population percentage change, 2016 to 2021",,-4.9,"",,"...",,"...",-4.9,"",,"...",,"..."
2021,"2021S051610010001","10010001","Aggregate dissemination area","10010001",3.1,4.6,"00000",4,"Total private dwellings",2,5737,"",,"...",,"...",,"...",,"...",,"..."
2021,"2021S051610010001","10010001","Aggregate dissemination area","10010001",3.1,4.6,"00000",5,"Private dwellings occupied by usual residents",3,4121,"",,"...",,"...",,"...",,"...",,"..."
2021,"2021S051610010001","10010001","Aggregate dissemination area","10010001",3.1,4.6,"00000",6,"Population density per square kilometre",,9.4,"",,"...",,"...",9.4,"",,"...",,"..."
2021,"2021S051610010001","10010001","Aggregate dissemination area","10010001",3.1,4.6,"00000",7,"Land area in square kilometres",,941.33,"",,"...",,"...",,"...",,"...",,"..."
2021,"2021S051610010001","10010001","Aggregate dissemination area","10010001",3.1,4.6,"00000",8,"Total - Age groups of the population - 100% data",,8880,"",4365,"",4515,"",100,"",100,"",100,""
2021,"2021S051610010001","10010001","Aggregate dissemination area","10010001",3.1,4.6,"00000",9,"  0 to 14 years",,925,"",460,"",465,"",10.4,"",10.5,"",10.3,""

head -n 16 98-401-X2021012_English_meta.txt

** General Information

Census Profile, 2021 - Aggregate Dissemination Areas
Catalogue number: 98-401-X2021012 Release date: 2022-12-15
Note on data quality and disclosure control
For information on data quality, as well as random rounding, data suppression, and other methods of disclosure control, please consult the Guide to the Census of Population, 2021, Catalogue no. 98-304-X, and the reference guides for the individual topics.
Source: Statistics Canada, 2021 Census, Catalogue no. 98-316-X2021001.
How to cite: Statistics Canada. 2022. Census Profile. 2021 Census. Statistics Canada Catalogue no. 98-316-X2021001. Ottawa. Released December 15, 2022. https://www12.statcan.gc.ca/census-recensement/2021/dp-pd/prof/index.cfm?Lang=E

Also from 98-401-X2021012_English_meta.txt:

Aggregate dissemination area (ADA)

An aggregate dissemination area (ADA) is a dissemination geography created for the Census. ADAs cover the entire country and, where possible, have a population between 5,000 and 15,000 based on the previous census population counts. ADAs are created by grouping existing dissemination geographic areas, including census tracts (CTs), census subdivisions (CSDs) or dissemination areas (DAs). ADA boundaries respect provincial, territorial, census division (CD), census metropolitan area (CMA) and census agglomeration (CA) boundaries.

The intent of the ADA geography is to ensure the availability of census data, where possible, across all regions of Canada.

Putting CSV into HDFS

Reference: https://hadoop.apache.org/docs/stable/hadoop-project-dist/hadoop-common/FileSystemShell.html

cd ~
hadoop fs -mkdir /census2021
hadoop fs -mkdir /census2021/ada
hadoop fs -put ./98-401-X2021012_English_CSV_data.csv /census2021/ada/ada.csv

hadoop fs -ls -R /census2021

drwxr-xr-x   - ubuntu supergroup          0 2024-02-09 09:50 /census2021/ada
-rw-r--r--   1 ubuntu supergroup 2415200889 2024-02-09 09:50 /census2021/ada/ada.csv

hadoop fs -head /census2021/ada/ada.csv

CENSUS_YEAR,DGUID,ALT_GEO_CODE,GEO_LEVEL,GEO_NAME,TNR_SF,TNR_LF,DATA_QUALITY_FLAG,CHARACTERISTIC_ID,CHARACTERISTIC_NAME,CHARACTERISTIC_NOTE,C1_COUNT_TOTAL,SYMBOL,C2_COUNT_MEN+,SYMBOL,C3_COUNT_WOMEN+,SYMBOL,C10_RATE_TOTAL,SYMBOL,C11_RATE_MEN+,SYMBOL,C12_RATE_WOMEN+,SYMBOL
2021,"2021S051610010001","10010001","Aggregate dissemination area","10010001",3.1,4.6,"00000",1,"Population, 2021",1,8881,"",,"...",,"...",,"...",,"...",,"..."
2021,"2021S051610010001","10010001","Aggregate dissemination area","10010001",3.1,4.6,"00000",2,"Population, 2016",1,9334,"",,"...",,"...",,"...",,"...",,"..."
2021,"2021S051610010001","10010001","Aggregate dissemination area","10010001",3.1,4.6,"00000",3,"Population percentage change, 2016 to 2021",,-4.9,"",,"...",,"...",-4.9,"",,"...",,"..."
2021,"2021S051610010001","10010001","Aggregate dissemination area","10010001",3.1,4.6,"00000",4,"Total private dwellings",2,5737,"",,"...",,"...",,"...",,"...",,"..."
2021,"2021S051610010001","10010001","Aggregate dissemination area","100100

hadoop fs -cat /census2021/ada/ada.csv | wc -l

14294224

(We remind ourselves that hadoop fs -head displays the first kilobyte of a file, not a fixed number of lines.)

MapReduce code for line count

From https://statinfer.com/301-2-3-map-reduce-code-for-line-count/

cf. https://hadoop.apache.org/docs/stable/hadoop-mapreduce-client/hadoop-mapreduce-client-core/MapReduceTutorial.html

We make new file /home/ubuntu/hadoop-3.3.6/bin/LineCount.java with the following contents:

import java.io.IOException;
import org.apache.hadoop.io.IntWritable;
import org.apache.hadoop.io.LongWritable;
import org.apache.hadoop.io.Text;
import org.apache.hadoop.conf.Configuration;
import org.apache.hadoop.conf.Configured;
import org.apache.hadoop.fs.FileSystem;
import org.apache.hadoop.fs.Path;
import org.apache.hadoop.io.IntWritable;
import org.apache.hadoop.mapreduce.lib.input.TextInputFormat;
import org.apache.hadoop.mapreduce.lib.output.TextOutputFormat;
import org.apache.hadoop.util.Tool;
import org.apache.hadoop.util.ToolRunner;
import org.apache.hadoop.mapreduce.Job;
import org.apache.hadoop.mapreduce.Mapper;
import org.apache.hadoop.mapreduce.Reducer;
import org.apache.hadoop.mapreduce.lib.input.FileInputFormat;
import org.apache.hadoop.mapreduce.lib.output.FileOutputFormat;

public class LineCount{
    public static class LineCntMapper extends Mapper<LongWritable, Text, Text, IntWritable> {
    Text keyEmit = new Text("Total Lines");
    private final static IntWritable one = new IntWritable(1);

    public void map(LongWritable key, Text value, Context context){
        try {
            context.write(keyEmit, one);
        }
        catch (IOException e) {
            e.printStackTrace();
            System.exit(0);
        }
        catch (InterruptedException e) {
            e.printStackTrace();
            System.exit(0);
        }
    }
}

    public static class LineCntReducer extends Reducer<Text, IntWritable, Text, IntWritable> {
    public void reduce(Text key, Iterable<IntWritable> values, Context context){
        int sum = 0;
        for (IntWritable val : values) {
            sum += val.get();
        }
        try {
            context.write(key, new IntWritable(sum));
        }
        catch (IOException e) {
            e.printStackTrace();
            System.exit(0);
        }
        catch (InterruptedException e) {
            e.printStackTrace();
            System.exit(0);
        }
    }
}

    public static void main(String[] args) throws Exception {
    Configuration conf = new Configuration();
    Job job = Job.getInstance(conf, "line count2");
    job.setJarByClass(LineCount.class);
    job.setMapperClass(LineCntMapper.class);
    job.setCombinerClass(LineCntReducer.class);
    job.setReducerClass(LineCntReducer.class);
    job.setOutputKeyClass(Text.class);
    job.setOutputValueClass(IntWritable.class);
    FileInputFormat.addInputPath(job, new Path(args[0]));
    FileOutputFormat.setOutputPath(job, new Path(args[1]));
    System.exit(job.waitForCompletion(true) ? 0 : 1);
    }
}

cd /home/ubuntu/hadoop-3.3.6/bin/
hadoop com.sun.tools.javac.Main LineCount.java
jar cf lc.jar LineCount*.class
hadoop jar lc.jar LineCount /census2021/ada/ada.csv /LineCount

2024-02-09 09:53:23,449 INFO client.DefaultNoHARMFailoverProxyProvider: Connecting to ResourceManager at /0.0.0.0:8032
2024-02-09 09:53:23,793 WARN mapreduce.JobResourceUploader: Hadoop command-line option parsing not performed. Implement the Tool interface and execute your application with ToolRunner to remedy this.
2024-02-09 09:53:23,804 INFO mapreduce.JobResourceUploader: Disabling Erasure Coding for path: /tmp/hadoop-yarn/staging/ubuntu/.staging/job_1707010689166_0003
2024-02-09 09:53:23,975 INFO input.FileInputFormat: Total input files to process : 1
2024-02-09 09:53:24,022 INFO mapreduce.JobSubmitter: number of splits:18
2024-02-09 09:53:24,127 INFO mapreduce.JobSubmitter: Submitting tokens for job: job_1707010689166_0003
2024-02-09 09:53:24,127 INFO mapreduce.JobSubmitter: Executing with tokens: []
2024-02-09 09:53:24,248 INFO conf.Configuration: resource-types.xml not found
2024-02-09 09:53:24,248 INFO resource.ResourceUtils: Unable to find 'resource-types.xml'.
2024-02-09 09:53:24,593 INFO impl.YarnClientImpl: Submitted application application_1707010689166_0003
2024-02-09 09:53:24,625 INFO mapreduce.Job: The url to track the job: http://LAPTOP-JBell.:8088/proxy/application_1707010689166_0003/
2024-02-09 09:53:24,626 INFO mapreduce.Job: Running job: job_1707010689166_0003
2024-02-09 09:53:29,730 INFO mapreduce.Job: Job job_1707010689166_0003 running in uber mode : false
2024-02-09 09:53:29,731 INFO mapreduce.Job:  map 0% reduce 0%
2024-02-09 09:54:17,390 INFO mapreduce.Job:  map 33% reduce 0%
2024-02-09 09:54:24,509 INFO mapreduce.Job:  map 67% reduce 0%
2024-02-09 09:54:31,629 INFO mapreduce.Job:  map 94% reduce 0%
2024-02-09 09:54:36,652 INFO mapreduce.Job:  map 100% reduce 100%
2024-02-09 09:54:37,663 INFO mapreduce.Job: Job job_1707010689166_0003 completed successfully
2024-02-09 09:54:37,776 INFO mapreduce.Job: Counters: 54
        File System Counters
                FILE: Number of bytes read=330
                FILE: Number of bytes written=5240478
                FILE: Number of read operations=0
                FILE: Number of large read operations=0
                FILE: Number of write operations=0
                HDFS: Number of bytes read=2415272483
                HDFS: Number of bytes written=21
                HDFS: Number of read operations=59
                HDFS: Number of large read operations=0
                HDFS: Number of write operations=2
                HDFS: Number of bytes read erasure-coded=0
        Job Counters
                Launched map tasks=18
                Launched reduce tasks=1
                Data-local map tasks=18
                Total time spent by all maps in occupied slots (ms)=346566
                Total time spent by all reduces in occupied slots (ms)=11939
                Total time spent by all map tasks (ms)=346566
                Total time spent by all reduce tasks (ms)=11939
                Total vcore-milliseconds taken by all map tasks=346566
                Total vcore-milliseconds taken by all reduce tasks=11939
                Total megabyte-milliseconds taken by all map tasks=354883584
                Total megabyte-milliseconds taken by all reduce tasks=12225536
        Map-Reduce Framework
                Map input records=14294224
                Map output records=14294224
                Map output bytes=228707584
                Map output materialized bytes=432
                Input split bytes=1962
                Combine input records=14294224
                Combine output records=18
                Reduce input groups=1
                Reduce shuffle bytes=432
                Reduce input records=18
                Reduce output records=1
                Spilled Records=36
                Shuffled Maps =18
                Failed Shuffles=0
                Merged Map outputs=18
                GC time elapsed (ms)=50184
                CPU time spent (ms)=88410
                Physical memory (bytes) snapshot=8216698880
                Virtual memory (bytes) snapshot=49281294336
                Total committed heap usage (bytes)=8874622976
                Peak Map Physical memory (bytes)=475123712
                Peak Map Virtual memory (bytes)=2596765696
                Peak Reduce Physical memory (bytes)=242405376
                Peak Reduce Virtual memory (bytes)=2599333888
        Shuffle Errors
                BAD_ID=0
                CONNECTION=0
                IO_ERROR=0
                WRONG_LENGTH=0
                WRONG_MAP=0
                WRONG_REDUCE=0
        File Input Format Counters
                Bytes Read=2415270521
        File Output Format Counters
                Bytes Written=21

hadoop fs -ls -R /LineCount

-rw-r--r--   1 ubuntu supergroup          0 2024-02-09 09:54 /LineCount/_SUCCESS
-rw-r--r--   1 ubuntu supergroup         21 2024-02-09 09:54 /LineCount/part-r-00000

hadoop fs -cat /LineCount/part-r-00000

Total Lines 14294224

Total lines are 14294224. This includes the header, thus there are 14294223 records.

Spark

References: https://spark.apache.org/docs/latest/spark-standalone.html and https://spark.apache.org/docs/3.3.4/configuration.html

https://spark.apache.org/downloads.html

cd ~
wget https://dlcdn.apache.org/spark/spark-3.3.4/spark-3.3.4-bin-hadoop3.tgz
tar -xvzf spark-3.3.4-bin-hadoop3.tgz

export SPARK_HOME=/home/ubuntu/spark-3.3.4-bin-hadoop3
export PATH=$SPARK_HOME/bin:$PATH
export SPARK_DIST_CLASSPATH=$(hadoop classpath)

cp $SPARK_HOME/conf/spark-defaults.conf.template $SPARK_HOME/conf/spark-defaults.conf

Include line

spark.driver.host localhost

run-example SparkPi 10

24/02/08 17:42:53 INFO DAGScheduler: Job 0 finished: reduce at SparkPi.scala:38, took 0.770038 s
Pi is roughly 3.1433631433631435
24/02/08 17:42:53 INFO SparkUI: Stopped Spark web UI at http://localhost:4040

spark-shell

SLF4J: Class path contains multiple SLF4J bindings.
SLF4J: Found binding in [jar:file:/home/ubuntu/spark-3.3.4-bin-hadoop3/jars/log4j-slf4j-impl-2.17.2.jar!/org/slf4j/impl/StaticLoggerBinder.class]
SLF4J: Found binding in [jar:file:/home/ubuntu/hadoop-3.3.6/share/hadoop/common/lib/slf4j-reload4j-1.7.36.jar!/org/slf4j/impl/StaticLoggerBinder.class]
SLF4J: See http://www.slf4j.org/codes.html#multiple_bindings for an explanation.
SLF4J: Actual binding is of type [org.apache.logging.slf4j.Log4jLoggerFactory]
24/02/09 16:54:19 WARN Utils: Your hostname, LAPTOP-JBell resolves to a loopback address: 127.0.1.1; using 172.17.194.43 instead (on interface eth0)
24/02/09 16:54:19 WARN Utils: Set SPARK_LOCAL_IP if you need to bind to another address
Setting default log level to "WARN".
To adjust logging level use sc.setLogLevel(newLevel). For SparkR, use setLogLevel(newLevel).
24/02/09 16:54:24 WARN NativeCodeLoader: Unable to load native-hadoop library for your platform... using builtin-java classes where applicable
Spark context Web UI available at http://localhost:4040
Spark context available as 'sc' (master = local[*], app id = local-1707515665178).
Spark session available as 'spark'.
Welcome to
      ____              __
     / __/__  ___ _____/ /__
    _\ \/ _ \/ _ `/ __/  '_/
   /___/ .__/\_,_/_/ /_/\_\   version 3.3.4
      /_/

Using Scala version 2.12.15 (OpenJDK 64-Bit Server VM, Java 1.8.0_392)
Type in expressions to have them evaluated.
Type :help for more information.

scala> val df = spark.read.option("header", "true").csv("/census2021/ada/ada.csv")
df: org.apache.spark.sql.DataFrame = [CENSUS_YEAR: string, DGUID: string ... 21 more fields]

scala> val distinctGeoNameCount = df.select("GEO_NAME").distinct().count()
distinctGeoNameCount: Long = 5433

That is, the number of ADA is 5433.

At the end of our tutorial, we calculate the number of ADA another way.

For WSL2, like the above with Hadoop, we can view the Spark Web UI using http://172.17.194.43:4040/ and http://localhost:4040/

spark-shell

SLF4J: Class path contains multiple SLF4J bindings.
SLF4J: Found binding in [jar:file:/home/ubuntu/spark-3.3.4-bin-hadoop3/jars/log4j-slf4j-impl-2.17.2.jar!/org/slf4j/impl/StaticLoggerBinder.class]
SLF4J: Found binding in [jar:file:/home/ubuntu/hadoop-3.3.6/share/hadoop/common/lib/slf4j-reload4j-1.7.36.jar!/org/slf4j/impl/StaticLoggerBinder.class]
SLF4J: See http://www.slf4j.org/codes.html#multiple_bindings for an explanation.
SLF4J: Actual binding is of type [org.apache.logging.slf4j.Log4jLoggerFactory]
24/02/09 09:58:14 WARN Utils: Your hostname, LAPTOP-JBell resolves to a loopback address: 127.0.1.1; using 172.17.194.43 instead (on interface eth0)
24/02/09 09:58:14 WARN Utils: Set SPARK_LOCAL_IP if you need to bind to another address
Setting default log level to "WARN".
To adjust logging level use sc.setLogLevel(newLevel). For SparkR, use setLogLevel(newLevel).
24/02/09 09:58:19 WARN NativeCodeLoader: Unable to load native-hadoop library for your platform... using builtin-java classes where applicable
Spark context Web UI available at http://localhost:4040
Spark context available as 'sc' (master = local[*], app id = local-1707490699967).
Spark session available as 'spark'.
Welcome to
      ____              __
     / __/__  ___ _____/ /__
    _\ \/ _ \/ _ `/ __/  '_/
   /___/ .__/\_,_/_/ /_/\_\   version 3.3.4
      /_/

Using Scala version 2.12.15 (OpenJDK 64-Bit Server VM, Java 1.8.0_392)
Type in expressions to have them evaluated.
Type :help for more information.

scala> val df = spark.read.option("header", "true").csv("/census2021/ada/ada.csv")
df: org.apache.spark.sql.DataFrame = [CENSUS_YEAR: string, DGUID: string ... 21 more fields]

scala> val rowCount = df.count()
rowCount: Long = 14294223

scala> df.columns.foreach(println)
CENSUS_YEAR
DGUID
ALT_GEO_CODE
GEO_LEVEL
GEO_NAME
TNR_SF
TNR_LF
DATA_QUALITY_FLAG
CHARACTERISTIC_ID
CHARACTERISTIC_NAME
CHARACTERISTIC_NOTE
C1_COUNT_TOTAL
SYMBOL12
C2_COUNT_MEN+
SYMBOL14
C3_COUNT_WOMEN+
SYMBOL16
C10_RATE_TOTAL
SYMBOL18
C11_RATE_MEN+
SYMBOL20
C12_RATE_WOMEN+
SYMBOL22

scala> val df = spark.read.option("header", "true").csv("/census2021/ada/ada.csv")
df: org.apache.spark.sql.DataFrame = [CENSUS_YEAR: string, DGUID: string ... 21 more fields]

scala> val dfRenamed = df.withColumnRenamed("C2_COUNT_MEN+", "C2_COUNT_MEN").withColumnRenamed("C3_COUNT_WOMEN+", "C3_COUNT_WOMEN").withColumnRenamed("C11_RATE_MEN+", "C11_RATE_MEN").withColumnRenamed("C12_RATE_WOMEN+", "C12_RATE_WOMEN")
dfRenamed: org.apache.spark.sql.DataFrame = [CENSUS_YEAR: string, DGUID: string ... 21 more fields]

scala> dfRenamed.columns.foreach(println)
CENSUS_YEAR
DGUID
ALT_GEO_CODE
GEO_LEVEL
GEO_NAME
TNR_SF
TNR_LF
DATA_QUALITY_FLAG
CHARACTERISTIC_ID
CHARACTERISTIC_NAME
CHARACTERISTIC_NOTE
C1_COUNT_TOTAL
SYMBOL12
C2_COUNT_MEN
SYMBOL14
C3_COUNT_WOMEN
SYMBOL16
C10_RATE_TOTAL
SYMBOL18
C11_RATE_MEN
SYMBOL20
C12_RATE_WOMEN
SYMBOL22

scala> val df = spark.read.option("header", "true").csv("/census2021/ada/ada.csv")
df: org.apache.spark.sql.DataFrame = [CENSUS_YEAR: string, DGUID: string ... 21 more fields]

scala> val dfRenamed = df.withColumnRenamed("C2_COUNT_MEN+", "C2_COUNT_MEN").withColumnRenamed("C3_COUNT_WOMEN+", "C3_COUNT_WOMEN").withColumnRenamed("C11_RATE_MEN+", "C11_RATE_MEN").withColumnRenamed("C12_RATE_WOMEN+", "C12_RATE_WOMEN")
dfRenamed: org.apache.spark.sql.DataFrame = [CENSUS_YEAR: string, DGUID: string ... 21 more fields]

scala> val df_reduced = dfRenamed.select("GEO_NAME", "CHARACTERISTIC_ID", "CHARACTERISTIC_NAME", "C1_COUNT_TOTAL", "C2_COUNT_MEN", "C3_COUNT_WOMEN", "C10_RATE_TOTAL", "C11_RATE_MEN", "C12_RATE_WOMEN")
df_reduced: org.apache.spark.sql.DataFrame = [GEO_NAME: string, CHARACTERISTIC_ID: string ... 7 more fields]

scala> df_reduced.coalesce(1).write.option("header", "true").csv("/census2021/ada_reduced")
[Stage 1:>                                                          (0 + 1) / 1]

hadoop fs -ls -R /census2021

drwxr-xr-x   - ubuntu supergroup          0 2024-02-09 10:29 /census2021/ada
-rw-r--r--   1 ubuntu supergroup 2415200889 2024-02-09 09:50 /census2021/ada/ada.csv
drwxr-xr-x   - ubuntu supergroup          0 2024-02-09 10:32 /census2021/ada_reduced
-rw-r--r--   1 ubuntu supergroup          0 2024-02-09 10:32 /census2021/ada_reduced/_SUCCESS
-rw-r--r--   1 ubuntu supergroup  763087510 2024-02-09 10:32 /census2021/ada_reduced/part-00000-f87a37db-597e-42df-bd05-bab7766f0c41-c000.csv

hadoop fs -head /census2021/ada_reduced/part-00000-f87a37db-597e-42df-bd05-bab7766f0c41-c000.csv

GEO_NAME,CHARACTERISTIC_ID,CHARACTERISTIC_NAME,C1_COUNT_TOTAL,C2_COUNT_MEN,C3_COUNT_WOMEN,C10_RATE_TOTAL,C11_RATE_MEN,C12_RATE_WOMEN
10010001,1,"Population, 2021",8881,,,,,
10010001,2,"Population, 2016",9334,,,,,
10010001,3,"Population percentage change, 2016 to 2021",-4.9,,,-4.9,,
10010001,4,Total private dwellings,5737,,,,,
10010001,5,Private dwellings occupied by usual residents,4121,,,,,
10010001,6,Population density per square kilometre,9.4,,,9.4,,
10010001,7,Land area in square kilometres,941.33,,,,,
10010001,8,Total - Age groups of the population - 100% data,8880,4365,4515,100,100,100
10010001,9,0 to 14 years,925,460,465,10.4,10.5,10.3
10010001,10,0 to 4 years,265,130,135,3,3,3
10010001,11,5 to 9 years,295,135,150,3.3,3.1,3.3
10010001,12,10 to 14 years,370,190,180,4.2,4.4,4
10010001,13,15 to 64 years,5010,2500,2510,56.4,57.3,55.6
10010001,14,15 to 19 years,370,185,180,4.2,4.2,4
10010001,15,20 to 24 years,280,155,120,3.2,3.6,2.7
10010001,16,25 to 29 years,250,130,120,2.8,3,2.7
10010001,17,30 to 34 years

(We remind ourselves that hadoop fs -head displays the first kilobyte of a file, not a fixed number of lines.)

scala> val df = spark.read.option("header", "true").csv("/census2021/ada_reduced/part-00000-f87a37db-597e-42df-bd05-bab7766f0c41-c000.csv")
df: org.apache.spark.sql.DataFrame = [GEO_NAME: string, CHARACTERISTIC_ID: string ... 7 more fields]

scala> val df_id = df.select($"CHARACTERISTIC_ID".cast("integer").alias("CHARACTERISTIC_ID"), $"CHARACTERISTIC_NAME").distinct().orderBy("CHARACTERISTIC_ID")
df_id: org.apache.spark.sql.Dataset[org.apache.spark.sql.Row] = [CHARACTERISTIC_ID: int, CHARACTERISTIC_NAME: string]

scala> df_id.coalesce(1).write.option("header", "true").csv("/census2021/ada_characteristic")

hadoop fs -ls -R /census2021/ada_characteristic

-rw-r--r--   1 ubuntu supergroup          0 2024-02-09 13:31 /census2021/ada_characteristic/_SUCCESS
-rw-r--r--   1 ubuntu supergroup      78032 2024-02-09 13:31 /census2021/ada_characteristic/part-00000-3912cf78-b4fe-4986-9e70-9a62cff91ab6-c000.csv

hadoop fs -head /census2021/ada_characteristic/part-00000-3912cf78-b4fe-4986-9e70-9a62cff91ab6-c000.csv

CHARACTERISTIC_ID,CHARACTERISTIC_NAME
1,"Population, 2021"
2,"Population, 2016"
3,"Population percentage change, 2016 to 2021"
4,Total private dwellings
5,Private dwellings occupied by usual residents
6,Population density per square kilometre
7,Land area in square kilometres
8,Total - Age groups of the population - 100% data
9,0 to 14 years
10,0 to 4 years
11,5 to 9 years
12,10 to 14 years
13,15 to 64 years
14,15 to 19 years
15,20 to 24 years
16,25 to 29 years
17,30 to 34 years
18,35 to 39 years
19,40 to 44 years
20,45 to 49 years
21,50 to 54 years
22,55 to 59 years
23,60 to 64 years
24,65 years and over
25,65 to 69 years
26,70 to 74 years
27,75 to 79 years
28,80 to 84 years
29,85 years and over
30,85 to 89 years
31,90 to 94 years
32,95 to 99 years
33,100 years and over
34,Total - Distribution (%) of the population by broad age groups - 100% data
35,0 to 14 years
36,15 to 64 years
37,65 years and over
38,85 years and over
39,Average age of the population
40,Median age of the population
41,Total - Occupied pri

We can determine the number of CHARACTERISTIC_ID:

hadoop fs -cat /census2021/ada_characteristic/part-00000-3912cf78-b4fe-4986-9e70-9a62c
ff91ab6-c000.csv | wc -l

2632

Thus there are 2631 CHARACTERISTIC_ID.

There are 14294223 records in the datasets. Thus there are ¹⁴²⁹⁴²²³⁄₂₆₃₁ = 5433 ADA.

gpsdecode and ogr2ogr for AIS messages

Jordan Bell — Wed, 14 Feb 2024 19:00:35 +0000

Access to AIS data

Kystverket (Norwegian Coastal Administration)

AIS data are also available at the following IP address: 153.44.253.27 port 5631

The AIS data stream is in standard IEC format IEC 62320-1.

See https://www.confluent.io/blog/streaming-etl-and-analytics-for-real-time-location-tracking/

OpenCPN

We start by viewing the AIS data from 153.44.253.27 port 5631 using OpenCPN.

AIS and NMEA message format

https://documentation.spire.com/tcp-stream-v2/

https://www.itu.int/rec/R-REC-M.1371-5-201402-I/en

AIVDM/AIVDO protocol decoding. Eric S. Raymond. Version 1.58, 24 June 2023.

ncat

Ncat Users’ Guide. Nmap

Ncat is a feature-packed networking utility which reads and writes data across networks from the command line. Ncat was written for the Nmap Project as a much-improved reimplementation of the venerable Netcat. It uses both TCP and UDP for communication and is designed to be a reliable back-end tool to instantly provide network connectivity to other applications and users. Ncat will not only work with IPv4 and IPv6 but provides the user with a virtually limitless number of potential uses.

timeout 3600s ncat 153.44.253.27 5631 > ais_2024_02_14_01_15_UTC

The messages look like this:

head -n 5 ais_2024_02_14_01_15_UTC

\s:2573555,c:1707873398*05\!BSVDM,1,1,,A,13mB7d001E1dDH0`KlKBjB?<0490,0*34
\s:2573105,c:1707873398*04\!B2VDM,2,1,4,B,53oHa402;QilhU1D00085=@v0TT000000000000t6`p>440Ht:j3lU1CcCCl,0*7A
\s:2573105,c:1707873398*04\!B2VDM,2,2,4,B,lh000000000,2*5F
\s:2573415,c:1707873398*00\!BSVDM,1,1,,B,13m9ERPP00Pp5opUivG8P?w:0<12,0*6F
\s:2573210,c:1707873398*03\!BSVDM,1,1,,B,13oQg60000PJMu`QhUB6D26v0D18,0*40

Inspecting format of messages with xxd

We use xxd and sed.

ncat 153.44.253.27 5631 | xxd -p -c 1024 | sed 's/0d0a/&\n/g'

Uses ncat to receive the data.
Pipes the data into xxd -p -c 1024 to convert it into a plain hex dump. The -c 1024 sets the number of columns to 1024, ensuring that each message up to 1024 bytes is processed as a single unit.
Pipes to sed 's/0d0a/&\n\n/g' to search for 0d0a (the hex representation of \r\n) in the hex output and replace it with itself followed by a newline.

ncat 153.44.253.27 5631 | xxd -p -c 1024 | sed 's/0d0a/&\n/g'

5c733a323537333438352c633a313730383138383437352a30465c21425356444d2c312c312c2c412c42315777743f5030354843733d3939756565416b57777634335030362c302a36460d0a
21425356444d2c312c312c2c422c31336d486751303030305049554f6a5264305250493251563044303b2c302a34330d0a
5c733a323537333333352c633a313730383138383437352a30335c21425356444d2c312c312c2c412c31336d3a6e6e3730303c30614256465477736d623d706d68303850652c302a34380d0a
5c733a323537333134352c633a313730383138383437352a30365c21425356444d2c312c312c2c412c31334d43306b4030303d303a446f4c504e44733864563f68303850652c302a33310d0a
5c733a323537333532352c633a313730383138383437352a30345c21425356444d2c312c312c2c412c48315777743f50606e3038754556333e32323232323167373e456c2c322a35360d0a

This is a fun roundabout way of seeing that the messages are terminated with \r\n. We could use ncat -x to directly output the hex dump, but I didn't figure out a way nicely to parse the separate messages.

gpsdecode

gpsdecode tool is a batch-mode decoder for NMEA and various binary packet formats associated with GPS, AIS, and differential-correction services. It produces a JSON dump on standard output from binary on standard input. The JSON is the same format documented by gpsd; this tool uses the same decoding logic as gpsd, but with a simpler interface intended for batch processing of data files.

All sensor-input formats known to the GPSD project can be decoded by this tool. These include: NMEA, AIVDM (the NMEA-derived sentence format used by AIS, the marine Automatic Identification System), RTCM2, and all supported GPS binary formats (notably including SiRF). See gpsd(8) for applicable standards and known limitations of the decoding logic.

You can use this tool with nc(1) to examine AIS feeds from AIS pooling services, RTCM feeds from RTCM receivers or NTRIP broadcasters.

The messages parsed using gpsdecode look like this:

head -n 5 ais_2024_02_14_01_15_UTC | gpsdecode

{"class":"AIS","device":"stdin","type":1,"repeat":0,"mmsi":257198000,"scaled":true,"status":0,"status_text":"Under way using engine","turn":0,"speed":8.5,"accuracy":false,"lon":23.662507,"lat":70.664715,"course":71.3,"heading":71,"second":38,"maneuver":0,"raim":false,"radio":16960}
{"class":"AIS","device":"stdin","type":1,"repeat":0,"mmsi":257054090,"scaled":true,"status":0,"status_text":"Under way using engine","turn":"nan","speed":0.0,"accuracy":true,"lon":12.253433,"lat":66.026820,"course":217.6,"heading":511,"second":37,"maneuver":0,"raim":false,"radio":49218}
{"class":"AIS","device":"stdin","type":1,"repeat":0,"mmsi":259551000,"scaled":true,"status":0,"status_text":"Under way using engine","turn":0,"speed":0.0,"accuracy":true,"lon":5.782060,"lat":58.998307,"course":161.6,"heading":67,"second":31,"maneuver":0,"raim":false,"radio":81992}

jq

cat ais_2024_02_14_01_15_UTC_decoded_with_timestamps| jq --slurp '.' > ais_2024_02_14_01_15_UTC_decoded_with_timestamps.json

jq '.[0]' ais_2024_02_14_01_15_UTC_decoded_with_timestamps.json

{
  "class": "AIS",
  "device": "stdin",
  "type": 1,
  "repeat": 0,
  "mmsi": 257198000,
  "scaled": true,
  "status": 0,
  "status_text": "Under way using engine",
  "turn": 0,
  "speed": 8.5,
  "accuracy": false,
  "lon": 23.662507,
  "lat": 70.664715,
  "course": 71.3,
  "heading": 71,
  "second": 38,
  "maneuver": 0,
  "raim": false,
  "radio": 16960,
  "timestamp": ""
}

Longitude and latitude

head -n 5 ais_2024_02_14_01_15_UTC |
gpsdecode |
jq -r '. | select(.lon != null and .lat != null) | "\(.lon) \(.lat)"'

23.662507 70.664715
12.253433 66.02682
5.78206 58.998307

cat ais_2024_02_14_01_15_UTC |
gpsdecode |
jq -r '. | select(.lon != null and .lat != null) | "\(.lon) \(.lat)"' > \
ais_2024_02_14_01_15_UTC.dat

It is worth noting that | does not require \ to split across lines while > does.

Bounding box with awk

awk '$1 >= 0 && $1 <= 35 && $2 >= 55 && $2 <= 80' ais_2024_02_14_01_15_UTC.dat > ais_2024_02_14_01_15_UTC_box.dat

wc -l ais_2024_02_14_01_15_UTC.dat
72999 ais_2024_02_14_01_15_UTC.dat
wc -l ais_2024_02_14_01_15_UTC_box.dat
72585 ais_2024_02_14_01_15_UTC_box.dat

Gnuplot

gnuplot> plot 'ais_2024_02_14_01_15_UTC_box.dat' using 1:2 with points

GDAL ogr2ogr

echo "Longitude,Latitude" > ais_2024_02_14_01_15_UTC_box.csv

awk '{print $1 "," $2}' ais_2024_02_14_01_15_UTC_box.dat >> ais_2024_02_14_01_15_UTC_box.csv

head -n 5 ais_2024_02_14_01_15_UTC_box.csv

Longitude,Latitude
23.662507,70.664715
12.253433,66.02682
5.78206,58.998307
2.966662,62.485987

ogr2ogr

ogr2ogr -f "CSV" -lco GEOMETRY=AS_XY -t_srs "EPSG:32633" -s_srs "EPSG:4326" projected_output/ ais_2024_02_14_01_15_UTC_box.csv -oo X_POSSIBLE_NAMES=Longitude -oo Y_POSSIBLE_NAMES=Latitude -overwrite

head -n 5 projected_output/ais_2024_02_14_01_15_UTC_box.csv
X,Y,Longitude,Latitude

819151.812467396,7862821.57947704,23.662507,70.664715
375505.625008375,7325626.69347061,12.253433,66.02682
-28538.7387103818,6576427.56565985,5.78206,58.998307
-117577.293012299,6986142.22730701,2.966662,62.485987

plot_ais.gp

set terminal pngcairo enhanced
set output 'ais_plot.png'
set title 'Ship Positions in Norwegian Waters'
set xlabel 'UTM Easting (kilometers)'
set ylabel 'UTM Northing (kilometers)'
set grid
set datafile separator ','
plot 'projected_output/ais_2024_02_14_01_15_UTC_box.csv' using ($1/1000):($2/1000) with points pt 7 ps 1 lc rgb 'blue' title 'Ships'

gnuplot plot_ais.gp

Including ISO 8601 timestamps

https://www.w3.org/TR/NOTE-datetime

ISO 8601 UTC timestamp format YYYY-MM-DDThh:mm:ssZ

To pipe the output of ncat, we do the following:

ncat 153.44.253.27 5631 | tr -d '\r'

Now we can compose this pipe with gpsdecode:

ncat 153.44.253.27 5631 | tr -d '\r' | gpsdecode

{"class":"AIS","device":"stdin","type":3,"repeat":0,"mmsi":257878000,"scaled":true,"status":0,"status_text":"Under way using engine","turn":"fastright","speed":6.0,"accuracy":false,"lon":9.639210,"lat":59.107997,"course":278.4,"heading":284,"second":30,"maneuver":0,"raim":false,"radio":7792}
{"class":"AIS","device":"stdin","type":1,"repeat":0,"mmsi":257021970,"scaled":true,"status":7,"status_text":"Engaged in fishing","turn":"nan","speed":0.0,"accuracy":false,"lon":17.902577,"lat":69.509332,"course":360.0,"heading":511,"second":32,"maneuver":0,"raim":false,"radio":27224}
{"class":"AIS","device":"stdin","type":1,"repeat":0,"mmsi":235103057,"scaled":true,"status":0,"status_text":"Under way using engine","turn":0,"speed":11.3,"accuracy":false,"lon":5.716548,"lat":63.704728,"course":206.7,"heading":204,"second":30,"maneuver":0,"raim":false,"radio":27224}

Now we compose with awk and date (https://www.gnu.org/software/coreutils/manual/html_node/Examples-of-date.html):

ncat 153.44.253.27 5631 | tr -d '\r' | gpsdecode | awk -v date="$(date -u +%Y-%m-%dT%H:%M:%SZ)" '{print "{\"timestamp_ISO_8601\":\"" date "\"," substr($0, 2)}'

{"timestamp_ISO_8601":"2024-02-16T21:38:38Z","class":"AIS","device":"stdin","type":18,"repeat":0,"mmsi":257129980,"scaled":true,"reserved":0,"speed":0.0,"accuracy":true,"lon":23.329870,"lat":69.977463,"course":360.0,"heading":511,"second":39,"regional":0,"cs":false,"display":false,"dsc":true,"band":true,"msg22":true,"raim":true,"radio":609914}
{"timestamp_ISO_8601":"2024-02-16T21:38:38Z","class":"AIS","device":"stdin","type":1,"repeat":0,"mmsi":257482600,"scaled":true,"status":0,"status_text":"Under way using engine","turn":"nan","speed":0.1,"accuracy":false,"lon":5.209840,"lat":59.298867,"course":27.4,"heading":511,"second":39,"maneuver":0,"raim":false,"radio":49221}
{"timestamp_ISO_8601":"2024-02-16T21:38:38Z","class":"AIS","device":"stdin","type":5,"repeat":0,"mmsi":258033200,"scaled":true,"imo":0,"ais_version":1,"callsign":"JWCH","shipname":"NORDNES","shiptype":0,"shiptype_text":"Not available","to_bow":9,"to_stern":6,"to_port":2,"to_starboard":2,"epfd":1,"epfd_text":"GPS","eta":"02-14T20:00Z","draught":0.9,"destination":"NAVIGATION TRAINING","dte":0}

timeout 3600s ncat 153.44.253.27 5631 | tr -d '\r' | gpsdecode | awk -v date="$(date -u +%Y-%m-%dT%H:%M:%SZ)" '{print "{\"timestamp_ISO_8601\":\"" date "\"," substr($0, 2)}' > ais_2024_02_16T21_42_57Z

We examine the contents of the file we created above:

head -n 2 ais_2024_02_16T21_42_57Z

{"timestamp_ISO_8601":"2024-02-16T21:43:35Z","class":"AIS","device":"stdin","type":1,"repeat":0,"mmsi":257827700,"scaled":true,"status":15,"status_text":"Not defined","turn":0,"speed":0.0,"accuracy":false,"lon":14.376283,"lat":67.284670,"course":342.6,"heading":48,"second":35,"maneuver":0,"raim":false,"radio":81943}
{"timestamp_ISO_8601":"2024-02-16T21:43:35Z","class":"AIS","device":"stdin","type":21,"repeat":1,"mmsi":992576156,"scaled":true,"aid_type":30,"aid_type_text":"Special Mark","name":"AQUACULTURE 4","lon":18.480670,"lat":69.895130,"accuracy":false,"to_bow":1,"to_stern":1,"to_port":1,"to_starboard":1,"epfd":7,"epfd_text":"Surveyed","second":60,"regional":0,"off_position":false,"raim":false,"virtual_aid":true}

Because the interrupt from timeout 3600s did not coincide with a complete message from ncat, the tail of the file looks like:

tail -n 1 ais_2024_02_16T21_42_57Z

{"timestamp_ISO_8601":"2024-02-16T21:43:35Z","class":"AIS","devi

Thus we omit the last line of the file, and then parse as JSON using jq --slurp:

head -n -1 ais_2024_02_16T21_42_57Z | jq --slurp '.' > ais_2024_02_16T21_42_57Z.json

We inspect the first element:

jq '.[0]' ais_2024_02_16T21_42_57Z.json

{
  "timestamp_ISO_8601": "2024-02-16T21:43:35Z",
  "class": "AIS",
  "device": "stdin",
  "type": 1,
  "repeat": 0,
  "mmsi": 257827700,
  "scaled": true,
  "status": 15,
  "status_text": "Not defined",
  "turn": 0,
  "speed": 0,
  "accuracy": false,
  "lon": 14.376283,
  "lat": 67.28467,
  "course": 342.6,
  "heading": 48,
  "second": 35,
  "maneuver": 0,
  "raim": false,
  "radio": 81943
}

Now we use jq to convert this to the GeoJSON format (RFC 7946).

jq '{type: "FeatureCollection", features: map({type: "Feature", geometry: {type: "Point", coordinates: [.lon, .lat]}, properties: .})}' ais_2024_02_16T21_42_57Z.json > ais_2024_02_16T21_42_57Z.geojson

ais_2024_02_16T21_42_57Z.geojson on ArcGIS Online

Map view on ArcGIS Online

ais_2024_02_16T21_42_57Z.geojson on GitHub

Ais Decoder by Neal Arundale

Ais Decoder Options:

Ais Decoder main menu:

Ais Decoder Nmea input:

Ais Decoder parsed Summary:

Ais Decoder NMEA sentence Detail:

Compare to gpsdecode:

echo '!BSVDM,1,1,,B,13mL5v0001PEjC6Rsu:Pca`L0@9W,0*59,2/14/2024 6:36:17 PM' | gpsdecode | jq --slurp '.' | jq '.[0]'

{
  "class": "AIS",
  "device": "stdin",
  "type": 1,
  "repeat": 0,
  "mmsi": 257361400,
  "scaled": true,
  "status": 0,
  "status_text": "Under way using engine",
  "turn": 0,
  "speed": 0.1,
  "accuracy": true,
  "lon": 4.759205,
  "lat": 61.056497,
  "course": 17.4,
  "heading": 308,
  "second": 14,
  "maneuver": 0,
  "raim": false,
  "radio": 66151
}

Sea surface currents derived from OSCAR

Jordan Bell — Tue, 20 Dec 2022 00:25:23 +0000

import xarray as xr
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

with xr.open_dataset('oscar_vel2020.nc.gz') as ds:
    print(ds)

    plt.figure(figsize=(18, 9))

    ax = plt.axes(projection=ccrs.PlateCarree()) # plate carrée projection

    dec = 10

    lon = ds.longitude.values[::dec]

    lon[lon > 180] = lon[lon > 180] - 360

    mymap = plt.streamplot(
        lon,
        ds.latitude.values[::dec],
        ds.u.values[0, 0, ::dec, ::dec],
        ds.v.values[0, 0, ::dec, ::dec],
        8,
        transform = ccrs.PlateCarree()
    )

    ax.coastlines()

    plt.title('Sea surface currents derived from OSCAR')

    plt.savefig("currents18x9.png", dpi=150)

    plt.show()

Lightly modified from "How to plot Ocean Currents with Cartopy". Reply by Jody Klymak. Oct 20, 2019

Ocean Surface Current Analyses Real-time (OSCAR)

NASA Earth Science Data Systems (ESDS)