Tag Archives: gis

Openstreetmap Notes in Ottawa

Posted on by

kicking off the long weekend here in canada in style. friday with some openstreetmap! openstreetmap added a new feature called notes that allows anyone to create an annotated message of a feature/item that needs to be added/fixed etc without needing an account. basically someone with an account attempts to add the note. great way to lower the barrier to reach a larger community that might not want to invest in learning how to use an editor & tagging. i can certainly see this being the case for small business owners that just want to get their location on the map.

anyhow, there’s a great stats page listing contributors here. i’ve wanted to query ottawa & see how many notes there have been.

bbox_ottawa_notes

used to qgis dev (which is amazing! mac dev builds) to import openstreetmap, creating a bounding box polygon & extracting the bbox coordinates. this could definitely be refined. in a nut shell you pass the bbox to the osm notes api & it returns xml (default) or several other formats. the start represents the only note in the region in gpx.

api call:

curl http://api.openstreetmap.org/api/0.6/notes?bbox=-76.30,45.09,-74.96,45.62?closed=-1

single note details — red star in graphic:

<?xml version=”1.0″ encoding=”UTF-8″?> <osm version=”0.6″ generator=”OpenStreetMap server”> <note lon=”-75.6967664″ lat=”45.421061″> <id>3125</id> <url>http://api.openstreetmap.org/api/0.6/notes/3125</url&gt; <comment_url>http://api.openstreetmap.org/api/0.6/notes/3125/comment</comment_url&gt; <close_url>http://api.openstreetmap.org/api/0.6/notes/3125/close</close_url&gt; <date_created>2013-05-14 22:33:23 UTC</date_created> <status>open</status> <comments> <comment> <date>2013-05-14 22:33:24 UTC</date> <uid>1679</uid> <user>andrewpmk</user> <user_url>http://api.openstreetmap.org/user/andrewpmk</user_url&gt; <text>Ottawa light rail, construction will begin soon. Also missing major bus routes on transitway like #95.</text> <html>&lt;p&gt;Ottawa light rail, construction will begin soon. Also missing major bus routes on transitway like #95.&lt;/p&gt;</html> </comment> </comments> </note> </osm>

really need to prompt the notes features in my neck of the woods. it’s a great feature. next steps would be to take this into python & clean it up, gives me a good excuse to keep learning python!

enjoy,

richard

Installing QGIS Master (1.9) on Mac OSX

Posted on by

just a few links that made the process of installing qgis master aka 1.9 development (soon to be 2.0; sometime this summer) on a mac osx. i’m running lion 10.7 with mac book pro. honestly, these folks do such an amazing job & free to provide these frameworks & builds; so kudos!!!

from qgis page, you can install either current qgis 1.8 or dev 1.9.  regardless you will need the frameworks provided by KyngChaos, install the GDAL & GSL complete; add GRASS if needed. i installed the nightly dev (1.9) builds provided by larry shaffer of dakota cartography, you will need to install GDAL 1.10 complete in order for the build to fire up on your box.

seriously, it’s that easy. install the frameworks & download nightly zip & unpack & double click the dmg. to start qgis from command line, you will need to add an alias, i’ll let you figure that out:)

 

enjoy,

richard

 

newest free semi-national 1 arc sec srtm derived elevation data -> cdsm

Posted on by

don’t like integer cded? probably not alone in that regard but there’s a partial answer. at least if you need data below 60 deg north. enter the  canadian digital surface model (cdsm), produced by the folks out at centre for topographic information (cits) in sherbrooke quebec. it’s floating point data to start with & is the best publicly available elevation data in canada. i’m currently using it for lake basin delineation & have been very happy with the results.

technical details below:

0.75 Second Canadian Digital Surface Model (CDSM) Version 1.0:

The national 0.75-second (~20 m) CDSM consists of a derived product from the original 1-second (30 m) Shuttle Radar Topographic Mission (SRTM) digital surface model (DSM). The SRTM data were reprocessed from their original form as follows: 1) the data was void-filled using the Canadian Digital Elevation Model (CDEM) with the Delta Fill Surface method from Grohman et al. (2006); 2) the Vertical Datum was changed from EGM96 to CGVD28; 3) the data was projected to UTM using cubic convolution; 4) the data was smoothed using the Denoise algorithm from Sun et al. (2010); 5) the data was re-projected from WGS84 to NAD83 and aligned to the 0.75-second CDEM grid resolution using cubic convolution; and, 6) the waterbodies for which new elevations were calculated due to the change of vertical datum were reflattened.

