Welcome to my Open Notebook

This is an Open Notebook with Selected Content - Delayed. All content is licenced with CC-BY. Find out more Here.

ONS-SCD.png

Ripping data from georeferenced images with QGIS

Introduction

Drought Declarations are made by government in Australia to determine when and where a drought is occuring. The Declarations are linked to agricultural and social support policies. We are working on a climatic drought index that should correlate with these Drought Declarations. The New South Wales Government has provided us with spatial data from 1986 but they also have a graphical visualisation available for earlier times, especially noteworthy is the extreme 1982-1983 drought. This post is a document of the process I will try to derive spatial data from the image.

Methods

Data

Download the drought maps. The current maps are at: http://www.dpi.nsw.gov.au/content/agriculture/emergency/drought/planning-archive/climate/advance-retreat

But there is a higher resolution map archived at: http://pandora.nla.gov.au/pan/120345/20120529-0000/advance-retreat-drought-map-april-2011.pdf

drt-nswdpi-raw

Results

Step 1 - load the reference layer

  • Load up the NSW boundaries by selecting "Add PostGIS Layers" > Connect > gislibrary
  • Then select a map layer of NSW

Step 2 - set up the raster data into the Georeferencer

  • Launch the Georeferencer from Raster > Georeferencer > Georeferencer.
  • If you do not see that menu item, you will need to enable the Georeferencer GDAL plugin from Plugins > Manage and install Plugins > Installed.
  • Zoom in on the month of interest
  • Now click on the Add Point button on the toolbar and select an easily identifiable location on the image. Corners, intersections, poles etc. make good control points.
  • Once you click on the image at a control point location, you will see a pop-up asking you to enter map coordinates. Click the button From map canvas.
  • Find the same location in your reference layer and click there. The coordinates are auto-populated from your click on the map canvas. Click Ok. Similarly, choose at least 4 points on the image and add their coordinates from the reference layer.
  • Now go to Settings > Transformation settings (yellow cog wheel).
  • Choose the settings as shown below.
  • Transformation type = Thin Plate Spline
  • output raster = someting
  • target EPS = 4283 (match what your reference layer is)
  • Make sure the Load in QGIS when done button is checked. Click OK.
  • Back in the Georeferencer window, go to File > Start georeferencing. This will start the process of warping the image using the GCPs and creating the target raster.

Step 3 - rip the data to a shapefile

  • create a new shp (Layer > New)
  • type = ploygon
  • Specify CRS = EPSG:4283
  • attribute = drought, text
  • save to a shp like 'drt198610rip' and OK
  • right-click, toggle editing, add new feature
  • save layer edits

drt-edit-shapefile

Results

This work

drt198610rip

Compared to NSW DPI drought declarations data

drt198610rip2

Conclusion

The process produces data that is very prone to imprecision. We may be able to use some extra data from before 1986 though.

Posted in  drought awap grids data operation


drought-awap-grids-notes-re-extract-point

Background

  • The old DROUGHT-BOM-GRIDS data came from a Barnes IDW algorithm, at 25km (actually 0.25 decimal degrees) resolution
  • 1890-2008
  • Going forward the AWAP grids are Preferred
  • this is a quick and dirty approach

Methods

  • AWAP grids are 5km (0.05 dd) but this will take longer to process and is probably overkill for monthly total rainfall when focus is on drought
  • For the centroid point of each 0.25 dd extract the data value (might be better to take mean of sub-cells within each, plan for future)
  • compare this with old BOM GRIDS

Results

images/qc_awap_totals_200001.png

compare with bom grids for Jan 2000

images/qc_awap_200001_vs_bomgrids.png

For a pixel (west wyalong my old home town)

images/qc_awap_1900_1908_vs_bomgrids_west_wyalong.png

Posted in  drought awap grids


colours-to-use-for-maps

Purpose

Selected subset of all colour palettes

  • This is the agreed set of colour palettes to be used for maps and figures
  • They are all 1) Colour blind safe, 2) Computer screen safe and 3) printer safe
### Load the package or install if not present
if (!require("RColorBrewer")) {
install.packages("RColorBrewer")
library(RColorBrewer)
}
## Loading required package: RColorBrewer
par(mfrow = c(3,3))
for(col_i in c('YlGn','RdPu', 'PuRd', 'BrBG', 'RdBu', 'RdYlBu', 'Set3', 'Set1')){
  display.brewer.pal(n = 5, name = col_i)
}

  • A machine-readable approach to color specification is as hexadecimal triplets.
  • Here is how the RColorBrewer RdYlBu palette is actually stored:
