Tag Archives: spatial

get more done

Posted on by

currently working on multiple contracts has prompted me to reevaluate how i work. by that not the time spent doing technical work (thats important) but the way to manage time, define high level tasks & accompanying sub tasks that must be done to achieve completion.
time is the currency of a contractor. more of it = either more money and/or improved quality of life.

just a few of my current practices & ideas.
-dont check email first thing in morning. use that time for the biggest/most avoidable task. remember checking email puts you at someone elses needs.

-setup a spreadsheet to manage each contract. i have several for both work & other that have a few simple headings: time critical, todays priority,followup, nice to do & review.

- no distractions. only use 1 browser tab & turn phone upside down. i attempt to work in burst of 60 min & then a 5 min reward break. put away the social media & only check email a few times a day.

- use dropbox across all platforms. this way you will always have the critical files & pdfs etc to review when extra time permits. its easy & easy again.

- ask why are you doing this. ask repeatedly to ensure you make the best use of your time.

- plan & think about what you are working on. turn the monitor off for a bit & whiteboard etc what the core goal is & what actually needs to be done. this saves so much time.

- automate where you can.

these are a few of my recent work implementations & i will post others & resources ive used in the near future.
remember, its all about you figuring out what works for you & a system/practice you will follow for the long term.

cheers,

richard

Creating Wind Rose Diagram with R — Data Provided by QGIS

Posted on by

R is much more than just statistics. It offers very powerful & flexible graphics to be produced.
For a project I needed a way to display graphically the different azimuths that polygon objects contained. A wind rose diagram is used extensively in meteorology to show direction & magnitude of wind, you can think of it as a circular histogram with 0-360 limits. In this case it would be also provide a visual exploratory method for drumlin fields.

Note:

  • the polygon objects  in question must have an azimuth field. Using QGIS to do the initial spatial work, you can have R do the same.
  • in examples “std_bear” is reference. i have standardized bearing to 0-180. this means degree > 180, then subtract 180.

Basic workflow:

– QGIS

  1. Ftools to create polygon centroids
  2. mmqgis — transfer — export attributes to csv. ensure that you select the attribute with azimuth prior to export.

– R

  1. download http://www.exposurescience.org/pub/software/heR.Misc.tar.gz
  2. install package – change directory where file saved. >> install.packages(“C:\\1_WorkStation_Folders\\Downloads\\heR.Misc.tar.gz”,repos=NULL,type=”source”)
  3. library(heR.Misc)
  4. import csv of centroids. >> dataIn<-read.csv(file=”c://tmp//1_analysis_dh_dl/fullarea_bearing.csv”,header=T,sep=”,”)
  5. create vector containing only the azimuth data. >>dir_stdBear<-inData$std_bear
  6. plot wind rose with 9 bins (review). >> rose2(dir_stdBear,nplumes=9)

 

Output:

 

This is default output & there is many options to improve appearance.

What is interesting is that there are 5 field groupings & the predominate azimuth is easterly. From the drumlin data, I know there are approximately 5 large fields on the island. The percentages indicate amount of data.

 

Cheers,

Richard

 

 

 

pick your flavor of gvsig: both are good

Posted on by

another high quality spatial desktop software is gvSig. there are two versions, gvsig & a forked gvsig ce (community edition). both are really great but i have a preference for the ce edition (more tools). through the ce edition is a bit stale with older release, hoping this will change. the topology functions are really outstanding as is the sextante toolbox (victor olay is a very talented dev!). there are many functions that are not found else where ie. mean direction of lines, minimun enclosing shapes (the rectangle is min area based as opposed to saga’s min rect perimeter) & the geometric properties of lines that will calc bearing from north:)

from website:
gvSIG CE is a community driven GIS project fork of gvSIG that will be bundled with SEXTANTE, GRASS GIS and SAGA. This project is not supported by the gvSIG Association. gvSIG CE is not an official project of gvSIG.

you’ll need to increase the memory of the default settings to handle large projects. java based program, so set the -Xmx higher (look in the bin dir for the *.ini file).

it’s context based workflow, meaning if the attribute table is highlighted only those tools will be shown & likewise for the active layer.

note: when defining an initial layer, you must set the projection. select layer & click properties in menu tab to right.

in terms of capabilities, ive set -Xms to 1500M & processed geometries of lines of ~ 6 million line segments. pretty good!

R & PostGIS Together – Brief Code Snippet to get you going

Posted on by

The more I dive into using R for spatial & plotting the more it strikes me that it is a wonderful & powerful software toolset. I’ve been using Postgis (coupled with Postgresql) to hold both my vector & raster data. It is much easier to have a single repository for all your spatial data & not have to search amongst folders for that one GeoTiff or the dreaded directory full of shapefiles.

