Re: [R-sig-Geo] projectRaster across the antimeridian

2016-07-01 Thread Bacou, Melanie 

Tony,
I can confirm this seems like an issue with `projectRaster()`. Using 
`gdalwarp` instead on the same tile seems to work.


# Download geoTIFF
BS.terra <- 
curl::curl_download("https://lance.modaps.eosdis.nasa.gov/imagery/subsets/?subset=BeringSea.2016182.terra.721.2km.tif;,

  "BeringSea.2016182.terra.721.2km.tif")
BS.terra <- brick("BeringSea.2016182.terra.721.2km.tif")

# Reproject using `gdalwarp`
system("gdalwarp -overwrite -s_srs EPSG:4326 -t_srs EPSG:3035 -multi -of 
GTiff BeringSea.20e16182.terra.721.2km.tif 
BeringSea.2016182.terra.721.2km.azi.tif")


# Load reprojected raster
BS.terra.sa1 <- brick("BeringSea.2016182.terra.721.2km.azi.tif")

# Reproject with projectRaster()
BS.terra.sa2 <- projectRaster(from=BS.terra, crs=CRS("+init=epsg:3035"))

# Plot
par(mfrow=c(2,1))
plotRGB(BS.terra)
plotRGB(BS.terra.sa1)

# This one errors out
plotRGB(BS.terra.sa2)
> Error in grDevices::rgb(RGB[, 1], RGB[, 2], RGB[, 3], alpha = alpha, 
max = scale) :

  color intensity -0.00130533, not in [0,1]

See https://dl.dropboxusercontent.com/u/30925475/Capture.PNG

--Mel.

On 7/1/2016 2:48 PM, Fischbach, Anthony wrote:

The solution to the previous post addressing projection of rasters across
the antimeridian
(
http://r-sig-geo.2731867.n2.nabble.com/coercing-marmap-bathy-object-to-raster-spanning-the-antimeridian-td7589828.html#a7589841
)
no-longer functions
(as of 24 June 2016).

Below I offer permutations of the solution offered in the previous post,
to no avail.

Is this a bug in the raster package, or may my code be adjusted to produce
a projected raster that includes both the western and eastern hemisphere of
the original image?

-Tony

require(sp)
require(raster)
setwd('D:/Temp')
## Define projections
prj.geo<-CRS('+proj=longlat +datum=WGS84')
prj.StudyArea <- CRS("+proj=aeqd +lat_0=70 +lon_0=-170")
prj.StudyAreaPositive <- CRS("+proj=aeqd +lat_0=70 +lon_0=190")
## Download geoJPEG
Yesterday <- Sys.Date() -1
cat(URL<- format(Yesterday, '
http://lance-modis.eosdis.nasa.gov/subsets/index.php?subset=BeringSea.%Y%j.terra.721.2km.jpg'),
fill=T)
cat(DestFile<-format(Yesterday,'BeringSea.%Y%j.terra.721.2km.jpg' ,
sep=""), fill=T)##
try(ds<-download.file(url=URL, destfile=DestFile, method="libcurl", quiet =
F, mode = "wb", cacheOK = F))
## Cast as raster stack
BS.terra<-stack(DestFile)
## Assign projection
projection(BS.terra)<- prj.geo
ymax(BS.terra) <- 70
ymin(BS.terra) <- 58
xmax(BS.terra) = 360-155
xmin(BS.terra) = 360-190
## project the raster stack
(BS.terra.sa1<-projectRaster(from=BS.terra,  res=5000,
crs=prj.StudyAreaPositive, method='ngb', over=TRUE))
(BS.terra.sa2<-projectRaster(from=BS.terra,  res=5000,
crs=prj.StudyAreaPositive, method='ngb', over=F))
(BS.terra.sa3<-projectRaster(from=BS.terra,  res=5000, crs=prj.StudyArea,
method='ngb', over=TRUE))
(BS.terra.sa4<-projectRaster(from=BS.terra,  res=5000, crs=prj.StudyArea,
method='ngb', over=F))
## plot
par(mfrow=c(5,1))
plotRGB(BS.terra)
plotRGB(BS.terra.sa1)
plotRGB(BS.terra.sa2)
plotRGB(BS.terra.sa3)
plotRGB(BS.terra.sa4)
#

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo


___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] determining values for nmax, nmin and maxdist

2016-07-01 Thread Tim Peterson 
Hi Daniel,

I've never had much success with rules of thumb (eg variogram range). I 
suggest you look into more objective approaches - eg see paper below. 
When I've played with global methods similar to Abedini I've found that 
cross-validation errors (with and without replacement) show a search 
distance that gives a minimum error. Re local  search distance 
estimation (see Abedini), I suspect that for any moderate to large 
problem the computations becomes prohibitive.

M.J. Abedini, M. Nasseri, D.H. Burn, The use of a genetic 
algorithm-based search strategy in geostatistics: application to a set 
of anisotropic piezometric head data, Computers & Geosciences, Volume 
41, April 2012, Pages 136-146, ISSN 0098-3004, 
http://dx.doi.org/10.1016/j.cageo.2011.08.024

Cheers

Tim

On 07/02/2016 07:47 AM, Dan Turenne wrote:
> Hello listers,
>
>
> As part of the research I am doing I'm attempting to fit some kriging models, 
> and I was wondering if there are any "rule of thumb" guidelines for selecting 
> the values for the nmax, nmin and maxdist parameters?  Intuitively I would 
> set maxdist equal to the range of the variogram but I'm not sure about the 
> other parameters.  Any advice on this subject would be appreciated.
>
>
> Thanks,
>
>
> Daniel Turenne
>
> University of Manitoba
>
>   [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo


-- 

Dr. Tim Peterson

The Department of Infrastructure Engineering
The University of Melbourne, 3010 Australia
T: +61 3 8344 9950 , M: +61 0438 385 937 


Dept. profile : 
http://www.ie.unimelb.edu.au/people/staff.php?person_ID=141135
Research Gate : https://www.researchgate.net/profile/Tim_Peterson7
Google Scholar: 
http://scholar.google.com.au/citations?user=kkYJLF4J=en=ao

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] projectRaster across the antimeridian

2016-07-01 Thread Fischbach, Anthony 
The solution to the previous post addressing projection of rasters across
the antimeridian
(
http://r-sig-geo.2731867.n2.nabble.com/coercing-marmap-bathy-object-to-raster-spanning-the-antimeridian-td7589828.html#a7589841
)
no-longer functions
(as of 24 June 2016).

Below I offer permutations of the solution offered in the previous post,
to no avail.

Is this a bug in the raster package, or may my code be adjusted to produce
a projected raster that includes both the western and eastern hemisphere of
the original image?

-Tony

require(sp)
require(raster)
setwd('D:/Temp')
## Define projections
prj.geo<-CRS('+proj=longlat +datum=WGS84')
prj.StudyArea <- CRS("+proj=aeqd +lat_0=70 +lon_0=-170")
prj.StudyAreaPositive <- CRS("+proj=aeqd +lat_0=70 +lon_0=190")
## Download geoJPEG
Yesterday <- Sys.Date() -1
cat(URL<- format(Yesterday, '
http://lance-modis.eosdis.nasa.gov/subsets/index.php?subset=BeringSea.%Y%j.terra.721.2km.jpg'),
fill=T)
cat(DestFile<-format(Yesterday,'BeringSea.%Y%j.terra.721.2km.jpg' ,
sep=""), fill=T)##
try(ds<-download.file(url=URL, destfile=DestFile, method="libcurl", quiet =
F, mode = "wb", cacheOK = F))
## Cast as raster stack
BS.terra<-stack(DestFile)
## Assign projection
projection(BS.terra)<- prj.geo
ymax(BS.terra) <- 70
ymin(BS.terra) <- 58
xmax(BS.terra) = 360-155
xmin(BS.terra) = 360-190
## project the raster stack
(BS.terra.sa1<-projectRaster(from=BS.terra,  res=5000,
crs=prj.StudyAreaPositive, method='ngb', over=TRUE))
(BS.terra.sa2<-projectRaster(from=BS.terra,  res=5000,
crs=prj.StudyAreaPositive, method='ngb', over=F))
(BS.terra.sa3<-projectRaster(from=BS.terra,  res=5000, crs=prj.StudyArea,
method='ngb', over=TRUE))
(BS.terra.sa4<-projectRaster(from=BS.terra,  res=5000, crs=prj.StudyArea,
method='ngb', over=F))
## plot
par(mfrow=c(5,1))
plotRGB(BS.terra)
plotRGB(BS.terra.sa1)
plotRGB(BS.terra.sa2)
plotRGB(BS.terra.sa3)
plotRGB(BS.terra.sa4)
#

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] spatial anova

2016-07-01 Thread Dexter Locke 
Normally I would agree with Steve, but this is a topic I'm interested in
and have searched extensively. I have not yet found a useable
implementation. If there is some R code that can do this, I would love to
hear about it.

In April Babak Naimi  posted on this list a very similar
question, but I did not see any replies.

For the underlying maths, see these papers:

Griffith, D. A. (1978). A spatially adjusted ANOVA model. *Geographical
Analysis*, *X*(3), 298–301. doi:10.1016/0166-0462(92)90034-X

Griffith, D. A. (1992). A spatially adjusted N-way ANOVA model. *Regional
Science and Urban Economics*, *22*(3), 347–369.
doi:10.1016/0166-0462(92)90034-X

These are also related:

Long, J. a., Nelson, T. a., & Wulder, M. a. (2010). Local indicators for
categorical data: Impacts of scaling decisions. *Canadian Geographer*, *54*(1),
15–28. doi:10./j.1541-0064.2009.00300.x

Pietrzak, M. B., Wilk, J., Kossowski, T., & Bivand, R. (2014). The
application of local indicators for categorical data (LICD) in the spatial
analysis of economic development, (14). Retrieved from
http://ideas.repec.org/p/pes/wpaper/2014no14.html

The old SpaceStat used to have a spatial anova option (white paper by
Anselin):

http://siteresources.worldbank.org/INTPGI/Resources/342674-1092157888460/Anselin.spacestatTutorial.pdf

SpaceStat 4 appears to be able to perform a spatial anova:

https://www.biomedware.com/files/SpaceStat_4.0_Documentation.pdf


Just to be sure I hadn't missed anything I double checked and found this:

library(spdep); help(joincount.multi) and joincount.test(). So it looks
like Roger Bivand provided the needed code once again! Thanks so much.

HTH,
Dexter




On Tue, Jun 28, 2016 at 12:20 PM, Steve Friedman 
wrote:

> I would start by reading google and R posts regarding spatial anova. There
> is plenty of help among that literature.  Formulate your model, share it,
> and then you might have better questions.
>
> Steve
>
> On Tuesday, June 28, 2016, Csany22 csany  wrote:
>
> > Hi R-sig-geo people,
> >
> >
> > what would be the best way to run spatial anova in R? I have a dependent
> > (continuous), independent (categorical) variables and lat and long fo
> reach
> > of the cell(pixel) based observations. I need to run anova. as would like
> > to compare the means (Tukey) of the independent variable.
> >
> > [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org 
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
>
> --
> Steve Friedman
> Cell: 517-648-6290
> Web:  www.envisioningecosystemsdecisions.com
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo