Airplanes above Mainz

Quiet a while ago I bought a SDR-RTL USB dongle. Originally to capture the signal from my outdoor temperature sensor. However, no one had decipherd the signal for that device so far and I was not motivated enough to invest too much time in it. I almost forgot the dongle and remembered it recently, when I read about Flight radar and dump1090, which logs the ID, position and elevation signal of airplanes. Almost all airplanes transmit this unencrypted signal about themselves at 1090 MHz the ‘Automatic Dependent Surveillance-Broadcast System (ADS-B)’. And I was curious how I could log and visualize those flight tracks from Frankfurt/Main airport towards Mainz.

I’ve build an R package to process the data from dump1090 and present some function for flat and 3D visualization for it.

Continue reading “Airplanes above Mainz”

Advertisements

Accelerate NetCDF-4 read access

For simulating the (global) terrestrial biosphere I need time series of climate data, which is mostly provided in a 3-dimesional NetCDF file with the dimensions longitude, latitude and time. If everything is in a single file, this is often several GB large and extracting a timeseries along one position can take quiet some time, since by default NetCDF files are optimized to extract time slices. In the following I show how I speed up reading from a large NetCDF file along the time axis.

Continue reading “Accelerate NetCDF-4 read access”

Flying bumblebees

Spring arrived and pollinators started their work. In 2013 I took a picture of a bumblebee, just leaving a flower and flying towards me by chance. If it was a sharp picture it would have been really great picture. Since then the same happened several times and in the following I show those photos, which are at least so (un-)sharp, that the insect can be recognized.

IMG_0682_01

Continue reading “Flying bumblebees”

Twitter friends and (un-)followers with R

A few weeks ago  the number of my twitter followers exceeded 100 and was going slightly down again and I was interested in who unfollowed me. Twitter does not offer that and I would have had to register at a third party page, which even wanted write access to my twitter account.  I did not think about it for a while and then I found the R package twitteR, which lets you fetch all required information and I used it in the following to evaluate my “friends” and “(un-)followers”.
Continue reading “Twitter friends and (un-)followers with R”

Broken axis with ggplot2

For visualizing my data I use R and the library ggplot2. And just lately I made some sensitivity simulations with out dynamic global vegetation model (DGVM) LPJ-GUESS. While summarizing the data per ecosystem and having a first look at the data I realized, that one ecosystem has up to 10 times higher values than all the others. That made me searching for “broken axis” and I didn’t find a satisfying solution, so I had to create my own.

Continue reading “Broken axis with ggplot2”

Wasp variety in our community garden

Last summer I realized that there is quiet a variety of wasps is our community garden. And I know that not only bees are pollenizer, but I was not aware of wasps doing the same job. I took several pictures of them, of which I now publish the nicest here.

UPDATED SPECIES NAMES (20.5.2016) Thanks to Alexandra Stevens (BUND Mainz) who contacted an expert.

Continue reading “Wasp variety in our community garden”

Mapping traceroutes with Perl

I was playing around on my smartphone, looking for interesting apps and found one that displays the trace of a request to a server (e.g. www.google.com). Nevertheless, since I didn’t really need it, it was to expensive and I thought: “ok, let’s do it myself!” and I wrote a little perl script and uploaded it to GitHub. The script traces the route of a request and either produces a static google map as png or writes a GPX file, which can be used in GIS-software or Google maps.

Continue reading “Mapping traceroutes with Perl”

Retrieving and reprojecting MODIS satellite products

Data from the MODIS satellites are a great source for comparison to vegetation models and also for teaching. However, they come in a format I am not really used to and in a strange projection. The format is HDF, of which even several versions exists. But there are also several tools to convert it to a “better” file format for me. HDF needs much less space compared to other formats. Therefore it makes sense for satellite data to be stored as HDF.

First I need to know which data I want and in which region. Just recently a colleague asked me to get the fire product (MCD45A1) for Vietnam for him. Therefore I show it as example here.

On the map provided on the MODIS overview webpage one can try to select the correct row (v) and column (h) values. However, I already said it is a strange projection. I reprojected the grid to a standard equirectangular projection:

modis sinusoidal

That’s also not easy to select the right h/v values. However with both maps it worked. I set up a perl script using my module WebService::MODIS, which is also available on CPAN. Here is the example script, how to download the desired satellite data:

#!/usr/bin/env perl
use warnings;
use strict;
use WebService::MODIS;
readCache; # or "initCache;" if not done yet
my $fire = WebService::MODIS->new();
$fire->product("MCD45A1");
$fire->version('051');
$fire->h([27,28]);
$fire->v([6,7,8]);
$fire->dates(["2007.01.01", "2014.12.31"]);
$fire->createUrl;
#print "$_\n" foreach $fire->url;
$fire->download();

Now the files are in the current working directory. I do the mosaicing and reprojecting with ModisDownload.R. You may now ask: “Why not also downloading the data with the R Skript?” I did it in the beginning and it was very frustrating, since it crashed or did not start downloading at all and continued download was not possible. Therefore the perl module. To make use of different directories, where the downloaded files were and the reprojected files should be saved, I modified the original ModisDownload.R. Here is the diff of line 184 and 185 (Version 3.2, 13th Jan. 2014):

