Subscribe to R-sig-geo feed
This is an archive for R-sig-geo, not a forum. It is not possible to post through Nabble - you may not start a new thread nor follow up an existing thread. If you wish to post, but are not subscribed to the list through the list homepage, subscribe first (through the list homepage, not via Nabble) and post from your subscribed email address. Until 2015-06-20, subscribers could post through Nabble, but policy has been changed as too many non-subscribers misunderstood the interface.
Updated: 2 hours 39 min ago

Re: Reproject MODIS data using R (results in NAs or no spatial extent)

Thu, 11/29/2018 - 13:38
Fwiw there shouldn't be any need to convert from hdf to tif - could you
please try this?

x <- readGDAL( ".../GLASS02B05.V04.A1990161.2018062.hdf")

If that works it at least removes some steps from your process.

Cheers, Mike.

On Fri, Nov 30, 2018, 05:03 Elizabeth Webb <[hidden email]> wrote:

> I am using GLASS albedo data stored here<
> ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/AVHRR/> for pre-2000 (AVHRR) data
> and here<ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/MODIS/0.05D/> for
> post-2000 data (MODIS). My end goal is to create a raster stack of each
> month that contains white sky albedo data from 1982-2015. The problem I
> have run into is that the MODIS and AVHRR data are in different spatial
> reference systems and I can't seem to reproject them to be in the same
> system.
>
> I convert from hdf to tif using R like this:
>
> fileavhrr <- ".../GLASS02B05.V04.A1990161.2018062.hdf"
> filemodis<-".../GLASS02B06.V04.A2013169.2017128.hdf"
> gdal_translate(get_subdatasets(filemodis)[10], dst_dataset =
>         ".../modis.tif")
> gdal_translate(get_subdatasets(fileavhrr)[8], projwin = c(-180,90,180,50),
> dst_dataset = ".../avhrr.tif") #ideally I'd only like data north of 50
> degrees
>
> avhrr<- raster(".../avhrr.tif")
>
> #class       : RasterLayer
> #dimensions  : 800, 7200, 5760000  (nrow, ncol, ncell)
> #resolution  : 0.05, 0.05  (x, y)
> #extent      : -180, 180, 50, 90  (xmin, xmax, ymin, ymax)
> #coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
> #values      : -32768, 32767  (min, max)
>
> modis<- raster(".../modis.tif")
>
> #class       : RasterLayer
> #dimensions  : 3600, 7200, 25920000  (nrow, ncol, ncell)
> #resolution  : 154.4376, 308.8751  (x, y)
> #extent   : -20015109, -18903159, 8895604, 10007555  (xmin, xmax, ymin,
> ymax)
> #coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
>     +b=6371007.181 +units=m +no_defs
> #values      : -32768, 32767  (min, max)
>
> Here are things I have tried:
>
> 1.) Use the MODIS Reprojection Tool<
> https://lpdaac.usgs.gov/tools/modis_reprojection_tool>. For whatever
> reason, this tool seems to think the subdatasets of the MODIS .hdf files
> are only one tile (the upper left most tile, tile 0,0) and not the global
> dataset. My understanding is that the MODIS data are global (not in
> tiles?), so I do not know why the MRT is doing this.
>
>
> 2.) Use the raster package in R.
>
> projectedMODIS <- projectRaster(modis,avhrr,method="bilinear")
>
> This returns a raster with values that are all NA:
>
> class       : RasterLayer
> dimensions  : 800, 7200, 5760000  (nrow,> ncol, ncell)
> resolution  : 0.05, 0.05  (x, y)
> extent      : -180, 180,> 50, 90  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
> values      : NA, NA  (min, max)
>
> 3.) Use the gdalUtils package in R:
>
> gdalwarp(srcfile=get_subdatasets(filemodis)[10], dstfile=
> ".../gdalMODIS_avhrr.tif", s_srs = crs(modis), t_srs =crs(avhrr) )
>
> This returns a raster with essentially no spatial extent.
>
> gdalMODISavhrr<-raster(".../gdalMODIS_avhrr.tif")
> #class       : RasterLayer
> #dimensions  : 357, 12850, 4587450  (nrow, ncol, ncell)
> #resolution  : 0.02801551, 0.02801573  (x, y)
> #extent      : -180, 179.9993, 79.99838, 90  (xmin, xmax, ymin, ymax)
> #coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
> #values      : -32768, 32767  (min, max)
>
> Any ideas on why reprojecting this MODIS data is so difficult?
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Reproject MODIS data using R (results in NAs or no spatial extent)

Thu, 11/29/2018 - 11:53
I am using GLASS albedo data stored here<ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/AVHRR/> for pre-2000 (AVHRR) data and here<ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/MODIS/0.05D/> for post-2000 data (MODIS). My end goal is to create a raster stack of each month that contains white sky albedo data from 1982-2015. The problem I have run into is that the MODIS and AVHRR data are in different spatial reference systems and I can't seem to reproject them to be in the same system.

I convert from hdf to tif using R like this:

fileavhrr <- ".../GLASS02B05.V04.A1990161.2018062.hdf"
filemodis<-".../GLASS02B06.V04.A2013169.2017128.hdf"
gdal_translate(get_subdatasets(filemodis)[10], dst_dataset =
        ".../modis.tif")