for(col_i in c('YlGn','RdPu', 'PuRd', 'BrBG', 'RdBu', 'RdYlBu', 'Set3', 'Set1')){
  print(col_i)
  print(brewer.pal(n = 5, name = col_i))
}
## [1] "YlGn"
## [1] "#FFFFCC" "#C2E699" "#78C679" "#31A354" "#006837"
## [1] "RdPu"
## [1] "#FEEBE2" "#FBB4B9" "#F768A1" "#C51B8A" "#7A0177"
## [1] "PuRd"
## [1] "#F1EEF6" "#D7B5D8" "#DF65B0" "#DD1C77" "#980043"
## [1] "BrBG"
## [1] "#A6611A" "#DFC27D" "#F5F5F5" "#80CDC1" "#018571"
## [1] "RdBu"
## [1] "#CA0020" "#F4A582" "#F7F7F7" "#92C5DE" "#0571B0"
## [1] "RdYlBu"
## [1] "#D7191C" "#FDAE61" "#FFFFBF" "#ABD9E9" "#2C7BB6"
## [1] "Set3"
## [1] "#8DD3C7" "#FFFFB3" "#BEBADA" "#FB8072" "#80B1D3"
## [1] "Set1"
## [1] "#E41A1C" "#377EB8" "#4DAF4A" "#984EA3" "#FF7F00"
# or for more levels
brewer.pal(n = 10, name = "RdYlBu")
##  [1] "#A50026" "#D73027" "#F46D43" "#FDAE61" "#FEE090" "#E0F3F8" "#ABD9E9"
##  [8] "#74ADD1" "#4575B4" "#313695"
The leading # is just there by convention. Parse the hexadecimal string like so: #rrggbb, where rr, gg, and bb refer to color intensity in the red, green, and blue channels, respectively. Each is specified as a two-digit base 16 number, which is the meaning of "hexadecimal" (or "hex" for short). Here's a table relating base 16 numbers to the beloved base 10 system.

All colour palettes

### Show all the colour schemes available
par(cex = .6)
display.brewer.all()

### Set the display a 2 by 2 grid
par(mfrow=c(2,2))


### Generate random data matrix
rand.data <- replicate(8,rnorm(100,100,sd=1.5))

### Draw a box plot, with each box coloured by the 'Set3' palette
boxplot(rand.data,col=brewer.pal(8,"Set3"))

### Draw plot of counts coloured by the 'Set3' pallatte
br.range <- seq(min(rand.data),max(rand.data),length.out=10)
results <- sapply(1:ncol(rand.data),function(x) hist(rand.data[,x],plot=F,br=br.range)$counts)
plot(x=br.range,ylim=range(results),type="n",ylab="Counts")
cols <- brewer.pal(8,"Set3")
lapply(1:ncol(results),function(x) lines(results[,x],col=cols[x],lwd=3))

### Draw a bar chart
table.data <- table(round(rand.data))
cols <- colorRampPalette(brewer.pal(8,"Dark2"))(length(table.data))
barplot(table.data,col=cols)

Other reference material

Posted in  exploratory data analysis


Australian Postal Areas In Geographical Vs Projected Coordinates

Download

'name:poa06-area-lambert'
setwd("~/projects/POA_centroids/POA2006_centroids")
library(swishdbtools)
ch <- connect2postgres2("delphe")

fout_geo=dbGetQuery(ch,
'select poa_2006,
  st_area(st_transform(the_geom, 3112))/1000000 as Geoscience_Australia_Lambert_area_km2,
st_x(st_centroid(st_transform(the_geom,3112))) as geocentx,
st_y(st_centroid(st_transform(the_geom,3112))) as geocenty
from abs_poa.auspoa06')
str(fout_geo)
sum(fout_geo$geoscience_australia_lambert_area_km2)
write.table(fout_geo,'data_derived/auspoa06_geocentroids_lambert_20160624.csv',
            row.names=F, sep=',')

plot(fout_geo[,3:4])
head(fout_geo)
nrow(fout_geo)
2507

Posted in  spatial


spatial-lag-and-timeseries-model-with-nmmaps-UPDATE

Posted in  spatial dependence