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.

Wednesday, October 14, 2009

My take on commercial vs. open source GIS

Was alerted to a posting on the LinkedIn OGC group that seriously chapped my hide. A person had asked for advice regarding whether they should consider Open Source as a viable alternative to commercial software. You can read the original post at this link.

Here's my response.

"GIS" is now a very broad topic, and means different things to different people. You didn't specify what type of work you needed to do, so I'll try to address a few major topics that I encounter routinely:
- data conversion and creation,
- spatial analysis,
- server based data storage and access,
- desktop visualization,
- paper map production,
- web-based static maps,
- web-based dynamic maps,
and finally a weird one,
- reprojection to and from local projections from around the world.

- Data conversion and creation:
Hands down I prefer OS for this. Although I admit that I have encountered difficulties with CAD data that was created in bleeding edge commercial software. Usually this is because the OS development team has had to reverse engineer a new data format. I feel that there is far greater ability to work with data programatically using OS tools, and in my experience these tools suffer fewer failures when working with very large data sets. Having said that, if all I ever did was convert data, all day, every day, I probably would buy FME by Safe Software. They support more formats than I'm even aware of, and make it (relatively) easy to setup automated conversion pipelines. I don't need that degree of interoperability, and I use GDAL for practically every task that falls in this category.

- Spatial Analysis:
I feel that OS can do most spatial analysis quite capably. I'm a huge fan of doing vector based spatial analysis using SQL in PostGIS. I do 90% of my analysis this way. For more complicated topological analysis, as well as raster analysis, GRASS works extremely well. Having said that, my hands down favorite tool for analysis remains ArcMap. I am still much, much, faster in ArcMap, and I find that it is much easier to see the results of a particular operation with it. However, I rarely use it, as I can't run it on my laptop, and we can't afford another license anyhow.

- Server-based data storage and access:
Unless you have already heavily invested in commercial systems, or require database storage of raster data, OS (PostGIS specifically) wins in my book. If your entire workflow is already based on using many ArcGIS-equipped workstations which all connect to an Oracle Enterprise DB, I would be hesitant about deciding to convert to OS. Not because OS can't do it (except for database raster storage), but because there will be so much work in converting both the workflow, and the storage system itself. I understand from others that this is getting better - ESRI now supports PostGIS via SDE for example - but it's not without its headaches. Database storage of rasters is essentially non-existent in OS. There are all sorts of "solutions" talked about, but none of them are great, and frankly it's a glaring hole in OS.

- Desktop visualization:
Depends. If you want to quickly see what a data set looks like, or style several layers together, OS is great. It meets all of my needs on a daily basis. Then again, I think the grayscale display of 32-bit data in OpenEV is fantastic. There definitely is a lack of refinement in this area, but I'm hesitant to say OS isn't as good, mostly because I wonder just how sophisticated it needs to be?

- Paper Map Productionn:
For the frequent production of "casual maps", I think commercial software is still much better. I say "casual" because if you're making professional quality maps, I doubt you're going to just use a commercial GIS package. More than likely, you are also going to invest in a commercial graphics application as well. For "casual" paper maps, it's alot easier to use ArcGIS. Most of the time the results look pretty good, although the PDF engine is terrible and anything can happen. Personally, I use MapServer and GIMP to create paper maps with OS tools. I think they look pretty good, but they take more work than they should.

- Static web maps:
You're absolutely silly if you use anything but one of the OS tools for this. They can create lovely images - better in some cases than commercial products - and the only expense is the setup and configuration time.

- Dynamic and "slippy" web maps:
And now we get to what has everyone up in arms these days - everyone wants a mashup. Everyone is getting into this market, and everyone touts the capabilities of their solution at the expense of all others. I think that again it depends on what your goals are. If you want to create a web-based tool for spatial analysis, then I think the commercial offerings provide more "all in one" capability. If on the other hand, your goal is to display single or multiple styled layers, to allow for feature Identification and Selection (ESRI definition intended) - either by map click or attribute search, then there is no reason to use commercial software. In fact, I actually think that for these situations OS can both look and perform better. There are even OS base maps now that rival any of the commercial ones. I think there is a HUGE (note use of CAPS) market for, "3 vector layers and one image" on a web map, and I'm ecstatic to hear that a major vendor of commercial software doesn't think so.

- Reprojection support:
This frankly is an edge category, but's it's caused me problems numerous times. I suspect most people won't care though. OS software supports hundreds of common (and some not so common) projections via the PROJ library. Not only that, but it is easy to create your own custom projection as well, to display a specific area exactly how'you'd like. However, it does not deal with the reprojection of certain ones very well, especially those which use old, grid-based datums. Commercial software definitely still does a far better job at doing this.

Finally, a couple general observations about OS vs. commercial software. One of the major differences between commercial and OS applications is that in general OS applications do not try to address every conceivable situation you might possibly encounter in your work. Instead, they tend to focus on doing a few things extremely well, then rely on the ability to access other tools when the need arises. For example, Quantum GIS does a good job of displaying raster and vector layers, but relies on the ability to access GRASS in order to do more complex spatial analysis. If you are looking for a single solution for everything, commercial software with its multitude of extensions might be for you.

Support is also a good topic worth mentioning. In my experience, the people who like and insist on support contracts (aside from the vendors, of course), are the ones who never actually use the software. It's the ones buying the software that want it, and they view it as insurance that if something goes wrong with their very expensive purchase, they will get support. And usually there is a Service Level Agreement (SLA) that specifies they will get support in a timely manner, the amount of time depending on whether they are "Gold" or "Platinum" contract holders. But an SLA doesn't guarantee the quality of support - it can't. It merely ensures that the customer's problem will be dealt with in whatever support process the vendor has in place. Some companies are great, some are not, and it has nothing to do with how much money you paid for your support contract. In general the support given to OS applications by the developers and users of the product is as good, OR BETTER (note use of CAPS for emphasis) than any I have received via contract. Support extends beyond just staffing a help desk and charging a monthly maintenance/support fee though. It includes documentation of the product, and rapid bug fixes for critical problems.

Documentation for OS applications can sometimes be tough to find, and the quality is variable. There definitely is documentation out there, but sometimes it takes more effort than it should to find, and then it's not always easily digestible. You won't find the Help system that's in ArcGIS, that's for sure. But you will find many well-written Wiki's that cover most topics, and very responsive user lists to answer specific questions. Bug fixes and feature enhancements are far more rapid in the OS world, and tend to be in response to problems and requests reported by users. I think it's Documentation can sometimes be tough to find, this is true. There definitely is good documentation out there, but it is sometimes hard to find, and not easily digestible. You won't find the Help system that's in ArcGIS, for sure. But you will find many well-written Wiki's that cover most topics, and very responsive user lists to answer specific questions.a function of the development process that enables this. The bug tracking systems are transparent, and at any time you can go take a look at what's being done to resolve your problem, or to see what problems exist in a specific version. And lets not forget that "money talks" in the OS world as well. You absolutely need to be able to connect to that new/old/weird/unsupported database with your OS app, and you need it done NOW, and are willing to pay $5000 for the priviledge? There's a really good chance someone in the development community will do it for you.

So, to wrap this thing up, Open Source GIS software is just like anything else in life, it has both good and bad points. The same is true of commercial software as well. I don't think you can generalize across the entire category of uses with simple blanket statements. Identify what specific uses you need to address, then ask the question again. You'll probably get different answers for each one.

Best of luck.