gdal_translate(get_subdatasets(fileavhrr)[8], projwin = c(-180,90,180,50), dst_dataset = ".../avhrr.tif") #ideally I'd only like data north of 50 degrees

avhrr<- raster(".../avhrr.tif")

#class       : RasterLayer
#dimensions  : 800, 7200, 5760000  (nrow, ncol, ncell)
#resolution  : 0.05, 0.05  (x, y)
#extent      : -180, 180, 50, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
#values      : -32768, 32767  (min, max)

modis<- raster(".../modis.tif")

#class       : RasterLayer
#dimensions  : 3600, 7200, 25920000  (nrow, ncol, ncell)
#resolution  : 154.4376, 308.8751  (x, y)
#extent   : -20015109, -18903159, 8895604, 10007555  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
    +b=6371007.181 +units=m +no_defs
#values      : -32768, 32767  (min, max)

Here are things I have tried:

1.) Use the MODIS Reprojection Tool<https://lpdaac.usgs.gov/tools/modis_reprojection_tool>. For whatever reason, this tool seems to think the subdatasets of the MODIS .hdf files are only one tile (the upper left most tile, tile 0,0) and not the global dataset. My understanding is that the MODIS data are global (not in tiles?), so I do not know why the MRT is doing this.


2.) Use the raster package in R.

projectedMODIS <- projectRaster(modis,avhrr,method="bilinear")

This returns a raster with values that are all NA:

class       : RasterLayer
dimensions  : 800, 7200, 5760000  (nrow,> ncol, ncell)
resolution  : 0.05, 0.05  (x, y)
extent      : -180, 180,> 50, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
values      : NA, NA  (min, max)

3.) Use the gdalUtils package in R:

gdalwarp(srcfile=get_subdatasets(filemodis)[10], dstfile= ".../gdalMODIS_avhrr.tif", s_srs = crs(modis), t_srs =crs(avhrr) )

This returns a raster with essentially no spatial extent.

gdalMODISavhrr<-raster(".../gdalMODIS_avhrr.tif")
#class       : RasterLayer
#dimensions  : 357, 12850, 4587450  (nrow, ncol, ncell)
#resolution  : 0.02801551, 0.02801573  (x, y)
#extent      : -180, 179.9993, 79.99838, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
#values      : -32768, 32767  (min, max)

Any ideas on why reprojecting this MODIS data is so difficult?


        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

How to manipulatechange the space between the numbers of axes (equal interval)? How to add a legend to a map (plot raster)?

Thu, 11/29/2018 - 10:42
Hello,
 I´m trying to plot a map in R. Is it possible to change the values of
x.axis and y.axis and have the same interval (space between the numbers)?
 How to add a legend for the map?

https://www.dropbox.com/s/60tavy6vez33mw3/Plot_map.tiff?dl=0

