Lotte GIS!

Python/Tkinter/ArcPy programming notes

This is what we learned (but not necessarily tested on) for using Python with ArcGIS software: Google Doc link

Map algebra, cont’d

I want to make an addendum to the map algebra chat from a previous post…

Apparently the syntax I used is a remnant of previous versions of ArcGIS. Our instructor had taught us the more clumsy one because it’s the way forward. The v10 documentation for the Con tool has a little look to the past when it shows the alternative (old) syntax in a sample expression.

I hope ESRI keeps the old syntax in some way. Maybe workflow bugs will be worked out in later version updates, but currently, I find workarounds are necessary to get some expressions working. For example, using the Con tool directly or called from the Raster Calculator to set false-values to NoData in the new syntax will output an error.

Con("%raster%", 1, NoData, "%raster% > 100")

No problem with this expression:

Con("%raster%" > 100, 1)

Not only do you not need to set the false-value because it’s already NoData by default, but you don’t have to bypass it at all when writing the where-clause. Otherwise, the false-value parameter has to be bypassed somehow, either explicitly assigning NoData or with "". This is what seems to produce an error. Is it a known bug?

Map of the Slesse Creek Watershed in Chilliwack, BC

This is another map, created recently for a class assignment on watershed spatial analysis. This time, I figured I commit fully to producing a proper cartographic map by use of a proper illustration tool. In comes Adobe Illustrator!

Here are some things I did to prepare for exporting to AI:

• set colour ramp for the hillshade basemap,
• set classified colours of the watershed raster,
• edited and extracted water and road features from given data to highlight,
• in ArcMap Layout View, set up the map scale and and positioning.

It felt great using Illustrator to manipulate type, colour, and other graphical elements. No more redundant clicking just to change a font style! I could add and manipulate labels to the roads and streams. Perhaps more importantly, the legend could be customized to make it more intuitive.

GIS as an Information Visualization

There was a presentation assignment we had to do last semester where we each select a GIS topic to present to the rest of the class. I decided to share this since I still receive positive feedback on it from classmates. Click here for a copy of my presentation as well as my talking points.

I think every map maker should have a background in information visualization. Cartographers have some of the thought processes and discipline as information visualizers and essentially, both try to communicate by turning data into information; however, the latter have been advancing map design in leaps and bounds, by the likes of New York Times, with free online data visualization tools, etc. Sometimes it’s fine to put out a map from a template, stick it in a report, and hope the work brings the message across, but other times, there should be emphasis on making the case by rethinking how best to present spatial information. I think applying some principles from information visualization will help do that.

Some cartographic models I made + a little chat about map algebra

Here are some cartographic models I made for my classes. I didn’t like using ModelBuilder from ArcGIS (mainly because I couldn’t figure out how to change the symbol shape for tools), and created these models from Google Docs instead. The UI, aesthetics, and ease of use for Drawing documents actually make them a joy to create.

The first model deals with remote sensed imagery.

In the second model, all the tools use ArcGIS’s implementation of map algebra for spatial analysis, entered via the Raster Calculator tool. Learning its unique syntax and discovering quirks were like taming a beast, oy vey… Such that I came out of it learning a better way of writing statements (at least, to me). Instead of writing:

Con("%input%", 20, "%output%", "%input%" == 1)

one could write it simply like this:

Con("%input%" == 1, 20, "%output%")

The latter format is more understandable to me, where the conditional is stated first, then you determine what happens if it’s true. Also, less typing!

SQLPlus: Using substitute column as variable in conditional expression

I have a query that I thought I should optimize. It uses an arithmetic operation to calculate a difference in values and a conditional statement to filter out the calculation. The query is this:

SELECT citycode, (pop2001-pop1996)
    FROM population 
    WHERE (pop2001-pop1996) > 10000;

I wondered if I could instead simplify the statement and output to something along the line of this:

SELECT citycode, (pop2001-pop1996) AS popincr
    FROM population 
    WHERE popincr > 10000;

so that the arithmetic operation (which can be more complex as queries can be) doesn’t have to be duplicated in the conditional statement. I looked for a possibility but encountered this explanation why that isn’t possible from a StackOverflow user comment:

“Because the column names are resolved at compile time not at run time for the SQL statement.”

Ah… Well, that’s out the window.

