

Dear RSIG colleagues,
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)
Thanks in advance,
Thiago Veloso.
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo


On Thu, 2 Jun 2011, Thiago Veloso wrote:
> 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.
> 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?
Thiago,
Please accept my suggestion that you're not using the best tool for this
job in the positive way I intend. This is a spatial analytical question best
answered by an analytical GIS rather than statistical software.
I've been using GRASS < http://grass.osgeo.org/> for a dozen or so years
and it will quickly answer your question as well as many others. Assuming
that your maps are grid (or raster) based a simple subtraction of one from
the other will show differences. You can also ask for differences in other
ways. The r.mapcalc() module is quite powerful. If your source data are in
vector format you can transform them to rasters for more powerful analysis.
HTH,
Rich
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo


On 06/02/2011 05:27 PM, Rich Shepard wrote:
>
> Please accept my suggestion that you're not using the best tool for this
> job in the positive way I intend. This is a spatial analytical question best
> answered by an analytical GIS rather than statistical software.
Well, yes... but the R raster package allows you to do map algebra and a host of
other raster functions as well: one may try that before switching to GIS.
It must be noted, however, that the raster package performance is poor when large
datasets are involved.
> I've been using GRASS < http://grass.osgeo.org/> for a dozen or so years
For a newbie I would suggest Quantum GIS coupled with its GRASS plugin instead:
same power, more friendly UI.
Regards,
Luca Morandini
http://www.lucamorandini.it
Regards,
Luca Morandini
http://www.lucamorandini.it_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo


In reply to this post by Thiago V. dos Santos
> 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 pvalues 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(xy)
plot(xy)
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


In reply to this post by Thiago V. dos Santos
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.433714e03 and etc) produced after agregation is the problem...
Thanks in advance and best wishes,
Thiago.
 On Fri, 3/6/11, Robert Hijmans < [hidden email]> wrote:
> From: Robert Hijmans < [hidden email]>
> Subject: Re: [RsigGeo] Methodology to compare crop maps
> To: [hidden email]
> 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
> pvalues 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(xy)
> plot(xy)
> 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://rsiggeo.2731867.n2.nabble.com/Methodologytocomparecropmapstp6431598p6435902.html> Sent from the Rsiggeo mailing list archive at
> Nabble.com.
>
> _______________________________________________
> RsigGeo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/rsiggeo>
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo


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 < [hidden email]> 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.433714e03 and etc) produced after agregation is the problem...
>
> Thanks in advance and best wishes,
>
> Thiago.
>
>  On Fri, 3/6/11, Robert Hijmans < [hidden email]> wrote:
>
>> From: Robert Hijmans < [hidden email]>
>> Subject: Re: [RsigGeo] Methodology to compare crop maps
>> To: [hidden email]
>> 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
>> pvalues 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(xy)
>> plot(xy)
>> 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://rsiggeo.2731867.n2.nabble.com/Methodologytocomparecropmapstp6431598p6435902.html>> Sent from the Rsiggeo mailing list archive at
>> Nabble.com.
>>
>> _______________________________________________
>> RsigGeo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/rsiggeo>>
>
>
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo


In reply to this post by Thiago V. dos Santos
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 < [hidden email]> wrote:
> From: Robert J. Hijmans < [hidden email]>
> Subject: Re: [RsigGeo] Methodology to compare crop maps
> To: "Thiago Veloso" < [hidden email]>
> Cc: [hidden email]
> 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 < [hidden email]>
> 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.433714e03 and etc) produced after agregation is the
> problem...
> >
> > Thanks in advance and best wishes,
> >
> > Thiago.
> >
> >  On Fri, 3/6/11, Robert Hijmans < [hidden email]>
> wrote:
> >
> >> From: Robert Hijmans < [hidden email]>
> >> Subject: Re: [RsigGeo] Methodology to compare
> crop maps
> >> To: [hidden email]
> >> 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
> >> pvalues 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(xy)
> >> plot(xy)
> >> 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://rsiggeo.2731867.n2.nabble.com/Methodologytocomparecropmapstp6431598p6435902.html> >> Sent from the Rsiggeo mailing list archive at
> >> Nabble.com.
> >>
> >> _______________________________________________
> >> RsigGeo mailing list
> >> [hidden email]
> >> https://stat.ethz.ch/mailman/listinfo/rsiggeo> >>
> >
> >
>
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo


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 < [hidden email]> 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 < [hidden email]> wrote:
>
>> From: Robert J. Hijmans < [hidden email]>
>> Subject: Re: [RsigGeo] Methodology to compare crop maps
>> To: "Thiago Veloso" < [hidden email]>
>> Cc: [hidden email]
>> 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 < [hidden email]>
>> 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.433714e03 and etc) produced after agregation is the
>> problem...
>> >
>> > Thanks in advance and best wishes,
>> >
>> > Thiago.
>> >
>> >  On Fri, 3/6/11, Robert Hijmans < [hidden email]>
>> wrote:
>> >
>> >> From: Robert Hijmans < [hidden email]>
>> >> Subject: Re: [RsigGeo] Methodology to compare
>> crop maps
>> >> To: [hidden email]
>> >> 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
>> >> pvalues 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(xy)
>> >> plot(xy)
>> >> 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://rsiggeo.2731867.n2.nabble.com/Methodologytocomparecropmapstp6431598p6435902.html>> >> Sent from the Rsiggeo mailing list archive at
>> >> Nabble.com.
>> >>
>> >> _______________________________________________
>> >> RsigGeo mailing list
>> >> [hidden email]
>> >> https://stat.ethz.ch/mailman/listinfo/rsiggeo>> >>
>> >
>> >
>>
>
>
>
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo


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 < [hidden email]> wrote:
> From: Robert J. Hijmans < [hidden email]>
> Subject: Re: [RsigGeo] Methodology to compare crop maps
> To: "Thiago Veloso" < [hidden email]>
> Cc: [hidden email]
> 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 < [hidden email]>
> 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 < [hidden email]>
> wrote:
> >
> >> From: Robert J. Hijmans < [hidden email]>
> >> Subject: Re: [RsigGeo] Methodology to compare
> crop maps
> >> To: "Thiago Veloso" < [hidden email]>
> >> Cc: [hidden email]
> >> 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
> < [hidden email]>
> >> 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.433714e03 and etc) produced after agregation is
> the
> >> problem...
> >> >
> >> > Thanks in advance and best wishes,
> >> >
> >> > Thiago.
> >> >
> >> >  On Fri, 3/6/11, Robert Hijmans < [hidden email]>
> >> wrote:
> >> >
> >> >> From: Robert Hijmans < [hidden email]>
> >> >> Subject: Re: [RsigGeo] Methodology to
> compare
> >> crop maps
> >> >> To: [hidden email]
> >> >> 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
> >> >> pvalues 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(xy)
> >> >> plot(xy)
> >> >> 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://rsiggeo.2731867.n2.nabble.com/Methodologytocomparecropmapstp6431598p6435902.html> >> >> Sent from the Rsiggeo mailing list
> archive at
> >> >> Nabble.com.
> >> >>
> >> >>
> _______________________________________________
> >> >> RsigGeo mailing list
> >> >> [hidden email]
> >> >> https://stat.ethz.ch/mailman/listinfo/rsiggeo> >> >>
> >> >
> >> >
> >>
> >
> >
> >
>
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo

