Thursday, February 6, 2020

R Basics

Tips and tricks in R captured in the last 6 months of work at Mazama Science. Most of this is not rocket science, but some of it took research to figure out and I’d like to keep it in my mind.

Basic Data Structures

Variables

foo <- "Hello World!"

Note the strange assignment operator “<-”. The equal sign works as well, but by convention is only used with function parameters.

Vectors

Basic unit of data in R. Pretty much the same as a Python list or an array in Perl (without the semantic stupidity). Vectors are ordered and must contain the same type of data. Elements in a vector can be accessed by index location, but using a negative index means “everything but that element”. Oh yeah, and R is not zero-indexed. Trying to access the [0] element yields no joy.

> foo <- c(1,2,3,4)

> foo[3] == 3
# TRUE

> foo[-1] == c(2,3,4)
# TRUE

> foo[-4] == c(1,2,3)
# TRUE

Lists

R lists are like Python dictionaries, or hashes in Perl, but with the important distinction that they retain their ordering. Lists can also store multiple data types.

1. Create a list using the list() function:
> contact_info <- list("first_name" = "Roger", 
                     "middle_initiual" = "T", 
                     "last_name" = "Shrubber", 
                     "address" = "123 Jump St.",
                     "age" = as.integer(100))

2. Access 1st element by index position (NOTE that this returns a list):
> contact_info[1]

$first_name = 'Roger'

# 3. Access list element by name:
> contact_info["first_name"]

$first_name = 'Roger'

4. Access first 3 elements in the list:
contact_info[1:3]

$first_name
'Roger'
$middle_initiual
'T'
$last_name
'Shrubber'

5. Access just the values, without their names (keys)
unlist(contact_info[1:3], use.names = FALSE)

'Roger' 'T' 'Shrubber'

6. Access just the value of an element using the $<name> syntax
contact_info$address

'123 Jump St.'

Lists are very useful as lookup tables and as catch-all containers for random crap. Lists are slower than vectors though, so should be used with discretion.

Matrices

R implements a 2-dimensional array using the matrix() function. Higher dimension arrays are possible as well, but since I haven’t used them yet, I’ll stick with what I know.

1. Create a 3 x 3 matrix using a sequence from [1,9]:
> threesome <- matrix(seq(1,9),
                    nrow = 3,
                    ncol = 3,
                    byrow = TRUE)
# yields
1	2	3
4	5	6
7	8	9

2. Access row 1 of matrix:
> threesome[1,]

# yields
1 2 3

3. Access column 1 of matrix:
threesome[,1]

# yields
1 4 7

4. Access a "window" of values. [row_start:row_end, col_start:col_end]
threesome[2:3, 1:2]

# yields
5	6
8	9

# NOTE: Setting "byrow = FALSE" when you create the matrix yields this:
1	4	7
2	5	8
3	6	9

Dataframes

Dataframes are probably the reason you started using R in the first place. They are the type of data structure that you seem to encounter the most when working with packages.

first_names <- c("Fred", "George", "Ronald", "Harry")
last_names <- c("Weasely", "Weasely", "Weasely", "Potter")
addresses <- c("27 Foobar Way", "27 Foobar Way", "27 Foobar Way", "4 Privet Dr")
counties <- c("Devon", "Devon", "Devon", "Surrey")
ages <- as.integer(c(17, 17, 15, 15))

potterTrivia <- data.frame("firstName" = first_names,
                           "lastName" = last_names,
                           "address" = addresses,
                           "county" = counties,
                           "age" = ages, 
                           stringsAsFactors = FALSE)

# yields
	firstName	lastName	address	county	age
    <chr>	<chr>	<chr>	        <chr>	<int>
1	Fred	Weasely	27 Foobar Way	Devon	17
2	George	Weasely	27 Foobar Way	Devon	17
3	Ronald	Weasely	27 Foobar Way	Devon	15
4	Harry	Potter	4 Privet Dr	Surrey	15

Like everything else, R dataframes are indexed and data inside them can be accessed by their index position.

  • By column:
potterTrivia[1]

# yields
firstName
<chr>
Fred
George
Ronald
Harry
  • By row:
potterTrivia[1,]

# yields
firstName	lastName	address	    county	age
    <chr>	<chr>	<chr>	        <chr>	<int>
1	Fred	Weasely	27 Foobar Way	Devon	17
  • By range of rows:
potterTrivia[1:3,]

# yields 
firstName	lastName	address	    county	age
    <chr>	<chr>	<chr>	        <chr>	<int>