On the other hand, in my search for an answer, I came across this very understandable article on how to use variables in SQLPlus: Spice up your SQL Scripts with Variables. It’s in my reading list.

Fundamentals of Multimedia

Is there a book I wish I hadn’t sold? Is there a book I miss right now? Yes…

As I was going through the past semester of BCIT’s accelerated GIS program, I kept trying to find a textbook from my BSc. studies. It was for a course that taught image data storage techniques and how to implement them. It was where theories of discrete signal processing were finally applied and I loved that class. But after getting my degree, I thought I sell my years’ worth of books. It managed to be one of the few to sell because it was still relevant by then.

It is apparently still relevant now as I went through the remote sensing course. As I ended up being a web developer, I never thought I would use signal processing again. Crikey, it all came back! Low-pass filters? High-pass filters? Look-up tables?

I remembered which book it was today. I want to get back “Fundamentals of Multimedia” by Ze-Nian Li and Mark S. Drew.

Numeric variable size types: SILFD

This is a reminder to self. I came up with this acronym to memorize the size order of numeric variables supported by ArcGIS. int is included for pronunciability of acronyms. It is actually a short int or long int.

SHORT: 2 bytes
INT: 4 bytes
LONG: 4 bytes
FLOAT: 4 bytes, fractional
DOUBLE: 8 bytes, fractional

Ref: http://help.arcgis.com/en/arcgisdesktop/10.0/help/003n/003n0000001m000000.htm

Problem: Attributes table doesn’t seem to show all the features as seen in the data frame

Let’s say a polygon was created while correcting topological errors. Double-checking if the polygon really was created, you opened the attribute table window - or had it opened already and refreshed the cache - and see that the number of objects is still the same. You’re sure there was a new polygon created because it’s selectable, you’re on the right layer, the table is set to show all records, etc.

Solution:

• Select one of the fields and apply a sort and maybe revert the sort as well. Doing this seems to get ArcMap to take a hard look at the features that feature class/attribute table has in order to update the table’s view. The new polygon should now appear in the table.

Possible reasons:

• You needed to save the edits in order to view changes to the attributes table.
• ArcMap buggered up the implementation of the “refresh cache” function.

wget a page

wget -v -e robots=off --wait 1 --ask-password --user=[username] -r -p -k -A [extensions list] [url]

• ignores robots.txt
• prompts password with HTTP Basic Authentication
• recursively crawl through accepted linked resources and saves onto local dir

FOUND: Photo base map of Vancouver.

Around a month ago, I thought I had to give up finding a suitable base map in which to edit my SHP files with. The problems were many-fold:

• free access to data
• data of raster layer of aerial, photographic Vancouver
• data that QGIS can import
• a GIS application that can edit my SHP file

Yesterday, I experimented on Google Map’s My Maps service to edit my SHP file. While I was able to edit it, still ended up not wanting to use Google Map or Google Earth based on 3 reasons:

• Google Earth is not free - to edit existing vector data, one must go Pro
• Editing in Google Maps is web-based and thus, cumbersome
• Saved work from Google Maps’ My Maps is not downloadable

In order to save to my PC, I clicked on the link to import into Google Earth. Shapes I had been editing become individual layers under the import folder. I could export that folder as KML or KMZ, but QGIS couldn’t import them. I then created a folder inside the import folder… and moved /every single shape layer/ into that folder, and exported that folder to KML. This, QGIS could import, but obviously all this work is far from ideal.

I was on the hunt again today and dug through GeoBC’s services more thoroughly (I also know more about GIS data than I did a month ago) and explored their WMS servers. Among the listed were imagery maps - hallelujah! Thank you, GeoBC!

Scrobbling radio tunes

This is amazing. A fellow last.fm scrobbler developed a little OSX app to scrobble radio tunes from iTunes: iTunes(Radio)Scrobbler.

The main radio I listen to is FM4 and it prefixes “[artist name] - [song title]” format with the name of their program, eg. “Morningshow: “. Obviously this skews the info to my last.fm library, so I made the following fix for my own use:

In iTunesRadioScrobbler.app/Contents/Resources/Perl/scrobbleTrack.pl:ln189

  if ($streaminfo ne $EMPTY) {
    # split streaminfo
    ($artist, $track, $leftover) = split /$inputSeparator/, $streaminfo, 3;
  }

to:

  if ($streaminfo ne $EMPTY) {
    # split streaminfo
    my $artistFull = $EMPTY;
    my $tmp = $EMPTY;
    ($artistFull, $track, $leftover) = split /$inputSeparator/, $streaminfo, 3;
    ($tmp, $artist) = split /\: /, $artistFull, 2;
  }

OSM PBF to PostgreSQL to QGIS

I have a SHP file of a city-level detail which I would like to edit. In order to do that, I’m juggling between uDig and QGIS to compare usability. I’m focusing on QGIS here, because it seems a lot of folks give it props. Here, from beginning to end, is what I did:

1. I have had no success finding satellite imagery of my city - Vancouver - through searching their list of WMS servers (the task does make me wish I was targeting Germany though).
2. I decided to look for a repository OSM/Mapnik maps… no luck there either.
3. I looked up “OSM + file” and found that there are snapshot dumps of OSM maps for download. Yes! However, the closest thing I find to Vancouver is at country-level. Alright, I’ll bite. It’s a 750meg PBF file, oh my…
4. I looked for a way to import that to QGIS and realized that while it supports OSM maps, it supports OSM maps in XML format. As of this date, OSM has already deprecated XML and switched to PBF. Great.
5. I browsed around looking for “QGIS + PBF” and came across osm2pgsql. This little tool takes in OSM PBF and converts/dumps the data to a PostgreSQL database, which QGIS can import!
6. I’m only versed in MySQL so I had to install PostgreSQL server and PostGIS. I forgot to mention that I’m on OSX so getting the database running took a couple of hours longer than usual.
7. Now finally I could get back to step 4 and did something with the PBF. This was my input:

   osm2pgsql -d canada_pbfimport -U postgres -H localhost canada.osm.pbf

   And the output:

   osm2pgsql SVN version 0.70.5
   
   Using projection SRS 900913 (Spherical Mercator)
   Setting up table: planet_osm_point
   NOTICE:  table "planet_osm_point" does not exist, skipping
   NOTICE:  table "planet_osm_point_tmp" does not exist, skipping
   Setting up table: planet_osm_line
   NOTICE:  table "planet_osm_line" does not exist, skipping
   NOTICE:  table "planet_osm_line_tmp" does not exist, skipping
   Setting up table: planet_osm_polygon
   NOTICE:  table "planet_osm_polygon" does not exist, skipping
   NOTICE:  table "planet_osm_polygon_tmp" does not exist, skipping
   Setting up table: planet_osm_roads
   NOTICE:  table "planet_osm_roads" does not exist, skipping
   NOTICE:  table "planet_osm_roads_tmp" does not exist, skipping
   Mid: Ram, scale=100
   
   Reading in file: canada.osm.pbf
   Processing: Node(71446k) Way(4608k) Relation(52k)  parse time: 1269s
   
   Node stats: total(71446261), max(1254505941)
   Way stats: total(4608725), max(109616802)
   Relation stats: total(52638), max(1553135)
   
   Writing way(4608k)
   
   Writing rel(52k)
   Committing transaction for planet_osm_roads
   Committing transaction for planet_osm_point
   Sorting data and creating indexes for planet_osm_roads
   Sorting data and creating indexes for planet_osm_point
   Committing transaction for planet_osm_line
   Sorting data and creating indexes for planet_osm_line
   Committing transaction for planet_osm_polygon
   Sorting data and creating indexes for planet_osm_polygon
   Completed planet_osm_roads
   Completed planet_osm_point
   Completed planet_osm_line
   Completed planet_osm_polygon
   
   

8. Now finally back in QGIS, upon quesstimating which table I should be importing, I picked one called public.planet_osm_polygon. After thinking a bit, it gave me this error:

   The table has no column suitable for use as a key.
   
   Quantum GIS requires that the table either has a column of type int4 with a unique constraint on it (which includes the primary key), has a PostgreSQL oid column or has a ctid column with a 16bit block number.

9. Argh!!!

Functional programming in Perl

I’ve been reading CGI::Ajax documentation and came across this example:

use strict;
use CGI::Ajax;
use CGI;

sub evenodd_func {
    my $input = shift;
    ...
}

where the function could also be written as:

my $evenodd_func = sub {
    ...
};

Awesome.

(Just want to put a disclaimer that I’m only working with Perl in order to fork an app I use.)