Warning: Declaration of syntax_plugin_alertbox::handle($match, $state, $pos, &$handler) should be compatible with DokuWiki_Syntax_Plugin::handle($match, $state, $pos, Doku_Handler $handler) in /homepages/38/d684553722/htdocs/clickandbuilds/spatialecology_wiki/lib/plugins/alertbox/syntax.php on line 22

Warning: Declaration of syntax_plugin_alertbox::render($mode, &$renderer, $data) should be compatible with DokuWiki_Syntax_Plugin::render($format, Doku_Renderer $renderer, $data) in /homepages/38/d684553722/htdocs/clickandbuilds/spatialecology_wiki/lib/plugins/alertbox/syntax.php on line 45

Warning: Declaration of syntax_plugin_bliki::handle($match, $state, $pos, &$handler) should be compatible with DokuWiki_Syntax_Plugin::handle($match, $state, $pos, Doku_Handler $handler) in /homepages/38/d684553722/htdocs/clickandbuilds/spatialecology_wiki/lib/plugins/bliki/syntax.php on line 75

Warning: Declaration of syntax_plugin_bliki::render($mode, &$renderer, $data) should be compatible with DokuWiki_Syntax_Plugin::render($format, Doku_Renderer $renderer, $data) in /homepages/38/d684553722/htdocs/clickandbuilds/spatialecology_wiki/lib/plugins/bliki/syntax.php on line 379
wikidk:dk10agrinut []

User Tools

Site Tools


Nutrient losses

Script done by Toke Panduro, KU, tepp@life.ku.dk, Maria Konrad, AU, mthk@dmu.dk, Line Block Hansen, AU, lbc@dmu.dk


Nutrient losses to the water environment is a major problem in many European and non-European countries. Excessive losses from agricultural areas enhance eutrophication processes with loss of biodiversity and groundwater contamination as a consequence. This again mean decreased economic values of fishery and recreational value.

General framework of this analysis

In this project we have two data files, an excel file containing information about nutrient losses from agricultural areas and a shape file with the same areas. The purpose of the project is to surf around getting to know a bit more about the different programs, GRASS, QGIS, R and the BASH language. We are going to explore our data, reading them into GRASS and R and carry out some simple statistical exercises.

Project objectives

To load data into GRASS and R. Carry out preliminary analysis of the data and plot different results in QGIS and R.


Excel file

We have an excel file containing information about nutrient losses, retention and application of nitrogen, application of fertilizer and manure, yield type and other variables for 10,682 smaller agricultural areas. These areas are aggregated into 25 catchment areas around Odense Fjord.

Vector file

We have two vector files in shape format. One containing the 25 catchment areas (polygons) and one with the 10,682 smaller agricultural areas (points)

Vectorfile 1: Layer name: 25oplinfo Geometry: Polygon Feature Count: 25 Extent: (568411.591663, 6106805.020500) - (604461.949080, 6159975.773157) Projection [“ETRS1989UTMZone_32N”, UNIT[“Meter”,1.0]]

Vectorfile 2: Layer name: Markblok Geometry: Point Feature Count: 10682 Extent: (568654.170939, 6106822.772875) - (604447.698678, 6159673.938541) Projection [“ETRS1989UTMZone32N”, UNIT[“Meter”,1.0]]

Finding info about the files: use ogrinfo
cd    ~/ost4sem/project/unidk2010/nutrients/input/
ogrinfo -al   ~/ost4sem/project/unidk2010/nutrients/input/25opl_info.shp | head -28
ogrinfo -al   ~/ost4sem/project/unidk2010/nutrients/input/Markblok.shp | head -28
Fancy loop for calling projection for all files 
for file in *.shp ; do ogrinfo -al $file    | grep PROJCS  ; done #-al list all features in all layers



We saved data from an USB stick to the desktop and loaded data from the desktop to input folder

INDIR=~/ost4sem/project/unidk2010/nutrients/input   # setting the input file directory which is necessary every time GRASS is opened 
OUTDIR=~/ost4sem/project/unidk2010/nutrients/output # setting the output file directory

Change directory to the directory where the file to move is located, we saved ours in desktop directory

