Tag Archives: R

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

 

 

 

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

 

 

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