1	Fred	Weasely	27 Foobar Way	Devon	17
2	George	Weasely	27 Foobar Way	Devon	17
3	Ronald	Weasely	27 Foobar Way	Devon	15

  • By range of rows and columns
potterTrivia[1:3, 1:2]

# yields
  firstName	lastName
    <chr>	<chr>
1	Fred	Weasely
2	George	Weasely
3	Ronald	Weasel

Column names can be viewed from a dataframe by using the names() function and $column_name can be used to extract it from the dataframe.

names(potterTrivia)

# yields
'firstName' 'lastName' 'address' 'county' 'age'

potterTrivia$firstName

# yields
'Fred' 'George' 'Ronald' 'Harry'

In many ways, a dataframe acts and feels like a database table. So it should come as no surprise that Pandas, which was developed to provide R-like functionality in Python, is described by some as, “an in-memory nosql database that has sql-like constructs”. I think similar things could be said about the R dataframe.

Written with StackEdit.

Friday, January 24, 2020

Creating Rasters in R and Python

Having the ability to create rasters yourself is useful, especially when working with raw sensor data that hasn’t been conveniently pre-packaged.

Basic Raster Concepts

Rasters are 2 dimensional data structures in which data is contained in a regularly spaced grid. Each cell in the grid contains a single data value which can be retrieved by specifying either the exact grid address (row and column), or the region (group of grid cells) in which the value is stored. Usually the spacing of the grid cells is uniform in both the X and Y directions, but this is not always the case. The formal definition for this type of data structure is a 2-dimensional array, or matrix.

Rasters are referenced from an origin that is either at the top or bottom left corner of the grid. The row and column numbers of each grid cell (or pixel) increase in count from the origin, but the Y position of the data may either go up or down as the count increases. When a raster is georeferenced, the origin corresponds to a geographic location expressed as a point in whatever coordinate system the raster has been defined in.

Fig 1. Raster grid with origin at top-left

Raster grid with top-left origin

Various types of data can be stored in rasters but commonly they are used to store images, in which the cell values are representations of RGD color values (in the case of color images), or light reflectance or brightness (in the case of B&W images). Another common use for rasters is to store scientific data such as land elevation, air temperature, smoke density or water depth. These are examples of continuous data, where each pixel can represent a different and continuously variable value. Rasters can also be used to store discrete values, where a group of pixels represent a discrete class of data, such as a land-use category, where each pixel in a specific class will have the same value. Finally, rasters can also be used to aggregate data from various point locations, like from a sensor network, where each sensor location is represented by a pixel in the raster. This type of point location data is often converted into a continuous data set later by extrapolating pixel values between the point locations.

Fig. 2 Continuous data in raster as a hillshade

cont_color_hillshade

Elements Required to Create a Raster

Very little information is needed to create a basic raster. All that is really needed is some data arranged in a matrix. (Note that in Python the term “array” is used instead of “matrix”.)

R

library(raster)

simple_matrix <- matrix(seq(1,9),
                        nrow = 3,
                        ncol = 3,
                        byrow = TRUE)

simple_raster <- raster(simple_matrix)

Fig 3. Simple 3x3 raster in R

R_simple_3x3

Python

from PIL import Image

simple_array = np.asarray(range(1,10)).reshape(3,3)

im = Image.fromarray(np.uint8(simple_array))
pyplot.imshow(im)

Fig. 4 Simple 3x3 raster in Python

python_simple_3x3

In both examples above, a simple 3x3 matrix was created from a sequence of numbers that range from [1, 9]. That matrix was then passed into an image creation library and the resulting image was then rendered. The important concept to take away from this is that we can control the content of a raster simply by controlling the content of the numeric matrix that is used to create it. Here are a few examples.

Creating a Raster with a User-Defined Matrix

R

simple_grid <- c(1,0,1, 0,1,0, 1,0,1)

simple_matrix <- matrix(simple_grid,
                 nrow = 3,
                 ncol = 3,
                 byrow = TRUE)

simple_raster <- raster(simple_matrix)

# plot(simple_raster, legend = FALSE, axes = FALSE)

plot(simple_raster)

Fig 5. User defined 3 x 3 raster in R
defined_R_simple_3x3

Python

user_array = np.array([[1,0,1],
                       [0,1,0],
                       [1,0,1]])
user_im = Image.fromarray(np.uint8(user_array))
pyplot.imshow(user_im)

Fig. 6 User defined 3x3 raster in Python
defined_Python_simple_3x3

