[R-sig-Geo] Methodology to compare crop maps
Thiago Veloso
thi_veloso at yahoo.com.br
Tue Jun 21 13:17:37 CEST 2011
Robert,
Thank you very much. Your tips helped me to achieve my goals! Thanks for your patience!
Best wishes,
Thiago.
--- On Thu, 16/6/11, Robert J. Hijmans <r.hijmans at gmail.com> wrote:
> From: Robert J. Hijmans <r.hijmans at gmail.com>
> Subject: Re: [R-sig-Geo] Methodology to compare crop maps
> To: "Thiago Veloso" <thi_veloso at yahoo.com.br>
> Cc: r-sig-geo at r-project.org
> Date: Thursday, 16 June, 2011, 8:53
> Dear Thiago,
>
> That's why I suggested:
>
> > cane_sa[cane_sa < 16] <- NA
> > veg_sa <- crop(vegmap,south_america)
> > newveg_sa <- cover(cane_sa, veg_sa)
> > new_vegmap <- merge(newveg_sa, vegmap)
>
> i.e. set all values you do not want to NA, and then use
> 'cover' to
> replace NA values in the new raster with the values of the
> old raster,
> and then use merge.
>
> Best, Robert
>
>
> On Thu, Jun 16, 2011 at 4:11 AM, Thiago Veloso <thi_veloso at yahoo.com.br>
> wrote:
> > Dear Robert,
> >
> > Sorry for not stating my question concisely. Still,
> you captured what was wrong in my code: that was exactly the
> threshold of the existence of the crop. Raising the limit
> helped to alleviate the crop existence map at a coarser
> resolution.
> >
> > Following the code, I experienced another
> difficulty. The merge function replaces the whole contents
> from one map to the other. Is there any way to impose
> conditions? Something like "merge if pixel value is greater
> than..."?
> >
> > Best wishes,
> >
> > Thiago Veloso.
> >
> > --- On Sat, 11/6/11, Robert J. Hijmans <r.hijmans at gmail.com>
> wrote:
> >
> >> From: Robert J. Hijmans <r.hijmans at gmail.com>
> >> Subject: Re: [R-sig-Geo] Methodology to compare
> crop maps
> >> To: "Thiago Veloso" <thi_veloso at yahoo.com.br>
> >> Cc: r-sig-geo at r-project.org
> >> Date: Saturday, 11 June, 2011, 1:27
> >> Hi Thiago,
> >>
> >> You are not saying _what_ goes wrong. That makes
> it
> >> difficult to help.
> >> I am guessing that you are missing some of the
> steps below.
> >> You do:
> >>
> >> cane_sa[cane_sa>0]<-16
> >>
> >> # I do not know your exact objectives but it seems
> strange
> >> to use 0 here,
> >> # as even a cell that is covered for 1/600 of its
> area with
> >> cane would
> >> be classified as cane.
> >> # but irrespective of the threshold you choose, I
> think you
> >> need:
> >>
> >> cane_sa[cane_sa < 16] <- NA
> >> veg_sa <- crop(vegmap,south_america)
> >> newveg_sa <- cover(cane_sa, veg_sa)
> >>
> >> # and then continue with your code again:
> >>
> >> new_vegmap <- merge(newveg_sa, vegmap)
> >>
> >> Robert
> >>
> >> On Thu, Jun 9, 2011 at 3:51 PM, Thiago Veloso
> <thi_veloso at yahoo.com.br>
> >> wrote:
> >> > 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