So why not combine R with Postgis? Below is a quick & complete code snippet that you need to get going. There are several ways to connect from R to Postgres, I have chosen to use the RPostgreSQL package. The code below assumes you have Postgis & R installed.

# load lib
library(RPostgreSQL)
# create db driver
dDriver<-dbDriver(“PostgreSQL”)
# open connection to db
conn<-dbConnect(dDriver,user=”enter here”,password=”enter here”,dbname=”enter here”)

# “dbSendQuery” allows SQL statement to be passed, in this case I am dropping tmp.dump table.
rs<-dbSendQuery(conn,”drop table tmp.dump;”)
# clearing returned object
dbClearResult(rs)

# create table tmp.dump
rs<-dbSendQuery(conn,”create table tmp.dump(ogc_fid integer,azimuth real,length real,width real, azimuthFrom text,azimuthTo text, centroidRect text);”)
dbClearResult(rs)

# fetching 1 column from db table
rs <- dbSendQuery(conn,”select ogc_fid from tmp.reduceftspace;”)
# “fetch” retrieves the returned db object from query & creates a dataframe
dbRows <-fetch(rs)
dbClearResult(rs)

# spatial query that converts polygon into pts feature
# query each geom using ogc_fid
ogc_fidin <- 1811
dbtable <- “tmp.reduceftspace”
sqlIndividualGeom <-paste(“select st_x(geom) as x, st_y(geom) as y from (select (st_dumppoints(the_geom)).* as g from”,dbtable,”where ogc_fid =”,ogc_fidin,”) as g;”,sep=” “)
rs <- dbSendQuery(conn,sqlIndividualGeom)
# fetch query results into dataframe
df_pts <-fetch(rs)
# clear results
dbClearResult(rs)

# disconnect from db
dbDisconnect(conn)

They should get you going to retrieve both tabular & spatial data from a Postgis/Postgresql database. Have fun!
Cheers,

Richard

 

 

Quick thoughts on transforming GeoTiff tags into parseable structure -> xml

Posted on by

Not a complete walk through here but some quick thoughts to a post for assist on twitter. What they wanted was to turn geotiff tags into csv or xml.

If you have GDAL installed & who doesn’t these days a quick way to turn the gdalinfo info (tags) into xml is to convert raster to Virtual Raster (VRT). This is an xml representation of raster.

>> gdal_translate -of VRT [destination name with .vrt extension]

This is what you end up with:

All that’s left to do is parse :)

<VRTDataset rasterXSize=”5529″ rasterYSize=”4607″>
<SRS>PROJCS[&quot;Lambert Conformal Conic&quot;,GEOGCS[&quot;NAD83&quot;,DATUM[&quot;North_American_Datum_1983&quot;,SPHEROID[&quot;GRS 1980&quot;,6378137,298.2572221010002,AUTHORITY[&quot;EPSG&quot;,&quot;7019&quot;]],AUTHORITY[&quot;EPSG&quot;,&quot;6269&quot;]],PRIMEM[&quot;Greenwich&quot;,0],UNIT[&quot;degree&quot;,0.0174532925199433],AUTHORITY[&quot;EPSG&quot;,&quot;4269&quot;]],PROJECTION[&quot;Lambert_Conformal_Conic_2SP&quot;],PARAMETER[&quot;standard_parallel_1&quot;,49],PARAMETER[&quot;standard_parallel_2&quot;,77],PARAMETER[&quot;latitude_of_origin&quot;,63.390675],PARAMETER[&quot;central_meridian&quot;,-91.86666666666666],PARAMETER[&quot;false_easting&quot;,6200000],PARAMETER[&quot;false_northing&quot;,3000000],UNIT[&quot;metre&quot;,1,AUTHORITY[&quot;EPSG&quot;,&quot;9001&quot;]]]</SRS>
<GeoTransform>  5.6058900000000000e+06,  3.0000000000000000e+01,  0.0000000000000000e+00,  3.8967900000000000e+06,  0.0000000000000000e+00, -3.0000000000000000e+01</GeoTransform>
<Metadata>
<MDI key=”AREA_OR_POINT”>Area</MDI>
</Metadata>
<VRTRasterBand dataType=”Float32″ band=”1″>
<Metadata>
<MDI key=”COLOR_TABLE_RULE_RGB_0″>0.000000e+00 4.821725e+01 255 255 0 0 255 0</MDI>
<MDI key=”COLOR_TABLE_RULE_RGB_1″>4.821725e+01 9.643450e+01 0 255 0 0 255 255</MDI>
<MDI key=”COLOR_TABLE_RULE_RGB_2″>9.643450e+01 1.446518e+02 0 255 255 0 0 255</MDI>
<MDI key=”COLOR_TABLE_RULE_RGB_3″>1.446518e+02 1.928690e+02 0 0 255 255 0 255</MDI>
<MDI key=”COLOR_TABLE_RULE_RGB_4″>1.928690e+02 2.410863e+02 255 0 255 255 0 0</MDI>
<MDI key=”COLOR_TABLE_RULES_COUNT”>5</MDI>
</Metadata>
<NoDataValue>-9.99900000000000E+03</NoDataValue>
<ColorInterp>Gray</ColorInterp>
<SimpleSource>
<SourceFilename relativeToVRT=”1″>cded_30m.tif</SourceFilename>
<SourceBand>1</SourceBand>
<SourceProperties RasterXSize=”5529″ RasterYSize=”4607″ DataType=”Float32″ BlockXSize=”5529″ BlockYSize=”1″ />
<SrcRect xOff=”0″ yOff=”0″ xSize=”5529″ ySize=”4607″ />
<DstRect xOff=”0″ yOff=”0″ xSize=”5529″ ySize=”4607″ />
</SimpleSource>
</VRTRasterBand>
</VRTDataset>