In these examples above, a user-defined 3 x 3 matrix is converted into a raster. Again, the important thing to remember is that we can control the raster content by modifying the numeric matrix that it’s created from. Below we can see how a single pixel value can be modified.

Modifying a Raster via Matrix Address

R

simple_matrix[2,2] <- 0
user_raster <- raster(simple_matrix)
plot(user_raster)

Fig. 7 Modifying a single pixel by matrix address in R
R_matrix_address

Python

user_array[1,1] = 0

user_im = Image.fromarray(np.uint8(user_array))
pyplot.imshow(user_im)

Fig. 8 Modifying a single pixel by matrix address in Python

python_matrix_address

Next

In my next post, I’ll cover how rasters are georeferenced and some of the ways in which having a spatial context allows us to use rasters to do interesting things.

Monday, January 20, 2020

Coding on a Chromebook

Chromebook Dev. Setup

Here's a blow-by-blow of setting up a Chromebook that can be used for coding in Python or R with Jupyter Notebooks.

Step 1. Get into dev mode

  • Hold ESC and Refresh at the same time, then hit Power.
  • Ctrl-D when screen with yellow "!" appears

Step 2. Update .bashrc file

While the aliases aren't needed, they're useful and once you get used to them, they're hard to live without. The exports are needed later on and you might as well add them now.

alias cls='clear'
alias gh='history | grep'
alias ll='ls -laF'
alias h='history | tail -n 10'

# Below are needed for R to work correctly
export EDITOR='vim'
export TMPDIR=/usr/local/tmp
export PAGER=/usr/local/bin/less

Step 3. Set developer passwd

This is recommended to help offset the vulnerability you create by running the machine in Developer mode.

$ sudo chromeos-setdevpasswd

Step 4. Install Chromebrew

Chromebrew is a package manager which allows you to install programs onto your Chromebook. It's similar to the Mac Homebrew package manager.

  • Download and install Chromebrew installer:

$ cd ~/Downloads # This is the only dir that is writable to you at the moment
$ curl https://raw.githubusercontent.com/skycocker/chromebrew/master/install.sh \
-o install.sh
$ bash ./install.sh

This will download and install a compiler toolchain which includes gcc and other useful items that are needed to compile source code on your machine.

  • At this point, you could have a look at what's available from Chromebrew here, https://github.com/skycocker/chromebrew/tree/master/packages.

NOTE: The list of installed packages is contained in /usr/local/etc/crew/device.json. It is often better to remove a package entry from here and then reinstall it, rather than 'crew remove <pkg", and="" can="" cause="" done="" gcc8="" level="" libiconv.<="" like="" low="" on="" p="" pkgs="" problems="" when="" which=""> </pkg",>

Step 5. Install some pre-requisites

NOTE: When installing from Chromebrew, I always use the "-s" flag to the install command. This builds the package locally from source, which I find to be much more reliable than installing from precompiled packages.

  • Install hdf5 library (needed for NetCDF)
 $ crew install -s hdf5
  • Reinstall libxml2 (used by ALOT of things. Best to reinstall it from source)
$ crew remove libxml2
$ crew install -s libxml2
  • Download and install NetCDF
    Go here, http://github.com/Unidata/netcdf-c/releases, and download the .tar.gz of the version you want. I'm using the latest (4.7.3) here. We need this before we build and install GDAL and R so that they have support for working with .nc files.
$ ll ~/Downloads/netcdf*
-rw-r--r--. 1 chronos chronos 18410112 Jan 18 11:06 netcdf-c-4.7.3.tar.gz

$ cd /usr/local/share
$ cp ~/Downloads/netcdf-c-4.7.3.tar.gz ./
$ tar xzvf netcdf-c-4.7.3.tar.gz 
$ cd netcdf-c-4.7.3
$ ./configure --prefix=/usr/local --libdir=/usr/local/lib64 --disable-dap
$ make check install
  • More prereqs that I'd prefer to install from source
$ crew install -s python27 openjpeg geos proj4 curl
# After removing from devices.json
$ crew install -s libiconv
$ crew install -s gettext
$ crew install -s cairo
  • Make sure numpy is installed before installing GDAL
$ pip install numpy
  • Install GDAL
$ crew install -s gdal

Step 7. Install R

$ crew install -s r

NOTE: This will add a ton of other packages as part of the dependency tree.

Step 8. Install Jupyter

This one is easy. And by all that is holy, do not infect your machine with Anaconda. It might work fine on a normal machine, but here it will create havoc.

