[R-sig-Geo] Methodology to compare crop maps
Thiago Veloso
thi_veloso at yahoo.com.br
Fri Jun 10 00:51:52 CEST 2011
Robert, Rich and Luca,
Thank you very much for the suggestions.
Robert, I still haven't managed to implement all of your directions due to the burocracy to obtain the sugarcane raster (or shape) files from INPE.
In the meanwhile, I am trying to convert mcgill data (you're right about the units) from fraction of land covered by a crop to presence/absence. This is because my ultimate goal is to add a new plant functional type (PFT) to a map where 15 PFTs already exist (http://www.sage.wisc.edu/download/potveg/global_potveg.html). This way, cells where cane is present should contain the value "16".
It seemed to be a simple task, but what I am doing is changing the map substantially. Please take a peek at the code:
#Loading required packages
library(raster)
library(maptools)
#Loading SAGE vegetation map and McGill sugarcane map
vegmap<-raster("data/vegtype.nc")
cane<-raster("data/sugarcane_5min.nc")
> vegmap
class : RasterLayer
dimensions : 360, 720, 259200 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84
values : data/vegtype.nc
min : ?
max : ?
> cane
class : RasterLayer
dimensions : 2160, 4320, 9331200 (nrow, ncol, ncell)
resolution : 0.08333333, 0.08333334 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84
values : data/sugarcane_5min.nc
min : ?
max : ?
#Resampling McGill map to match SAGE's resolution
cane_coarser<-aggregate(cane,fact=6.00000024,fun=mean)
> cane_coarser
class : RasterLayer
dimensions : 360, 720, 259200 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84
values : /tmp/R_raster_tmp/raster_tmp_39924224710.grd
min value : 0
max value : 0.5680773
#Loading SA shapefile
south_america<-readShapePoly("shapes/southamerica.shp")
#Cropping sugar cane map
cane_sa<-crop(cane_coarser,south_america)
> cane_sa
class : RasterLayer
dimensions : 137, 93, 12741 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -81.49999, -34.99999, -56, 12.5 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84
values : in memory
min value : 0
max value : 0.5237272
#Below is an attempt to replace fraction with presence
cane_sa[cane_sa>0]<-16
#Check result
plot(cane_coarser)
#Adding cane map to SAGE vegetation map
new_vegmap<merge(vegmap,cane_sa)
#Writing out updated map to a netcdf file
if(require(ncdf)){
writeRaster(new_vegmap,filename="/home/thiago/IBIS/data/inpue/new_vegmap.nc",format="CDF",overwrite=TRUE)
}
Everything goes fine until cropping the resampled map. However, replacing grater than zero values results in a totally different map.
What am I doing wrong? Maybe spatial data frames are harder to deal with (for newbies like me), or maybe the massive amount of NA's and tiny numbers (0.000000e+00, 2.433714e-03 and etc) produced after agregation is the problem...
Thanks in advance and best wishes,
Thiago.
--- On Fri, 3/6/11, Robert Hijmans <r.hijmans at gmail.com> wrote:
> From: Robert Hijmans <r.hijmans at gmail.com>
> Subject: Re: [R-sig-Geo] Methodology to compare crop maps
> To: r-sig-geo at r-project.org
> Date: Friday, 3 June, 2011, 13:17
> > I am working with crops
> planted area maps from two distinct sources.
> > One of the maps is based on a maximum NDVI
> composition, and the other
> > map uses joint information from satellite and census
> to estimate the
> > planted area.
> >
> > Although the sources employ different methodologies
> to map the area
> > where the crop exists, the results should be
> comparable.
> >
> > After downloading the datasets, I have performed a
> visual inpection,
> > and they show reasonable agreement. However, I need a
> more robust
> > comparison method. Could anybody point out a
> methodology which allows
> > me to show the difference between both maps?
> >
> > Here is an example of each one of the maps:
> > http://www.geog.mcgill.ca/landuse/pub/Data/175crops2000/NetCDF/sugarcane_5min.nc.gz
> > (in netcdf) and http://www.dsr.inpe.br/laf/canasat/en/map.html (not
> > available to download directly, but I can get it in
> shapefile)
> >
>
>
> Thiago,
>
> I assume that the Brazilian data has a much higher spatial
> resolution than
> the mcgill data (that I think I am familiar with), and it
> probably has a
> different CRS. And I assume that you can get it as a the
> original raster
> file (and not as shapefile) for the Brazilian data. If I am
> not mistaken,
> the mcgill data has the fraction of land area covered by a
> crop. I assume
> that the Brazilan data is presence/absence. If so I would
> use the raster
> package and aggregate the Brazilian data to a cell size
> that is similar to
> the mcgill data (~9 km), computing the fraction of cells
> that have sugarcane
> (sum divided by the number of cells, make sure to handle NA
> values). Then
> use function projectRaster to transform the mcgill data to
> the same
> extent/resolution as the aggregated Brazilian data. Now you
> have two layers
> that you can compare in different ways.
>
> You can make plots, compute correlation, etc. Of course the
> p-values are no
> good because of spatial autocorrelation.
>
> library(raster)
> x <- y <- raster(nc=100, nr=100)
> x[] <- runif(ncell(r))
> y[] <- runif(ncell(r))
> plot(x, y)
> m <- lm(values(x), values(y))
> summary(m)
> abline(m)
>
> hist(x-y)
> plot(x-y)
> cor(values(x), values(y))
>
>
> Perhaps you want to treat your data as presence/absence
> (with presence being
> > 0 or some another threshold). These can then be easily
> compared with the
> crosstab function which returns, in this case, a confusion
> matrix which can
> be directly interpreted or used to compute some statistics
> from.
> crosstab(x>0, y>0)
>
> crosstab(x>0.5, y>0.5)
>
> And there surely are many other approaches possible, which
> is why I think
> that R is the way to go in this case: it is easy, flexible
> and fast.
>
> Robert
>
>
> --
> View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Methodology-to-compare-crop-maps-tp6431598p6435902.html
> Sent from the R-sig-geo mailing list archive at
> Nabble.com.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
More information about the R-sig-Geo
mailing list