cd ~/Desktop 
mv basedata.csv ~/ost4sem/project/unidk2010/nutrients/input/basedata.csv 
mv 25opl_info.dbf $INDIR/25opl_info.dbf
mv 25opl_info.prj $INDIR/25opl_info.prj
mv 25opl_info.sbn $INDIR/25opl_info.sbn
mv 25opl_info.sbx $INDIR/25opl_info.sbx
mv 25opl_info.shp $INDIR/25opl_info.shp
mv 25opl_info.shx $INDIR/25opl_info.shx
mv Markblok_25opl.dbf $INDIR/Markblok.dbf
mv Markblok_25opl.prj $INDIR/Markblok.prj
mv Markblok_25opl.sbn $INDIR/Markblok.sbn
mv Markblok_25opl.sbx $INDIR/Markblok.sbx
mv Markblok_25opl.shp $INDIR/Markblok.shp
mv Markblok_25opl.shx $INDIR/Markblok.shx


We work a bit in GRASS, picking out a specific catchment area and calculate shortest distances between points and the catchment boarders. The result can be seen in QGIS

Open grass in a preexisting mapset with a projection - in order to copy the excising projection

grass -text ~/ost4sem/grassdb/europe/PERMANENT 

Create a new location in GRASS to save our file with the projection obtained later. GRASS is by default creating a file called PERMANENT inside the odense directory.

rm -r ~/ost4sem/grassdb/odense

Import the shape file 25oplinfo.shp into the new grass location <code bash> v.in.ogr -e dsn=$INDIR/25oplinfo.shp output=opl25 location=odense –overwrite </code>

Exit GRASS to save the file in our working directory - doing this we have projected the projection from the excising data (in Permanent) to our file

rm -r ~/ost4sem/project/unidk2010/nutrients/grassdb/odense #Be sure that the directory not already exists!
mv ~/ost4sem/grassdb/odense ~/ost4sem/project/unidk2010/nutrients/grassdb/. #Move Odense to our GRASS database in the nutrient directory

Enter GRASS again to clean file. When importing shape files into GRASS they have to be cleaned due to errors created when imported

grass -text ~/ost4sem/project/unidk2010/nutrients/grassdb/odense/PERMANENT

Create cleaned file and an error file. It will only work if we remember the -e, making us able to extend the data, allowing errors when splitting and coupling the shape file

v.in.ogr -e dsn=$INDIR/25opl_info.shp output=opl25 --overwrite 
v.clean input=opl25 output=oplc25 tool=rmdupl,bpol error=myerroropl
g.remove vect=opl25 #Remove file not cleaned
v.in.ogr -e dsn=$INDIR/Markblok.shp output=Markblok1 --overwrite
v.clean input=Markblok1 output=Markblokc tool=rmdupl,bpol error=myerrormarkblok
g.remove vect=Markblok1 
qgis & #Open Quantum GIS to explore data

Import csv file to dbf file in GRASS

db.in.ogr dsn=$INDIR/basedata.csv output=basedata key=BLOKNR

Now the two files, Markblokc (the shapefile) and basedata (info file) is joined in QGIS using Vector/join attributes, and choose the relevant files and variables to join upon. The resulting file is a new shapefile named joinbasedatamarkblokc.shp placed in: ost4sem/project/unidk2010/nutrients/output

Load the joined file into GRASS and clean it
OUTDIR=~/ost4sem/project/unidk2010/nutrients/output # setting the output file directory
v.in.ogr -e dsn=$OUTDIR/join_basedata_markblokc.shp output=join_basedata_markblokc --overwrite 
v.clean input=join_basedata_markblokc output=join_basedata_markblokc_c tool=rmdupl,bpol error=myerrorjoin
g.remove vect=join_basedata_markblokc #Remove file not cleaned

We would like to find the catchment area with the greatest N loss, that is sum variable Nudv(kg) per catchment area

and find the a maximum value. Then we calculate the distance from markbloks in catchment area 25 to the other catchments using v.distance

v.out.ascii convert the GRASS vector file into a GRASS txt file, use the joinbasedatamarkblokcc file, seperating the colums with space (fs), only choose the two colums BNInr and Nudvkg and save it (>) into joinbasedatamarkblokc_c_table.txt

v.out.ascii   input=join_basedata_markblokc_c  fs=" "   columns=BNI_nr,Nudv_kg_,Opl_NS  > join_basedata_markblokc_c_table.txt

Do a loop to caculate the sum of Nudvkg for every catchment area (BNI_nr) and print the file with three relevant columns.