> class(proj22PA$P_EMmeanByROC_mergedAlgo_mergedRun_mergedData)
[1] "RasterLayer"
attr(,"package")
[1] "raster"
>
> proj22PA$P_EMmeanByROC_mergedAlgo_mergedRun_mergedData
class       : RasterLayer
band        : 2  (of  2  bands)
dimensions  : 277, 412, 114124  (nrow, ncol, ncell)
resolution  : 100, 100  (x, y)
extent      : 461214, 502414, 4272892, 4300592  (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : grid
names       : P_EMmeanByROC_mergedAlgo_mergedRun_mergedData
values      : 42, 924  (min, max)

>

Code:

plot(proj22PA$P_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
axes=FALSE,legend=FALSE,box=FALSE,par(mar=c(3,3,0.5,0.5),mgp=c(2,1,0)))
axis(side = 1, at =c(470000,480000,500000),labels= c("470000", "480000",
"500000"))
axis(side = 2, at =c(4272892,4300592),labels= c("4272892", "4300592"))
box()

Regards,

Silva

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: How to change axis (plot map)?

Sat, 11/24/2018 - 11:50
Hi Hugo,
thank you for your help and advices.

Regards,

Silva

Hugo Costa <[hidden email]> escreveu no dia sexta, 23/11/2018 à(s)
20:27:

> Probably you're looking for arguments xaxp and yaxp, used like this:
> plot(1:100,xaxp=c(0,100,4), yaxp=c(0,100,2))
>
> Please look for more info in the help. Also, if you intend to make maps,
> possibly you should spend some time getting familiar with packages, such as
> ggplot2 or tmap.
>
> Cheers
> Hugo
>
>
> Lara Silva <[hidden email]> escreveu no dia sexta, 23/11/2018
> à(s) 14:31:
>
>> Hello,
>>
>>  I am trying to change the scale/number of x-axis and y-axis (map plot -
>> Figure 1).
>>
>> https://www.dropbox.com/sh/34nddj3405q3a6p/AACe-_23Wx7bZRETYIOYjsiQa?dl=0
>> For example:
>>
>> x-axis: 465000,480000, 500000 (min. med. max values)
>> The same idea for y-axis ....
>>
>> Code:
>>
>> plot(proj22PA$OB_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
>> main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis)
>>
>>
>> plot(proj22PA$PittOB_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
>> main= "", xlab ="x.coords", ylab="y.coords",
>> cex.axis=caxis,xlim=c(465000,500000,480000), ylim=c(4270000,4300000))
>>
>> Do not change ... The map has shifted. Probably, I need a new function or
>> modified the argumente for link the information.
>>
>>
>> Another hypothesis
>>
>> plot(proj22PA$Pittosporum_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
>> main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis,
>> xaxt="n",yaxt="n")
>>
>> axis(1, seq(450000,480000,500000)) - x.axis
>> axis(2, seq(4270000,425000,4300000)) - y.axis
>>
>> Output - Figure 2. I can not add the "numbers".
>>
>> https://www.dropbox.com/sh/34nddj3405q3a6p/AACe-_23Wx7bZRETYIOYjsiQa?dl=0
>>
>> Is it possible change the scale of legend (200, 400, 600, 800)?
>>
>> Any suggestion for my doubts/problems?
>>
>> Thanks,
>>
>> Silva
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: How to sum area of adjacent patches?

Fri, 11/23/2018 - 17:48
Hi Ben,
thank you for your help. Clump is a possibility, however I was looking for
something different. In my real data, 9 rasters have to be merged (1
central raster and the 8 surrounding rasters), and running clump after
merging 9 rasters is slow and memory consuming. And this has to be repeated
thousands of times (hopefully in parallel with the rasters in memory). I
didn't explain these details to make the message short.
I'll continue searching for more alternatives, and all ideas are welcome.
Thanks
Hugo

Ben Tupper <[hidden email]> escreveu no dia sexta, 23/11/2018 à(s)
23:58:

> Hi,
>
> I think you may need to use raster::clump() and table() to label your
> regions.  Check out an old message about clump...
>
> https://stat.ethz.ch/pipermail/r-sig-geo/2017-January/025346.html
>
> It's sort of cumbersome, but I think the following works...
>
>
> library(raster)
> r1<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=0, ymx=1)
> r2<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=1, ymx=2)
> r1[]<-c(rep(24,24), rep(0,76))
> r2[]<-c(rep(0,86), rep(14,14))
> par(mfrow=c(3,2)); plot(r1,main="r1"); plot(r2,main="r2")
>
>
> r3<-merge(r1,r2) ; plot(r3, main="merged rasters")
>
> r3 <- r3
> r4[r4>0]<-38     ; plot(r4, main="desired result")
>
> r5 <- clump(r3) ; plot(r5, main = 'clumped')
> r5[is.na(r5)] <- 0
>
> r5t <- table(getValues(r5))
> for (i in names(r5t)[-1]){
>         id <- as.numeric(i)
>         print(id)
>         r5[r5 == id] <- r5t[i]
> }
>
> plot(r5, main = 'relabeled with area')
>
>
> Cheers,
> Ben
>
> P.S. Years ago I used a language called IDL.  It's histogram function has
> an output keyword called reverse_indices (more here
> https://www.harrisgeospatial.com/docs/HISTOGRAM.html#H_835179117_677188
> and here http://www.idlcoyote.com/tips/histogram_tutorial.html )  - super
> useful!  Your situation just begs for that same solution in R.
>
>
>
>
>
> > On Nov 23, 2018, at 5:13 PM, Hugo Costa <[hidden email]> wrote:
> >
> > Deal all,
> >
> > I'm merging several rasters representing patches of forest, and the pixel
> > values represent the area of the patches. The patches' area should be
> > summed along the edges of the merged rasters. I provide a reproducible
> > example. Any help is welcome.
> > Hugo
> >
> > library(raster)
> > r1<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=0, ymx=1)
> > r2<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=1, ymx=2)
> > r1[]<-c(rep(24,24), rep(0,76))
> > r2[]<-c(rep(0,86), rep(14,14))
> > par(mfrow=c(2,2)); plot(r1,main="r1"); plot(r2,main="r2")
> >
> >
> > r3<-merge(r1,r2) ; plot(r3, main="merged rasters")
> >
> > # do something
> >
> > r3[r3>0]<-38     ; plot(r3, main="desired result")
> >
> >       [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> Ben Tupper
> Bigelow Laboratory for Ocean Sciences
> 60 Bigelow Drive, P.O. Box 380
> East Boothbay, Maine 04544
> http://www.bigelow.org
>
> Ecological Forecasting: https://eco.bigelow.org/
>
>
>
>
>
>
        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: How to sum area of adjacent patches?

Fri, 11/23/2018 - 16:58
Hi,

I think you may need to use raster::clump() and table() to label your regions.  Check out an old message about clump...

https://stat.ethz.ch/pipermail/r-sig-geo/2017-January/025346.html

It's sort of cumbersome, but I think the following works...


library(raster)
r1<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=0, ymx=1)
r2<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=1, ymx=2)
r1[]<-c(rep(24,24), rep(0,76))
r2[]<-c(rep(0,86), rep(14,14))
par(mfrow=c(3,2)); plot(r1,main="r1"); plot(r2,main="r2")


r3<-merge(r1,r2) ; plot(r3, main="merged rasters")

r3 <- r3
r4[r4>0]<-38     ; plot(r4, main="desired result")

r5 <- clump(r3) ; plot(r5, main = 'clumped')
r5[is.na(r5)] <- 0

r5t <- table(getValues(r5))
for (i in names(r5t)[-1]){
        id <- as.numeric(i)
        print(id)
        r5[r5 == id] <- r5t[i]
}

plot(r5, main = 'relabeled with area')


Cheers,
Ben

