in which ESRI does not die so much as expire from boredom

I've been loving using Quantum GIS lately.  Wow, look out Manifold.

Listing of Google imagery updates

This blog has a KML file showing areas of imagery update on specific dates:

http://www.gearthblog.com/blog/archives/2012/06/large_new_update_to_google_earth_im.html?utm_source=feedburner&utm_medium=feed&utm_campaign=Feed%3A+GoogleEarthBlog+%28Google+Earth+Blog%29

\dontask

## \dontask{


suppressMessages(library(mgcv))

invisible(capture.output(library(deldir)))

invisible(capture.output(library(spatstat)))

library(foreign)

library(lattice)

library(sp)

library(trip)

suppressMessages(library(rgdal))

suppressMessages(library(spatstat))

invisible(capture.output(library(maptools)))

suppressMessages(library(geosphere))

## }

read LAS data with R

Finding software to read LAS files can be difficult, support for 3D data is still hard to come by.

R's support for the "raw" data type provides a very neat way of reading data from binary files, such as LAS files. These files have a header containing metadata that is used to specify reading from the bulk data in the file. Essentially, you can read large sets of records at once into a raw matrix with readBin and then use the indexing tools to extract the relevant bytes for each element of the record.

This means that the bulk read is fast, and processing is also vectorized as the raw matrix can then be read in index pieces by readBin. I've done this in a simplistic way in the following functions.

• "publicHeaderDescription" - function returns a list generated to hold information about the header 
• "readLAS" - read file function

It works for some LAS files that I have and returns a matrix of record values. Certainly it works fast and efficiently enough for smallish files. I really don't have experience with the likely range of sizes of these files.  I've included "skip" and "nrows" arguments to allow for chunked/partial reading.


The remaining work I see involves:

• wrap the "chunked" read so large files can be processed in parts
•  return only specified attributes from each record (can be x, y, z, gpstime, intensity)
• extract all components of the records, and generalize for LAS 1.1, 1.2, ... (there is minor handling of gpstime presence or absence, but nothing else)
•  re-organize the header read and arrangement (what I've done was easy enough for me to understand, but it's a bit of a mess)
• (no doubt many more general improvements I haven't thought of)

To use:

source("http://staff.acecrc.org.au/~mdsumner/las/readLAS.R")
## sample LAS file from here: http://liblas.org/samples/
f <- "Lincoln.las"

file.info(f)$size/1e6
## [1] 185.5660 Mb
system.time({
  lasdata <- readLAS(f)
})
# user system elapsed

# 4.04 0.62 6.93
dim(lasdata)
# [1] 9278073 4

colnames(lasdata)
#[1] "x" "y" "z" "intensity"

Optionally the list of header components (some extracted as actual data, some left in their "raw" form) can be returned instead of data using "returnHeaderOnly=TRUE".

Read data directly from Manifold tables using R

There are (at least) two ways to do this:

1. Via the clipboard

Select the records you want from a table in Manifold, Copy

From R read the data from the clipboard using read.delim the default separator is a tab, so it works already - but we can specify that and other options via arguments.

>d <- read.delim("clipboard")

"d" is now a data frame in R with all the data from the Manifold table. This works for virtual tables of images and surfaces, but it doesn't work for the binary geometry columns (you will get a text summary of those).

2. Using RODBC.

This option requires R with the RODBC package installed, get these from a CRAN mirror via http://www.r-project.org, and a function defined for connecting to a save Manifold project file (see below).

>library(RODBC)

Open a connection to a .map file (I copied this from the function odbcConnectAccess)

>ch <- odbcConnectManifold("C:/temp/world.map")

Submit a query to the connection, and assign the results to variable "d"

>d <- sqlQuery(ch, "SELECT [Longitude (I)] AS lon, [Latitude (I)] AS lat, [Capital] AS name FROM [Countries] WHERE [Area (I)] > 150;")

(We convert the field names to sensible character strings for ease of use later).

"d" is now a data frame with the (centroid) location coordinates and city names of all objects from the "Countries" drawing in "world.map"

My original posting on Georeference, including the function required is here:
http://69.17.46.171/Site/Thread.aspx?id=29419&ti=632979299200630000

My first blog

Well this is it, obviously I should have something to say on my first blogoff so check out the Manifold toolbar:

http://www.manifold.net/toolbar

Manifold is a GIS program I use, it is colourfully loved and hated. For the love join us at

http://www.georeference.org

For some hate, well - who needs it?