PostGIS 2.0svn Build & Install for Raster Support on Virtual Machine install of Fedora 16

Posted on by

Just a quick overview of how to build & install the great raster features of PostGIS 2.0 svn. I am aiming at using this data type in my current work of dem interpolation & raster derivatives based on Victorial Island in Canada’s Arctic.

As always please let me know what you think.

Build targets Fedora 16 & Postgresql 9.1 on a 32 bit machine. I am running this as a virtual machine using Virtualbox. I am using a virtual machine for the features it offers of cloning & snap shoting.

Assumptions:

  • localhost connection desired for database work
  • fresh install Fedora 16 with no spatial packages installed

We will be using YUM to supply all spatial components. This might not be ideal for the latest packages but will serve to build a functional postgis 2.0 with Raster support.

  1. Open terminal & su root
  2. install spatial components
    1. QGIS will install many of the spatial dependencies & will be used to review rasters from postgis2.0
      1. yum install qgis.i686 qgis-python.i686 qgis-devel.i686
      2. yum install proj-devel.i686 geos-devel.i686 gdal-devel.i686 libxml2.i686 libxml2-devel.i686 dblatex.noarch ImageMagick.i686 ImageMagick-devel.i686
  3. install postgresql database
    1. yum install postgresql.i686 postgresql-server.i686 postgresql-libs.i686 postgresql-contrib.i686 postgresql-devel.i686 pgadmin3
  4. setup database
    1. postgresql-setup initdb
  5. start postgresql
    1. service postgresql start
  6. enable service to start when machine starts
    1. systemctl enable postgresql.service
  7. su postgres -c psql
    1. ALTER USER postgres WITH PASSWORD ‘your password for postgres‘;
    2. \q
  8. edit pg_hba.conf file for localhost connections
    1. leafpad /var/lib/pgsql/data/pg_hba.conf
      1. change “method” column values to “md5” & save
    2. service postgresql restart
  9. open pgAdmin 3
    1. File – Add Server
      1. name = test
      2. host = localhost
      3. username = postgres
      4. password = what was set in step 7.1
    2. Right click Login Roles – New Login Role ( create another user to connect to database)
      1. Properties Tab: role name = enter a name
      2. Definition Tab: password = enter password
      3. Role Privileges Tab: click Superuser
  10. create database to install postgis onto. (-U = user to connect to database. Use name from step 9.2)
    1. createdb -U richard test_postgis
  11. cd into directory to download PostGIS2.0svn in & extract
    1. wget http://postgis.refractions.net/download/postgis-2.0.0SVN.tar.gz
    2. tar -zxvf downloaded file from above
    3. cd into extracted directory
    4. ./configure –with-raster
      1. check for warnings with ./configure –with-raster | grep -i warnings
      2. install any dependencies that are required. Often the development packages are required.
    5. Make
      1. ensure this is output at end “PostGIS was built successfully. Ready to install.”
    6. make install
  12. enable postgis onto created database from step 10
    1. psql -d test_postgis -U richard -f /usr/share/pgsql/contrib/postgis-2.0/postgis.sql
    2. psql -d test_postgis -U richard -f /usr/share/pgsql/contrib/postgis-2.0/spatial_ref_sys.sql
    3. psql -d test_postgis -U richard -f /usr/share/pgsql/contrib/postgis-2.0/rtpostgis.sql

 

The spatial database has been created with the PostGIS 2.0svn Raster support. Please refer to the following documentation on the functions present for the raster type:

 

In a follow up article I will be talking about loading raster data into the postgis database using the GDAL toolkit.

 

 

Mapping Ottawa Transit Stops with GTFS data & Fusion Tables

Posted on by

Back to some spatial goodness for the holidays.