$ pip install --upgrade pip
$ pip install jupyter

Step 9. Install IR Kernel

This you do from within R, but first we need to install the system libs for ZeroMQ, and they were incompatible with the latest version of gcc 8 that I had.

Had to download Version 4.3.2 from https://github.com/zeromq/libzmq/releases/download/v4.3.2/zeromq-4.3.2.tar.gz

$ crew install -s zeromq 
$ R
> install.packages("Cairo")
> install.packages('IRkernel') 
> IRkernel::installspec()
  • Set cairo as default device in IRKernel
$ vi /usr/local/lib64/R/library/IRkernel/kernelspec/kernel.json
{"argv": ["R", "--slave", "-e", "options(bitmapType='cairo') ; 
"IRkernel::main()", "--args", "{connection_file}"],
 "display_name":"R",
 "language":"R"
}

Step 10. Install Python dependencies

$ pip install matplotlib

Step 11. Debug

  • R notebook won't open because it lacks a PNG device
$ crew remove cairo
$ crew install -s cairo
  • Fix Python ascii decode error
$ vi /usr/local/lib/python2.7/site-packages/nbformat/sign.py

# Add to the beginning of the script
reload(sys)
sys.setdefaultencoding('utf8')

Tuesday, December 10, 2019

The Beauty of Signed vs. Unsigned Integers

The Beauty of Signed and Unsigned Integers

"The time has come", the Walrus said, "to talk of many things, but mostly about signed and unsigned integers."

Well, he didn't really say that and I feel a bit guilty about butchering one of my favorite poems, but that's what I'm going to talk about here.

I think this is a topic that most modern coders don't think about a whole lot, but I'm sure that people who work in compiled languages, or who have been around the block a few times, will chuckle appreciatively about this. It recently bit me in the tuckus and I thought it would be useful to talk about it. But first, why does it matter?

RTFM

The data product I've been working with recently is produced by NOAA and is extremely well documented. It's one of the products generated by the GOES satellites named AOD (Aerosol Optical Depth) and a key point the NOAA User Guide calls out is this.

"The AOD and its valid range are stored as packed 16-bit unsigned integers in the product files. This may be important when reading the data as some software may interpret them as signed integers. The packed values of AOD should have a valid range of [-0.05, 5]. "

Imagine my surprise then, when I plotted some of this data and got something that looked like this.

old_AOD

It's very pretty, and it shows the smoke being produced from the Kincade Fire in Sonoma CA nicely, but notice the range of the color bar. It shows a range of roughly [-2.5, 2.5], which is just not right. (BTW, the square brackets indicate a closed interval.)

Look at the Data, Dummy!

The NOAA User Guide states that the AOD data are, "stored as packed 16-bit unsigned integers". I'll talk about "packing" in a separate post, but for now, if we just look at the "16-bit unsigned integer" spec, that means I should have raw data values in the file that range from [0, 65536].

However, when I looked at the raw data I had extracted from the file, I found that I had Min and Max values of -31128 and 32750, respectively, and a whole lot of other numbers in between. Here's a little bit of the raw data, dumped out of the AOD file.

<snip>
    _, _, _, _, _, _, _, _, _, -30111, _, _, _, _, _, _, _, _, _, _, _, _, _, 
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, 12079, 12714, _, _, _, _, _, _, 
    4837, 3861, 3509, _, _, _, _, _, 3167, _, 3519, 4368, 2789, 3478, 2406, 
    3080, _, 2784, 2407, 3194, 2592, _, _, _, _, _, _, _, _, _, _, _, 3845, 
    _, _, _, 3712, _, _, _, _, 3015, 2937, 2818, _, 2981, _, _, _, _, _, _, 
    7934, 7676, _, _, _, _, _, _, _, _, _, _, _, _, _, 16670, 12271, _, -3, 
<snip>

(NOTE presence of -30111 value in 1st line of output and -3 in last line)

Signedness

The Wikipedia page describes signedness better than I can, but in a nutshell, if you have "unsigned" integers, there should be no negative numbers present. If there are, it means that the data is being stored as a "signed" type. This is totally fine, and in fact, certain data types, like NetCDF for example, silently convert unsigned integers to signed when data is loaded into them. It's up to the data vendor to use metadata to notify the user that the data is unsigned, and in fact, the GOES data has this metadata tag,

AOD:`_Unsigned` = 'true'.