for (( clas=1 ; clas<=25 ; clas++)) ; do 
  awk  -v clas=$clas  '{ if ($4==clas) sum5=sum5+$5  }  END {print clas, sum5, $6 }' join_basedata_markblokc_c_table.txt   >>  join_basedata_markblokc_c_sum.txt 

Find the maximum value and print the maximum value and the BNInr ($1 = colum 1) and the OplNS (colum 3)

awk '{ if (NR>1) {if ($2>max) max= $2 }} END {print $1, $3, max }'  join_basedata_markblokc_c_sum.txt 

We found that catchment area no. 9 have the greatest N loss

Now we are going to estimate the distance between markbloks and catchment areas

First we clip the catchment area to focus on, catchment area 15 and the mark points from catchment area 9.

v.extract input=Markblokc output=Markblok9 where="BNI_nr=9"
v.extract input=oplc25 output=opl15 where="BNI_nr=15"

Now we calculate the distances from the markbloks in catchment area 9 to the boundaries in catchment area 15

v.db.addcol map=Markblok9@PERMANENT layer=1 'columns=dist double precision'#Add a new column for use in the v.distance 
v.info -c Markblok9
v.distance from=Markblok9 to=opl15 from_type=point to_type=boundary output=distance_9_15 upload=dist column=dist --overwrite

Now we do a loop to estimate distances from Markblok9 to all other catchment areas contained in oplc25

for (( clas=1 ; clas<=25 ; clas++)) ; do  # clas is the poligon ID
v.db.dropcol map=Markblok9  column=dist$clas
v.db.addcol map=Markblok9@PERMANENT layer=1 columns="dist$clas int"  #Add a new column for use in the v.distance
v.extract input=oplc25 output=opl$clas where="BNI_nr=$clas" --overwrite 
v.distance from=Markblok9 to=opl$clas from_type=point to_type=boundary upload=dist column=dist$clas --overwrite

Jubiiiiiiiiiiiiiiiiiiiiiiiiiiiii it worked…!!!!!!!!!!! :-P

Catchments areas and smaller agricultural areas

Model parametrization

Now we are changing to R, load data into R, make some regressions and plot different results. We start with making a new directory in Bash

mkdir ~/ost4sem/project/unidk2010/nutrients/R
mkdir ~/ost4sem/project/unidk2010/nutrients/R/output

Then we start up R


Read data into R and change values from factor to numeric values, finally saving data in the R directory

markdata <- read.dbf("~/ost4sem/project/unidk2010/nutrients/grassdb/odense/PERMANENT/dbf/join_basedata_markblokc_c.dbf", as.is=TRUE)

as.is=TRUE means that the function DO NOT convert characters to factors. This is important because characters can be converted to numeric numbers, factors can NOT

markdata$GLRareal_N <- as.numeric(markdata$GLRareal)
markdata$Nudv_kg_N <- as.numeric(markdata$Nudv_kg_)
markdata$BNI_nr_N <- as.numeric(markdata$BNI_nr)
markdata$jb_N <- as.numeric(markdata$jb)
markdata$N_han_N <- as.numeric(markdata$N_han)
markdata$N_hus_N <- as.numeric(markdata$N_hus)
markdata$HaGrovf_N <- as.numeric(markdata$HaGrovf)
markdata$HaGraesVed_N <- as.numeric(markdata$HaGraesVed)
markdata$HaKornVaar_N <- as.numeric(markdata$HaKornVaar)
markdata$HaKornVin_N <- as.numeric(markdata$HaKornVin)
markdata$HaAndet_N <- as.numeric(markdata$HaAndet)
markdata$HaBrakSkov_N <- as.numeric(markdata$HaBrakSkov)
save(markdata, file="~/ost4sem/project/unidk2010/nutrients/R/markdata.Rdata")

Linear regression

reg1 <- lm(markdata$Nudv_kg_N ~ markdata$GLRareal_N+markdata$jb_N+markdata$N_han_N+markdata$N_hus_N)
reg2 <- lm(markdata$Nudv_kg_N ~ markdata$GLRareal_N+markdata$N_han_N+markdata$N_hus_N)
reg8 <- lm(markdata$Nudv_kg_N ~ markdata$GLRareal_N+markdata$N_han_N+markdata$N_hus_N+markdata$HaGraesVed_N+markdata$HaKornVaar_N+markdata$HaBrakSkov_N+markdata$HaKornVin_N)