you can download it from the new geogratis portal., which by the way is a massive improvement over the older version. if you need it in larger swaths or by nts index, use the ftp site for download.

colored hillshaded cdsm 30m created in saga gis to give you an idea of what the data set looks like.

cdsm_30m_hillshade

enjoy,

richard

fme 2013 so far — python focused

Posted on by

working with fme desktop 2013 sp1 for the past few couple of weeks for a clients research project. first off, the more time i spend with it, the more it’s value is apparent. yes there is certainly a learning curve but well worth it.

my work involves calculating lake basin watersheds for an entire lake feature (on-demand from user & at national scale), not a single pour point. how do that, download saga & look under the terrain analysis – hydrology – upslope area, the target area attribute is the input rasterized lake.  saga offers an api using saga_cmd to access functions, so how to preprocess lakes ie buffer, generalize, reproject or any number of modifications & interface that newly lake feature to create basins? python is well intergrated in fme, i tried the pythoncaller but have not figured out how to write out single feature as temp shapefile to be read in by saga; had issues with the universalwriter & have filed case. will update when they get back in touch because the ability to write out a temp shapefile or other mid workflow stream, process with someothe software & reread back into the workflow stream is huge. so my approach has been to use scripted python parameters to setup the initial directory structure & covert the elevation geotiff to saga sgrd format. waterbodies get read in, processed & written out as individual lakes where a python shutdown script reads the lakes directory & calculates the lake basin for each file.

i’m sure there are better ways to accomplish this but i’ve only been digging into the python aspect for the last few days. it offers a lot of flexibility!

posted the initial python code for fme shutdown script to interface with saga. hope this gets you going in fme or you can use the example outside of fme, your choice:)

 

enjoy,

richard

Windows Install Postgresql + Postgis

Posted on by

Overview:
Installing Postgresql & Postgis in a Windows setting. These are some notes I’ve prepared for a recent Windows 7 x64 install.

Download installer (86-x64 9.1.*) from EnterpriseDB from http://www.enterprisedb.com/products-services-training/pgdownload#windows

Information

  1. http://www.postgresql.org/
  2. http://postgis.refractions.net/
  3. http://postgis.refractions.net/documentation/manual-2.0/

Installation

  1. Double click executable
  2. Next
  3. Accept default installation (C:\Program Files\PostgreSQL\9.1) directory unless guidelines oppose
  4. Accept default data directory (C:\Program Files\PostgreSQL\9.1\data)
  5. enter admin password for “postgres” user
  6. Next – Next
  7. Select “Postgresql 9.1 (x64) on port 5432
  8. Expand & Select Spatial Extensions – Postgis 2.0
  9. Next – Next
  10. Accept Postgis license
  11. Click “Create spatial database”
  12. Next
  13. user name = postgres
    1. password = pass word from step 5
    2. port = 5432
  14. database name = db_gis (used in example)
  15. Click yes to register GDAL_Data environmental variable

To verify installation

  1. Start Menu – pgAdmin III
  2. Double click on “PostgreSQL 9.1 (localhost:5432)
  3. Expand Databases
  4. select db_gis
  5. expand Extensions (if postgis & postgis_topology present db is installed correct with spatial extensions)

Import Shapefile into database

  1. Open pgAdmin III
  2. select database
  3. Plugins – PostGIS shapefile & dbf loader 2.0
  4. Click “Add File” & browse to shapefile location
  5. Click Import
  6. Review output log window

Interacting with Database & Data

  1. Open pgAdmin III
  2. select database
  3. Click “SQL” magnifier icon
  1. Sql query editor is used to write ad-hoc commands or load pre-written *.sql files
  2. Example ad-hoc commands Show geometry information of loaded data tables (data must be loaded prior)
    1. select * from geometry_columns
    2. select * from [table_name] limit 5 (show all columns but limit to 5 rows. Table name is imported shapefile name excluding the *.shp)
  3. File – Open. To browse for *.sql file to execute.