Knowing that I shouldn't be seeing negative values in my data, but that I was, and that the data dictionary explicitly stated what I should be seeing, I decided to do some testing. Sure enough, if I loaded the following values into a NetCDF I had defined as being "16-bit unsigned,

63896, 50000, 0, 32128

I got the following out when I read the file,

-1640, -15536, 0, 32128

And that was the key to figuring out the problem.

The Solution

As I mentioned above, the maximum value I should have been able to store as a 16-bit integer was 65536. (In actuality, it's 65535, but you can read about that here Any value that is negative and coming originally from an unsigned data set means that the original value was larger than the signed data maximum. (I know, I know, what the hell does that mean?)

16-bit signed integers have a range of [−32768, 32767]. So a value that is negative means that it was greater than 32767 originally.

You can see that in my example above, where 63896 signed got stored as -1640 unsigned. Basically the negative side of 0 is used to store the offset from the maximum unsigned value.

max_uint - offset = original value

65536 - 1640 = 63896

In any case, the simple solution to convert back is to add 65536 to every negative value in the raw data. Astute readers will notice that my original Min and Max values of -31128 and 32750 then convert to 34408 and 32750, which seems weird at first. However, recall that as the negative signed numbers get closer to zero, the closer they get to the maximum unsigned value. So in the little data snippet I dumped out of the file, there was a -3. If we convert that to unsigned, we have 65533, which is actually very close to the maximum unsigned value possible.

Moral of the Story

  • Look at your raw data very closely before assuming the code is incorrect
  • Better yet, write unit tests against test data you create, not operational data
  • Understand the data you're working with and read the documentation about it

We assumed initially that there was a problem in the code we'd written, or that we had somehow received scaled data from NOAA. Going back to look at the raw data, before ANY scaling was applied by the libraries was the key to discovering what was happening. This is why it's good to create test data that you can write unit-tests against, to validate that your code is performing "as expected". More importantly though, really, really understand the data you're working with. This was an amazing learning experience for me, one that I really enjoyed, but it could have been incredibly frustrating.

Next posts will be on working with the NetCDF data format and on the concept of data "packing".

Displaying Rasters in Jupyter Notebooks

Rasters in Jupyter

It's the easy stuff that trips you up the most, right? Here I was, wanting to talk about creating rasters in R and Python when I realized that I didn't have a great handle on how to display them. Turns out that it is very simple to do. (GeoTIFF generated from code taken from Jared's excellent gdal/ogr cookbook.)

In Python

Use GDAL and Matplotlib

from matplotlib import pyplot
%matplotlib inline
from osgeo import gdal

img = gdal.Open('test.tif').ReadAsArray()
im = pyplot.imshow(img)

python<em>jupyter</em>geotif

In R

Use the raster library's built-in plot() function

library(raster)
library(repr)

options(repr.plot.width=4, repr.plot.height=3)
img_file <- raster("test.tif")
plot(img_file)

NOTE: repr is used solely to control the output size of the plot. Without it, the image is rather large.

R<em>jupyter</em>geotif

Rebirth of the Blog

Rebirth of the Rants

Rebirth of the Rants

Well, maybe not so much the rants, but since going back to work, I've felt the need to capture a few of the things I've had to re-learn for posterity. Maybe one day I'll organize them and use them for tutorials, or something, but in the meantime, at least I'll have captured them somewhere that I can search and easily share them.

I had been using Github personal pages to do something similar, but I could never get into the swing of updating the pages there. The Pure.css templates that I used on that site looked great, but there was overhead in maintaining them AND, as my friend Mash pointed out to me years ago, it's stupid not to leverage some of the built-in tools that Google provides on Blogger.

This time around, I'm going to draft all my entries in Simplenote as Markdown and then cut-and-paste the HTML they generate into Blogger. I use SImplenote already and while it's not perfect, I like it. I did a little write-up on it a year ago and it's probably time to do another that captures what I've learned about using it since then, both good and bad.

There's a great deal of old and krufty stuff in the blog right now, but I'll whittle it down as I move forward. For now, I'll focus on getting new material up that I think is relevant.

So without further ado, let's get to it!

Friday, October 16, 2009

Using GRASS to fix topology

Started trying to use GRASS to solve some of the more difficult problems that I can't fix in PostGIS. Here's one of them. I've imported a shapefile of country borders and am trying to fix the topology. The screenshot shows a country border with a bunch of "x" marks that I *think* represent centroids. I'd like to take them out, and just have one, long, continuous line at the country border. Not sure how to do that yet.