Print residuals and save in dataset markdata as the variable residual

markdata$residual <- reg8$residual

Now we make some plots…… for fun :-)

We make a png-function which creates a directory OUTR where plots are stored, and close the directory with dev.off

par(mfrow=c(2,1))   #De to plots kommer i samme billede
qqnorm(markdata$residual,main="QQ plot") #Two plots in the same picture
hist(markdata$residual, freq=TRUE, main="Test for normal distribution",xlab="Residuals") #Plot residuals against the frequency
dev.copy(png,filename="~/ost4sem/project/unidk2010/nutrients/R/Residuals.png",width = 1500, height = 1000, bg="white")

Now we make a loop: a regression for the 25 different catchment areas and save the estimate values in a 1:25 vector The estimate value from the regression is called by '$coefficient', we see that by wrigting str(summary(reg8))

for (i in c(1:25)){
temp.lm <- lm (subtmp$Nudv_kg_N~subtmp$GLRareal_N+subtmp$N_han_N+subtmp$N_hus_N+subtmp$HaGraesVed_N+subtmp$HaKornVaar_N+subtmp$HaBrakSkov_N+subtmp$HaKornVin_N)   
estim[i]=temp.lm$coefficients[3] #coeff henviser til estimate væriden og 3 til variabel 3 der er manure og gemmer i es1

Now we make a matrice with three different variables, estimate, r squared and intercept for the 25 catchments

for (i in c(1:25)){
temp.lm <- lm (subtmp$Nudv_kg_N~subtmp$GLRareal_N+subtmp$N_han_N+subtmp$N_hus_N+subtmp$HaGraesVed_N+subtmp$HaKornVaar_N+subtmp$HaBrakSkov_N+subtmp$HaKornVin_N)   
colnames(full.est)=c("estimate","R squared","intercept")
write.table(full.est,file="~/ost4sem/project/unidk2010/nutrients/R/Matrix.csv",row.names=TRUE,col.names=TRUE,sep=" ")

Doing some more fancy plots… We more libraries


Now we rename our data because the structure of data are going to be changed with the coordinate system


Look up the 'coordinates' and define geographic units

coordinates(fields)<-c("xcoor", "ycoor")


pal <- brewer.pal(10,"RdBu")
q10 <- classIntervals(fields$residual, n = 20, style="quantile")
q10Colours <- findColours(q10, pal)
plot(fields, axes=TRUE, col = q10Colours, pch =20)
legend("topleft", fill = attr(q10Colours, "palette"), legend=names(attr(q10Colours, "table")), bty="n")
dev.copy(png,filename="~/ost4sem/project/unidk2010/nutrients/R/map_of_Residuals.png",width = 1500, height = 1000, bg="white")

Creating and saving results from loop models

plot(full.est[,1],main="Manure",xlab="Catchment areas", ylab="Parameter estimates")
plot(full.est[,2],main="R square",xlab="Catchment areas", ylab="R square")
plot(full.est[,3], main="Intercept",xlab="Catchment areas", ylab="Parameter estimates")
dev.copy(png,filename="~/ost4sem/project/unidk2010/nutrients/R/loop_values.png",width = 1500, height = 1000, bg="white")

Fun merge which didn't work

# v.to.db -p mark option=coor > "\$INPUT"coordinates.txt
v.db.addcol map=mark layer=1 columns="xcoor double,ycoor double"
v.to.db  map=mark layer=1 option=coor column=xcoor,ycoor
R --no-save -q  << EOF
mark<-read.dbf("~/ost4sem/project/unidk2010/nutrients/grassdb/odense/fjord/dbf/mark.dbf") # 29var
base<-read.dbf("~/ost4sem/project/unidk2010/nutrients/grassdb/odense/PERMANENT/dbf/basedata.dbf") # 38
intable=merge(mark, base, by.mark = "comm.id", by.base = "comm.id") #64 var  not 67
write.table(intable,file="~/ost4sem/project/unidk2010/nutrients/input/intable_1.txt",sep="|",col.names=TRUE, row.name=FALSE)
### making a map from a text file #
cat ~/ost4sem/project/unidk2010/nutrients/input/intable_1.txt | v.in.ascii out=mypoints x=30 y=31 cat=1 skip=1 --overwrite
wikidk/dk10agrinut.txt · Last modified: 2021/01/20 20:36 (external edit)