We are going to map transit stops using GTFS data & Fusion Tables for the Ottawa (Canada) OC Transpo Service . I didn’t look for duplicates or other data scrubbing within the data. This is just a quick work flow to get some data in, transform it to meet your needs and map it using Fusion Tables.

I’ve created it as a PDF, link below.

Mapping OC Transpo stops with GTFS and Fusion Tables + Some Linux

Have fun and let me know what you think!

Speaking with GIS students on Working for the Government

Posted on by

Had the opportunity this week to speak to a class of GIS Application Specialist students at Fleming College in Lindsay. A big thanks to my former professor Noreen for being so receptive. I had signed up for a student mentor ship program in the summer with the college and have yet to hear back but pitched my professor with a presentation on “getting into the government”. I’ve been very fortunate to have worked in spatial depts for the Provincial & Federal government for the past 3 years. My motivation for the talk was to give back the place that helped started my career. Presentation covers 3 key areas: i) reality of geospatial govt jobs and how to get “in” ii) the critical importance of networking & iii) how open source tools can be a competitive advantage (deeper understanding of tools and processes).

So a minimal of slides in lighting talk format, I wanted to give plenty of question time. The students were great, once the questions got going. Even through they have just started the program, many were eager to query the coming job search they’d face at the end of next June. All in all, it was a fantastic time to engage with the students & touch base with so many of professors.

Please feel free to comment on the slides here (pdf format). I’ve released them under a create commons license.

Get Working for the Government Talk

Upcoming items

Posted on by

I will be following with the “Understanding Your Data” series shortly.

Just two items below that have been taking up my spare time recently. Will be posting pics & short descriptions of them. Have found the Arduino to be such a neat platform, working on the Processing language and working along with the electronic tutorials. Very much into data loggers at the moment, huge growth coming in this sector. Think “internet of things”, but being able to capture “own” data and slice/dice how I want is really fascinating.

  • I’ve borrowed an Em50R data logger from Decagon Devices from a friend for the next bit to play around with. This is a 902 MHz ISM 5 channel logger with the ECH2O-TE (Water Content, EC & Temperature sensors). Planning on installing on Linux with Wine. Photos to follow soon.
  • For my Arduino, I’ve picked up the GPS Logger shield kit & 20 Channel GPS Engine from the great folks at the Robotshop in Montreal. It’s been a long time since I’ve soldered and it shows..lol…but I’ve managed to get the kit working and had GPS signals inside my 14th floor condo (by the window) after initializing outside. Will have to buy another logger shield kit, as I soldered the Arduino connection pins the wrong way,  no fault of the great documentation from ladyada.net

 

 

 

 

Understanding your data: Part 2

Posted on by

We left off discussing the effects that different filters can have on slope values & alter the elevation data for future derivatives.

Here we will be examining the outcomes of filling a raster.

This is the process of applying a filling algorithm to the raster to remove pits/depressions (a single cell or areas with no drainage). In most cases, pits are the product of data errors/artifacts and/or the way in which the raster was assembled. Most natural landscapes do not contain pits (there are many landscapes that do such as karst landscapes, but you would/should be aware of this). So removing pits is about cleaning up your data and again generating as smooth as surface as possible for flow routing. The D8 ( 8 direction) flow direction outputs will not be correct if pits/depressions are not filled. The D8 is the most commonly implemented flow direction algorithm and used in several main stream GIS’s & open source tools There are numerous other flow direction implementations that use multiple directions i.e FD8, D-Inf, DEMON.

Getting back to the filling process, there are numerous ways this can be accomplished from assigning a percent slope that each cell must maintain to route flow to its neighbor to specifying a Z threshold to fill pits up to.

It’s important to review the filling algorithms present in your software package; read the documentation & search on-line for the authors papers describing their process.

3 examples below showing nearly the same area as Part 1 (PEI, Lidar 1.5m) elevation rasters that have been filled using SAGA. The filled rasters where created using the default tool options. Graphics below created by subtracting the filled & original rasters.

What’s important is to note the spatial variability in location & extent of areas that where filled amongst the different algorithms. View this as a relative overview, I have not provided legends because I feel it is more important to view the differences between algorithms and not imply which is “better”.

As I’ve said before, you will need to test & test. Creating a difference map is a great start but viewing the flow direction raster should be the first step to ensure flow can be routed correctly & the output makes sense. In the case of the D8 flow dir, if you have more than 8 directions you have a problem & need to fix your surface. A single pass of a filling tool may not be enough but this will depend on the tool, data & goals. Again it is importance to create hill shades and other visual outputs (profile plots comparing the original to modified) to help understand the surface alterations.

 

Next up I will be showing how to create Kernel Density Plots in R to access how slope/curvatures change as you apply different filter sizes.

Please drop me a note to tell me what you think.

Happy Thanksgiving weekend!