P.S. Years ago I used a language called IDL.  It's histogram function has an output keyword called reverse_indices (more here https://www.harrisgeospatial.com/docs/HISTOGRAM.html#H_835179117_677188 and here http://www.idlcoyote.com/tips/histogram_tutorial.html )  - super useful!  Your situation just begs for that same solution in R.





> On Nov 23, 2018, at 5:13 PM, Hugo Costa <[hidden email]> wrote:
>
> Deal all,
>
> I'm merging several rasters representing patches of forest, and the pixel
> values represent the area of the patches. The patches' area should be
> summed along the edges of the merged rasters. I provide a reproducible
> example. Any help is welcome.
> Hugo
>
> library(raster)
> r1<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=0, ymx=1)
> r2<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=1, ymx=2)
> r1[]<-c(rep(24,24), rep(0,76))
> r2[]<-c(rep(0,86), rep(14,14))
> par(mfrow=c(2,2)); plot(r1,main="r1"); plot(r2,main="r2")
>
>
> r3<-merge(r1,r2) ; plot(r3, main="merged rasters")
>
> # do something
>
> r3[r3>0]<-38     ; plot(r3, main="desired result")
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org

Ecological Forecasting: https://eco.bigelow.org/

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

How to sum area of adjacent patches?

Fri, 11/23/2018 - 16:13
Deal all,

I'm merging several rasters representing patches of forest, and the pixel
values represent the area of the patches. The patches' area should be
summed along the edges of the merged rasters. I provide a reproducible
example. Any help is welcome.
Hugo

library(raster)
r1<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=0, ymx=1)
r2<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=1, ymx=2)
r1[]<-c(rep(24,24), rep(0,76))
r2[]<-c(rep(0,86), rep(14,14))
par(mfrow=c(2,2)); plot(r1,main="r1"); plot(r2,main="r2")


r3<-merge(r1,r2) ; plot(r3, main="merged rasters")

# do something

r3[r3>0]<-38     ; plot(r3, main="desired result")

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: How to change axis (plot map)?

Fri, 11/23/2018 - 15:27
Probably you're looking for arguments xaxp and yaxp, used like this:
plot(1:100,xaxp=c(0,100,4), yaxp=c(0,100,2))

Please look for more info in the help. Also, if you intend to make maps,
possibly you should spend some time getting familiar with packages, such as
ggplot2 or tmap.

Cheers
Hugo


Lara Silva <[hidden email]> escreveu no dia sexta, 23/11/2018
à(s) 14:31:

> Hello,
>
>  I am trying to change the scale/number of x-axis and y-axis (map plot -
> Figure 1).
>
> https://www.dropbox.com/sh/34nddj3405q3a6p/AACe-_23Wx7bZRETYIOYjsiQa?dl=0
> For example:
>
> x-axis: 465000,480000, 500000 (min. med. max values)
> The same idea for y-axis ....
>
> Code:
>
> plot(proj22PA$OB_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
> main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis)
>
>
> plot(proj22PA$PittOB_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
> main= "", xlab ="x.coords", ylab="y.coords",
> cex.axis=caxis,xlim=c(465000,500000,480000), ylim=c(4270000,4300000))
>
> Do not change ... The map has shifted. Probably, I need a new function or
> modified the argumente for link the information.
>
>
> Another hypothesis
>
> plot(proj22PA$Pittosporum_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
> main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis,
> xaxt="n",yaxt="n")
>
> axis(1, seq(450000,480000,500000)) - x.axis
> axis(2, seq(4270000,425000,4300000)) - y.axis
>
> Output - Figure 2. I can not add the "numbers".
>
> https://www.dropbox.com/sh/34nddj3405q3a6p/AACe-_23Wx7bZRETYIOYjsiQa?dl=0
>
> Is it possible change the scale of legend (200, 400, 600, 800)?
>
> Any suggestion for my doubts/problems?
>
> Thanks,
>
> Silva
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

How to change axis (plot map)?

Fri, 11/23/2018 - 07:31
Hello,

 I am trying to change the scale/number of x-axis and y-axis (map plot -
Figure 1).

https://www.dropbox.com/sh/34nddj3405q3a6p/AACe-_23Wx7bZRETYIOYjsiQa?dl=0
For example:

x-axis: 465000,480000, 500000 (min. med. max values)
The same idea for y-axis ....

Code:

plot(proj22PA$OB_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis)


plot(proj22PA$PittOB_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
main= "", xlab ="x.coords", ylab="y.coords",
cex.axis=caxis,xlim=c(465000,500000,480000), ylim=c(4270000,4300000))

Do not change ... The map has shifted. Probably, I need a new function or
modified the argumente for link the information.


Another hypothesis

plot(proj22PA$Pittosporum_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis,
xaxt="n",yaxt="n")

axis(1, seq(450000,480000,500000)) - x.axis
axis(2, seq(4270000,425000,4300000)) - y.axis

Output - Figure 2. I can not add the "numbers".

https://www.dropbox.com/sh/34nddj3405q3a6p/AACe-_23Wx7bZRETYIOYjsiQa?dl=0

Is it possible change the scale of legend (200, 400, 600, 800)?

Any suggestion for my doubts/problems?

Thanks,

Silva

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Problem with initGRASS

Fri, 11/23/2018 - 02:52
Jaime,

if the location of your GRASS installation is not known you may also
give the link2GI package a try.