- write(paste(getwd(),"/",hdfNames[1], sep=""), mosaicname)
- for (j in 2:length(hdfNames)) write(paste(getwd(),"/",hdfNames[j], sep=""),mosaicname,append=T)
+ write(hdfNames[1], mosaicname)
+ for (j in 2:length(hdfNames)) write(hdfNames[j],mosaicname,append=T)

Now the filenames given to the mosaicHDF function must contain the full path to the files if called from another directory. ModisDownload.R uses the Modis Reprojection Tool (MRT), which is written in Java. Last year we used it for teaching. However, the installation under Windows it is a pain in the … Therefore, we preprocessed the data this year and supplied the students with GeoTiffs. Here is now the R script, which processes the downloaded hdf files and creates one GeoTiff per timeslice in UTM projection in the current working directory:

MRT.PATH <- "/path/to/MRT/bin"
MODISDOWNLOAD.DIR <- "/path/to/ModisDownload.R"
source(paste(MODISDOWNLOAD.DIR, "ModisDownload.R", sep="/"))
MODIS <- list(product="MCD45A1", bands="1 0 0 0 0 0 0 0", pixel.size=500)
DOWNLOAD.DIR <- "/path/to/downloaded/MODIS/files"
OUTPUT.DIR <- "/path/to/large/disk/for/GeoTIFFs"
setwd(OUTPUT.DIR)
downloaded <- list.files(DOWNLOAD.DIR, pattern="hdf$", full.names=TRUE)
if (length(downloaded) != 0) {
  yyyyddd <- sub("^.*\\.A", "", downloaded)
  yyyyddd <- unique(sub("\\..*$", "", yyyyddd))
  year <- trunc(as.integer(yyyyddd) / 1000)
  day.of.year <- as.integer(yyyyddd) - year * 1000
  dates <- as.Date(paste(year,"-01-01", sep=""), format="%Y-%m-%d") + day.of.year - 1
  for (j in 1:length(dates)) {
    in.hdf <- downloaded[grepl(paste("A", yyyyddd[j], sep=""), downloaded)]
    out.hdf <- paste(MODIS[["product"]], "_", dates[j], ".hdf", sep="")
    out.tif <- paste(MODIS[["product"]], "_", dates[j], ".tif", sep="")
    if (length(in.hdf)>1) {
      mosaicHDF(in.hdf, out.hdf, MRTpath=MRT.PATH)
    } else {
      file.symlink(in.hdf, out.hdf)
    }
    reprojectHDF(out.hdf,
      out.tif,
      MRTpath=MRT.PATH,
      proj_type="UTM",
      proj_params="0 0 0 0 0 0 0 0 0 0 0 0 0 0 0",
      utm_zone=0,
      datum="WGS84",
      bands_subset=MODIS[["bands"]],
      pixel_size=MODIS[["pixel.size"]])
    ## unlink(out.hdf)
  }
}

Now I have a timeseries of Geotiffs in UTM projection I can work with.

Reading GPX tracks in R

I found another post on R-bloggers, which is very similar to my previous post. So I modified my script and it does not use the library plotKML anymore:

library(XML)
library(raster)

getTrack <- function(file, vars=c("time", "ele")) {
  pfile <- htmlTreeParse(file, error=function (...) {}, useInternalNodes=TRUE)

  trkpt <- xpathSApply(pfile, path = "//trkpt", xmlAttrs)
  create.df.str <- 'out.df <- data.frame(lon=as.numeric(trkpt["lon",]), lat=as.numeric(trkpt["lat",])'
  for (n in vars) {
    if (n=="ele") {
      ele <- as.numeric(xpathSApply(pfile, path = "//trkpt/ele", xmlValue))
      create.df.str <- paste(create.df.str, ", ele=ele", sep="")
    } else if (n=="time") {
      time <- xpathSApply(pfile, path = "//trkpt/time", xmlValue)
      time <- strptime(time, format = "%Y-%m-%dT%H:%M:%OS")
      create.df.str <- paste(create.df.str, ", time=time", sep="")
    } else {
      eval(parse(text=paste(n, ' <- xpathSApply(pfile, path = //trkpt/', n, ', xmlValue)', sep="")))
      create.df.str <- paste(create.df.str, ", ",n, "=", n, sep="")
    }
  }
  create.df.str <- paste(create.df.str, ")", sep="")
  eval(parse(text=create.df.str))
  return(out.df)
}

track <- getTrack("track.gpx")

track$delta.dist = 0
track$delta.dist[2:nrow(track)] = pointDistance(track[2:nrow(track), c("lon", "lat")], track[1:(nrow(track)-1), c("lon", "lat")], lonlat=TRUE)

track$delta.time = 0
track$delta.time[2:nrow(track)] = as.numeric(difftime(track$time[2:nrow(track)], track$time[1:(nrow(track)-1)], units="secs"))

track$speed = 3.6 * track$delta.dist / track$delta.time