Exporting Table Data to Shapefile

Note: This applies only to tables with geometry column.

  1. Open pgAdmin III
  2. select database
  3. Plugins – PostGIS shapefile & dbf loader 2.0
  4. Click Export tab
  5. Click Add Table & select spatial table
  6. Click Export & select location to save shapfile
  7. Review log window output

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>

Using GRASS within R: Part 1

Posted on by

Recently I’ve been doing work on determining oriented landscapes on Eastern Victoria Island. If you’re familiar with the Canadian NTS index I’m working primarily in NTS 77E for testing. Surficial geologists are very interested in processes that assist in explaining genesis. The landscape from a regional perspective can be broken down into two broad categories: 1)oriented & 2) non-oriented. Oriented really means stream line like landforms (think drumlins, flutes, eskers) that have an azimuth & can be grouped into fields. All other terrain, primarily hummocky are grouped into non-oriented.

Using the artefact rich CDED dem, I have a work-flow that identifies larger oriented features. Now the trick is to figure out the orientation that each feature exhibits & further begin to group fields that share similar orientations. I have not mentioned direction because that implies an angle & “flow” direction. By finding the streamline forms, the surficial geologist can determining from the drumlin fields what the direction of ice travel was. So the orientations I calculate will need to be adjusted after consulting with them. Not a big issue really, once you have the algorithm down for this.

 

I have been using GRASS within R. The spgrass6 package is great as it allows the creation of temporary grass location & mapsets for use. For this to work you need to have GRASS installed.

## load libraries
library(rgdal)
library(spgrass6)
library(raster)

## create tmp grass location
initGRASS(gisBase=”/usr/local/grass-6.4.3svn/”,home=tempdir(),override=T)
# gmeta<-gmeta6()

## change to permanent mapset to apply projection change for region
execGRASS(“g.mapset”,flags=c(“c”),parameters=list(mapset=”PERMANENT”))
## convert default coordinate system x,y to epsg code 3347 (can be altered later)
execGRASS(“g.proj”,flags=c(“c”),epsg=as.integer(3347))

## import clusters majority. “o” flag uses location proj. location & raster epsg must be same.
execGRASS(“r.in.gdal”,flags=c(“overwrite”,”o”),parameters=list(input=”../Tmp/saga_walk_through/Clusters_gdal.sdat”,output=”clustersInput”))

## review env settings
execGRASS(“g.gisenv”,flags=”s”)

# ## list dems
execGRASS(“g.list”,parameters=list(type=”rast”))

## set region to import clusters majority raster
execGRASS(“g.region”,flags=c(“p”,”a”),parameters=list(rast=”clustersInput”))

## apply mode 3*3 to clusters
execGRASS(“r.neighbors”,flags=c(“overwrite”),parameters=list(input=”clustersInput”,output=”clustersMode3″,method=”mode”))

### write vect to R obj for viewing
## read org clusters mode 3*3 into R
rastMode3Clusters<-readRAST6(“clustersMode3″)
## convert the spatialgriddataframe to raster object
rMode3Clusters<-raster(rastMode3Clusters)

## plot
plot(rMode3Clusters)

## write to disk for future reference. INT1U == BYTE
writeRaster(rMode3Clusters,filename=”Tmp/orgClustersMode3.tif”,format=”GTiff”,datatype=”INT1U”,NAFlag=as.integer(255),overwrite=TRUE)

 

I hope this gives you a quick idea of using GRASS within R.

I will be posting more on this matter in the coming days & will get a bit more technical with it.

 

Cheers,

Richard

OSGEO Ottawa Update

Posted on by

Hope everyone has had a great Holiday period!!

I’ve recently been asked to be the chair for our local Osgeo Ottawa group. Must say I am very honoured to have been asked:) Very excited to be working in this role and working with everyone in our community to keep the great momentum & friendly atmosphere that others before me have fostered. If you have any ideas, please post a comment.

Just proposed our first New Yrs meetings here. Really hope you can get out!

Looking forward to seeing everyone very soon.

Cheers,

Richard