try link2GI::linkGRASS7() for a "brute force" searching and linking.
This should setup your enviroment for using the newest GRASS version.

cheers Chris

On 23/11/2018 09:44, Roger Bivand wrote:
> On Fri, 23 Nov 2018, Jaime Burbano Girón wrote:
>
>> Hi everyone,
>>
>> I'm trying to use GRASS in R, but I can't get to initialize it. I'm
>> using
>> Ubuntu 18.04 and GRASS 7.4.2. This is the system information:
>>
>> Versión de GRASS: 7.4.2
>>
>> Revisión SVN de GRASS: exported
>>
>> Fecha de compilación: 2018-10-28
>>
>> Construir plataforma: x86_64-pc-linux-gnu
>>
>> GDAL: 2.2.3
>>
>> PROJ.4: 4.9.3
>>
>> GEOS: 3.6.2
>>
>> SQLite: 3.22.0
>>
>> Python: 2.7.15rc1
>>
>> wxPython: 3.0.2.0
>>
>> Plataforma: Linux-4.15.0-38-generic-x86_64-with-Ubuntu-18.04-bionic
>>
>>
>> I'm running: > initGRASS(gisBase="/usr/lib/grass74/bin",override=T)
>> And this is the error:
>
> No, ?initGRASS says:
>
> gisBase: The directory path to GRASS binaries and libraries
>
> and grass74/bin only contains (for me):
>
> $ ls bin
> grass74
>
> which is a Python script, ASCII text executable. Trying
>
> grass74/grass-7.4.2
>
> (no idea where your installer puts this):
>
> $ ls grass-7.4.2
> AUTHORS           contributors_extra.csv  fonts REQUIREMENTS.html
> bin               COPYING                 GPL.TXT  scripts
> CHANGES           demolocation            gui      share
> CITING            docs                    include  tools
> config.status     driver                  INSTALL  translators.csv
> contributors.csv  etc                     lib
>
> which contains the binaries in bin and the libraries in lib. The
> required string is what is set in the shell variable GISBASE when
> running GRASS interactively:
>
> echo $GISBASE
>
> hence the name of the argument. initGRASS() cannot discover this as
> GRASS is not running, and has not set its environment variables. I
> have added a test looking for a directory called bin in the directory
> given as the gisBase argument, check on
> https://r-forge.r-project.org/R/?group_id=2020 to see when rgrass7
> 0.1-13 (Rev: 68) has build status current, try that with your current
> setting, the correct setting and report back.
>
> Hope this clarifies,
>
> Roger
>
>
>>
>> sh: 1: g.gisenv: not foundsh: 1: g.gisenv: not foundsh: 1: g.gisenv:
>> not foundsh: 1: g.gisenv: not foundsh: 1: g.gisenv: not foundsh: 1:
>> g.version: not foundError in system(paste("g.version", get("addEXE",
>> envir = .GRASS_CACHE),  :
>>  error in running commandAdemás: Warning messages:1: In
>> system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
>> error in running command2: In system(paste(paste("g.gisenv",
>> get("addEXE", envir = .GRASS_CACHE),  :  error in running command3: In
>> system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
>> error in running command4: In system(paste(paste("g.gisenv",
>> get("addEXE", envir = .GRASS_CACHE),  :  error in running command5: In
>> system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
>> error in running command
>>
>>
>> I verified that g.gisenv and g.version are located in
>> /usr/lib/grass74/bin.
>> Any suggestion?
>>
>> Thanks in advance.
>>
>> Best,
>>
>> Jaime.
>>
>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Dr Christoph Reudenbach, Philipps-University of Marburg, Faculty of Geography, GIS and Environmental Modeling, Deutschhausstr. 10, D-35032 Marburg, fon: ++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web: gis-ma.org, giswerk.org, moc.environmentalinformatics-marburg.de


        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Problem with initGRASS

Fri, 11/23/2018 - 02:44
On Fri, 23 Nov 2018, Jaime Burbano Girón wrote:

> Hi everyone,
>
> I'm trying to use GRASS in R, but I can't get to initialize it. I'm using
> Ubuntu 18.04 and GRASS 7.4.2. This is the system information:
>
> Versión de GRASS: 7.4.2
>
> Revisión SVN de GRASS: exported
>
> Fecha de compilación: 2018-10-28
>
> Construir plataforma: x86_64-pc-linux-gnu
>
> GDAL: 2.2.3
>
> PROJ.4: 4.9.3
>
> GEOS: 3.6.2
>
> SQLite: 3.22.0
>
> Python: 2.7.15rc1
>
> wxPython: 3.0.2.0
>
> Plataforma: Linux-4.15.0-38-generic-x86_64-with-Ubuntu-18.04-bionic
>
>
> I'm running: > initGRASS(gisBase="/usr/lib/grass74/bin",override=T)
> And this is the error: No, ?initGRASS says:

gisBase: The directory path to GRASS binaries and libraries

and grass74/bin only contains (for me):

$ ls bin
grass74

which is a Python script, ASCII text executable. Trying

grass74/grass-7.4.2

(no idea where your installer puts this):

$ ls grass-7.4.2
AUTHORS           contributors_extra.csv  fonts    REQUIREMENTS.html
bin               COPYING                 GPL.TXT  scripts
CHANGES           demolocation            gui      share
CITING            docs                    include  tools
config.status     driver                  INSTALL  translators.csv
contributors.csv  etc                     lib

which contains the binaries in bin and the libraries in lib. The required
string is what is set in the shell variable GISBASE when running GRASS
interactively:

echo $GISBASE

hence the name of the argument. initGRASS() cannot discover this as GRASS
is not running, and has not set its environment variables. I have added a
test looking for a directory called bin in the directory given as the
gisBase argument, check on https://r-forge.r-project.org/R/?group_id=2020 
to see when rgrass7 0.1-13 (Rev: 68) has build status current, try that
with your current setting, the correct setting and report back.

Hope this clarifies,

Roger


>
> sh: 1: g.gisenv: not foundsh: 1: g.gisenv: not foundsh: 1: g.gisenv:
> not foundsh: 1: g.gisenv: not foundsh: 1: g.gisenv: not foundsh: 1:
> g.version: not foundError in system(paste("g.version", get("addEXE",
> envir = .GRASS_CACHE),  :
>  error in running commandAdemás: Warning messages:1: In
> system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
> error in running command2: In system(paste(paste("g.gisenv",
> get("addEXE", envir = .GRASS_CACHE),  :  error in running command3: In
> system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
> error in running command4: In system(paste(paste("g.gisenv",
> get("addEXE", envir = .GRASS_CACHE),  :  error in running command5: In
> system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
> error in running command
>
>
> I verified that g.gisenv and g.version are located in /usr/lib/grass74/bin.
> Any suggestion?
>
> Thanks in advance.
>
> Best,
>
> Jaime.
> --
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Problem with initGRASS

Thu, 11/22/2018 - 21:42
Hi everyone,

I'm trying to use GRASS in R, but I can't get to initialize it. I'm using
Ubuntu 18.04 and GRASS 7.4.2. This is the system information:

Versión de GRASS: 7.4.2

Revisión SVN de GRASS: exported

Fecha de compilación: 2018-10-28

Construir plataforma: x86_64-pc-linux-gnu

GDAL: 2.2.3

PROJ.4: 4.9.3

GEOS: 3.6.2

SQLite: 3.22.0

Python: 2.7.15rc1

wxPython: 3.0.2.0

Plataforma: Linux-4.15.0-38-generic-x86_64-with-Ubuntu-18.04-bionic


I'm running: > initGRASS(gisBase="/usr/lib/grass74/bin",override=T)
And this is the error:

sh: 1: g.gisenv: not foundsh: 1: g.gisenv: not foundsh: 1: g.gisenv:
not foundsh: 1: g.gisenv: not foundsh: 1: g.gisenv: not foundsh: 1:
g.version: not foundError in system(paste("g.version", get("addEXE",
envir = .GRASS_CACHE),  :
  error in running commandAdemás: Warning messages:1: In
system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
 error in running command2: In system(paste(paste("g.gisenv",
get("addEXE", envir = .GRASS_CACHE),  :  error in running command3: In
system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
 error in running command4: In system(paste(paste("g.gisenv",
get("addEXE", envir = .GRASS_CACHE),  :  error in running command5: In
system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
 error in running command


I verified that g.gisenv and g.version are located in /usr/lib/grass74/bin.
Any suggestion?

Thanks in advance.

Best,

Jaime.
--
Jaime Burbano Girón
Estudiante Doctoral
Doctorado en Estudios Ambientales y Rurales
Pontificia Universidad Javeriana
Bogotá, Colombia
(57-1) 3208320 Ext. 4844

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Doubt: add a boderline map

Wed, 11/21/2018 - 10:42
Yes, it worked.

Thank you very much!

Lara

Hugo Costa <[hidden email]> escreveu no dia terça, 20/11/2018 à(s)
18:15:

> Hi Lara,
>
> have you tried to run the following?
>
> map <-plot(proj$A_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
> main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis)
>
> plot(nc.sids, asp=1, main="", add=TRUE)
>
> Hugo
>
> Lara Silva <[hidden email]> escreveu no dia terça, 20/11/2018
> à(s) 16:04:
>
>> Hello,
>>
>> I got a map through the plot function.
>>
>> Code:
>>
>> map <-plot(proj$A_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
>> main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis)
>>
>> However I needed to add a border line in my study area to become more
>> perceptible (Figure 1).
>>
>> I tried the border(), rect () functions  but it gave error.
>>
>> So, I tried to create a border line
>>
>> Code:
>> nc.sids <- readShapePoly("PoligonoSM.shp")
>>
>> plot(nc.sids, asp=1, main="")
>>
>> nc.border <- unionSpatialPolygons(nc.sids, rep(1, nrow(nc.sids)))
>> plot(nc.border, asp=1, main="")
>>
>>
>> But I am having a lot of difficulties in building the script.
>>
>> How can cross the information of the 2 plots?
>>
>> https://www.dropbox.com/sh/y5by1sxjcvdl6rh/AAAIF4srkChW2FxeEfFG4d_ta?dl=0
>> Thanks.
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

job position on geostatistics at the London School of Hygiene & Tropical Medicine

Wed, 11/21/2018 - 06:57
Dear community,

We are looking for a geostatistician interested in spatial modelling of infectious diseases to join our group (London Applied and Spatial Epidemiology Research group) at the London School of Hygiene & Tropical Medicine. More details in the link below. Also, further information in our research projects in these links: www.thiswormyworld.org<http://www.thiswormyworld.org/> / http://laser.lshtm.ac.uk/

__________________________________________________________________________

Assistant Professor in Spatial Statistics and Epidemiology
Dept of Disease Control, LSHTM
Salary: �45,878 to �52,520 per annum, inclusive.
Closing Date:  Wednesday 05 December 2018

The successful applicant will be conducting research on the epidemiology and control of neglected tropical diseases (NTDs) in sub-Saharan Africa. This research is a component of the Expanded Special Project for the Elimination of NTDs within the Africa Region of the WHO. The post-holder will be responsible for leading the development and implementation of appropriate statistical methods to quantify the spatial distribution of NTDs.
Job Description and application: https://jobs.lshtm.ac.uk/vacancy.aspx?ref=ITD-DCD-2018-24


Thanks,

________________________________________________________
Jorge Cano | Assistant Professor | Department of Disease Control
London Applied & Spatial Epidemiology Research Group (LASER)
www.thiswormyworld.org<http://www.thiswormyworld.org/> / http://laser.lshtm.ac.uk/
T: +44 (0) 20 7927 2584



        [[alternative HTML version deleted]]


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Doubt: add a boderline map

Tue, 11/20/2018 - 13:15
Hi Lara,

have you tried to run the following?

map <-plot(proj$A_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis)

plot(nc.sids, asp=1, main="", add=TRUE)

Hugo

Lara Silva <[hidden email]> escreveu no dia terça, 20/11/2018
à(s) 16:04:

> Hello,
>
> I got a map through the plot function.
>
> Code:
>
> map <-plot(proj$A_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
> main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis)
>
> However I needed to add a border line in my study area to become more
> perceptible (Figure 1).
>
> I tried the border(), rect () functions  but it gave error.
>
> So, I tried to create a border line
>
> Code:
> nc.sids <- readShapePoly("PoligonoSM.shp")
>
> plot(nc.sids, asp=1, main="")
>
> nc.border <- unionSpatialPolygons(nc.sids, rep(1, nrow(nc.sids)))
> plot(nc.border, asp=1, main="")
>
>
> But I am having a lot of difficulties in building the script.
>
> How can cross the information of the 2 plots?
>
> https://www.dropbox.com/sh/y5by1sxjcvdl6rh/AAAIF4srkChW2FxeEfFG4d_ta?dl=0
> Thanks.
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Unable to use factors in raster::predict

Tue, 11/20/2018 - 13:11
This is not of much help, but I run the code, and the line
pc <- predict(logo, m, OOB=TRUE, factors=f)
did work fine. I got a Rasterlayer, which is the R logo with pixels ranging
from 0 to 1.
Hugo

Gonzalez-Mirelis, Genoveva via R-sig-Geo <[hidden email]> escreveu
no dia terça, 20/11/2018 à(s) 10:23:

> Dear list,
>
> I would like to use the 'predict' function in the 'raster' package in an
> implementation of species distribution modelling with a couple of factor
> variables. My case can be set up exactly as the cforest example listed in
> the help file. Unfortunately, I cannot get the example to work: it throws
> an error because of the unused argument 'factors', which I need. Here is
> the code I am referring to:
>
> # create a RasterStack or RasterBrick with with a set of predictor layers
> logo <- brick(system.file("external/rlogo.grd", package="raster"))
> names(logo)
>
> # known presence and absence points
> p <- matrix(c(48, 48, 48, 53, 50, 46, 54, 70, 84, 85, 74, 84, 95, 85,
>               66, 42, 26, 4, 19, 17, 7, 14, 26, 29, 39, 45, 51, 56, 46,
> 38, 31,
>               22, 34, 60, 70, 73, 63, 46, 43, 28), ncol=2)
>
> a <- matrix(c(22, 33, 64, 85, 92, 94, 59, 27, 30, 64, 60, 33, 31, 9,
>               99, 67, 15, 5, 4, 30, 8, 37, 42, 27, 19, 69, 60, 73, 3, 5,
> 21,
>               37, 52, 70, 74, 9, 13, 4, 17, 47), ncol=2)
>
> # extract values for points
> xy <- rbind(cbind(1, p), cbind(0, a))
> v <- data.frame(cbind(pa=xy[,1], extract(logo, xy[,2:3])))
>
> # cforest (other Random Forest implementation) example with factors
> argument
>
> v$red <- as.factor(round(v$red/100))
> logo$red <- round(logo[[1]]/100)
>
> library(party)
> m <- cforest(pa~., control=cforest_unbiased(mtry=3), data=v)
> f <- list(levels(v$red))
> names(f) <- 'red'
> pc <- predict( m,logo, OOB=TRUE, factors=f)
>
> Error in RET@prediction_weights(newdata = newdata, mincriterion =
> mincriterion,  :
>   unused argument (factors = f)
>
>
> Note that I needed to change the order in the arguments in the last line.
> If I run it the way it was in the example, like this:
>
> pc <- predict(logo, m, OOB=TRUE, factors=f)
>
>
> then the error is the following:
>
>
> Error in v[cells, ] <- predv :
>
>   number of items to replace is not a multiple of replacement length
>
> Also note that if I run the line without the 'factors' argument the
> resulting value (r) is a matrix, and not a raster object.
>
> I am using Package 'raster' version 2.6-7 on Windows 7 (where this worked
> fine in the past).
>
> Many thanks for any help,
>
> Genoveva
>
> Genoveva Gonzalez Mirelis, Scientist
> Institute of Marine Research
> Nordnesgaten 50
> 5005 Bergen, Norway
> Phone number +47 55238510
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Doubt: add a boderline map

Tue, 11/20/2018 - 09:04
Hello,

I got a map through the plot function.

Code:

map <-plot(proj$A_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis)

However I needed to add a border line in my study area to become more
perceptible (Figure 1).

I tried the border(), rect () functions  but it gave error.

So, I tried to create a border line

Code:
nc.sids <- readShapePoly("PoligonoSM.shp")

plot(nc.sids, asp=1, main="")

nc.border <- unionSpatialPolygons(nc.sids, rep(1, nrow(nc.sids)))
plot(nc.border, asp=1, main="")


But I am having a lot of difficulties in building the script.

How can cross the information of the 2 plots?

https://www.dropbox.com/sh/y5by1sxjcvdl6rh/AAAIF4srkChW2FxeEfFG4d_ta?dl=0
Thanks.

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Unable to use factors in raster::predict

Tue, 11/20/2018 - 03:22
Dear list,

I would like to use the 'predict' function in the 'raster' package in an implementation of species distribution modelling with a couple of factor variables. My case can be set up exactly as the cforest example listed in the help file. Unfortunately, I cannot get the example to work: it throws an error because of the unused argument 'factors', which I need. Here is the code I am referring to:

# create a RasterStack or RasterBrick with with a set of predictor layers
logo <- brick(system.file("external/rlogo.grd", package="raster"))
names(logo)

# known presence and absence points
p <- matrix(c(48, 48, 48, 53, 50, 46, 54, 70, 84, 85, 74, 84, 95, 85,
              66, 42, 26, 4, 19, 17, 7, 14, 26, 29, 39, 45, 51, 56, 46, 38, 31,
              22, 34, 60, 70, 73, 63, 46, 43, 28), ncol=2)

a <- matrix(c(22, 33, 64, 85, 92, 94, 59, 27, 30, 64, 60, 33, 31, 9,
              99, 67, 15, 5, 4, 30, 8, 37, 42, 27, 19, 69, 60, 73, 3, 5, 21,
              37, 52, 70, 74, 9, 13, 4, 17, 47), ncol=2)

# extract values for points
xy <- rbind(cbind(1, p), cbind(0, a))
v <- data.frame(cbind(pa=xy[,1], extract(logo, xy[,2:3])))

# cforest (other Random Forest implementation) example with factors argument

v$red <- as.factor(round(v$red/100))
logo$red <- round(logo[[1]]/100)

library(party)
m <- cforest(pa~., control=cforest_unbiased(mtry=3), data=v)
f <- list(levels(v$red))
names(f) <- 'red'
pc <- predict( m,logo, OOB=TRUE, factors=f)

Error in RET@prediction_weights(newdata = newdata, mincriterion = mincriterion,  :
  unused argument (factors = f)


Note that I needed to change the order in the arguments in the last line. If I run it the way it was in the example, like this:

pc <- predict(logo, m, OOB=TRUE, factors=f)


then the error is the following:


Error in v[cells, ] <- predv :

  number of items to replace is not a multiple of replacement length

Also note that if I run the line without the 'factors' argument the resulting value (r) is a matrix, and not a raster object.

I am using Package 'raster' version 2.6-7 on Windows 7 (where this worked fine in the past).

Many thanks for any help,

Genoveva

Genoveva Gonzalez Mirelis, Scientist
Institute of Marine Research
Nordnesgaten 50
5005 Bergen, Norway
Phone number +47 55238510


        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Dealing with non-nested geographic units using the RQGIS and tigris R packages

Mon, 11/19/2018 - 09:34
I recently came across this tutorial suggesting that the RQGIS and tigris R
packages could be used to overcome issues with non-nested geographies:

https://rpubs.com/corey_sparks/256942

I'm trying to apply this to some data I have where I would like to look at
Census Block Groups nested within Zip Code Tabulation Areas. Specifically,
I'd like to compute a segregation index using Census Block Groups as the
lower units and Zip Code Tabulation Areas as the higher level units.

I am confused about how these R packages deal with this issue of non-nested
units. Here's the explanation from the tutorial:

"In order to get around this problem [non-nested units], we can use the
QGIS application to perform a geographic intersection between two spatial
polygon layers that we want to use. This will allow us to merge the
population data at the census tract level to a higher level, so we can
calculate a segregation index."

Can anyone say more about how this is possible? I can run the code in the
tutorial but I don't understand what it's doing.

Thanks!

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

set up months to a time serie

Mon, 11/19/2018 - 05:26
Dear all,

I am trying to set a sequence of  48000 months to a simulation between
10200 - 6201  before present.

However I always get the following message
Error in charToDate(TS) : #character string is not in a standard
unambiguous format
Could anyone help me?
####
library(zoo)
TS1<-brick("TS1.nc")
# dimensions  : 12, 12, 144, 48000  (nrow, ncol, ncell, nlayers)
TS1 <- as.yearmon(10200 - seq(0, 48000)/12)

####

Moreover, I would like to detect changes inside my rasterbrick serie, I
read that package bfast can detect changes in such of time series and
files, anyone experienced with bfast could help me with it as well?

thank you all,

Best

Jackson

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Pages