Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

2016-02-26 Thread Jonathan Greenberg
Folks:

I went ahead and pushed a new version of gdalUtils to r-forge:

install.packages("gdalUtils", repos="http://R-Forge.R-project.org;)

gdal_translate has a new option to attempt to fix the issue above.  Can you
confirm if e.g.:

gdal_translate(...,config='GDAL_NETCDF_BOTTOMUP="YES"')
works for your application?

--j

On Fri, Feb 26, 2016 at 9:17 AM Jonathan Greenberg 
wrote:

> Folks:
>
> Just noticed this thread -- I see I didn't include a "--config" option
> with any of the gdalUtils functions (it isn't one of the documented
> parameters on the individual utility website, but it seems it would have
> allowed you to run the GDAL_NETCDF_BOTTOMUP without setting a system
> environment variable -- see
> http://lists.osgeo.org/pipermail/gdal-dev/2014-July/039452.html).  If I
> make this tweak and pushed it to r-forge, would one of you be willing to
> see if it solves the problem?
>
> I assume this will be something that would be needed for any gdal utility
> that allows an "of" to be set?
>
> --j
>
> On Wed, Feb 24, 2016 at 2:37 AM Chris Reudenbach <
> reudenb...@uni-marburg.de> wrote:
>
>> Agus,
>>
>> sorry for the addon but I think I have to provide a correction of the
>> corner coordinates (e.g. the extent values):
>>
>> In the example that you have posted below I did calculate the extent
>> using the domain center coordinate and the WRF grid resolution in meter
>> and the number of rows and cols.
>>
>> Since Dominik provides a link to the file description of the netcdf file
>> I think it is more accurate to reproject the  corner coordinates  as
>> given by the netcdf header variables (NC_GLOBAL#corner_lats,
>> NC_GLOBAL#corner_lons). Assuming that your variable "dctmx" (which I can
>> not identfy it in the nc file)  is of type "Mass" (staggered = M) the
>> correct corner coordinates are stored as the first 4 entrys of the dump
>> snippet below:
>>
>> lats<-"
>>
>> NC_GLOBAL#corner_lats={12.355667,50.26619,50.26619,12.355667,12.308136,50.18787,50.18787,12.308136,12.210785,50.403816,50.403816,12.210785,12.163345,50.325382,50.325382,12.163345}"
>>
>> lons<-"
>>
>> NC_GLOBAL#corner_lons={-131.43678,-151.29639,-48.703613,-68.563232,-131.5851,-151.51157,-48.488434,-68.414917,-131.38828,-151.41891,-48.581085,-68.611725,-131.53641,-151.63441,-48.36557,-68.463593}"
>>
>> after cleaning and converting the strings you may calculate the corner
>> coordinates:
>>
>>
>> library(proj4)
>> ## project mass corner coordinates to lambertian
>> llMass <- ptransform(cbind(clon[1],clat[1])/180*pi,'+proj=longlat
>> +datum=WGS84 +no_defs',proj)
>> ulMass <- ptransform(cbind(clon[2],clat[2])/180*pi,'+proj=longlat
>> +datum=WGS84 +no_defs',proj)
>> lrMass <- ptransform(cbind(clon[3],clat[3])/180*pi,'+proj=longlat
>> +datum=WGS84 +no_defs',proj)
>> urMass <- ptransform(cbind(clon[4],clat[4])/180*pi,'+proj=longlat
>> +datum=WGS84 +no_defs',proj)
>> wrfLccExtMass<-extent(ulMass[1],lrMass[1],llMass[2],ulMass[2])
>>
>>
>> According to this the correct extent for mass variables should be:
>>
>> extent(wrfLccExtMass)
>> class   : Extent
>> xmin: -3575343
>> xmax: 3575342
>> ymin: -2293058
>> ymax: 2306330
>>
>>
>> hope this is correct now
>>
>> cheers Chris
>>
>> Am 24.02.2016 um 05:12 schrieb Agus Camacho:
>> > With the help of everybody i finally got this to work, so here is a
>> script
>> > that does the job of reprojecting both, a raster layer obtained from a
>> .nc
>> > and some locations in order to overplot them using plot.raster or
>> mapview.
>> > Im using a combination of the advices of Dominik, Michael and Chris.
>> >
>> >
>> > require(ncdf4)
>> > require(raster)
>> >
>> >
>> > setwd("C:/~")
>> >
>> >
>> > r=raster("geo_em.d01.nc",
>> >varname="dctmx")# days of ctmax events
>> >
>> > # Set extent and projections of rasters for plotting
>> > # chris gave me the orig data fom the nc file because i could not
>> install
>> > gdal
>> > xmin= -3545999
>> > xmax= 3546000
>> > ymin = -2286000
>> > ymax=2286000
>> > pr<- "+proj=lcc +lat_1=25 +lat_2=45 +lat_0=38.08 +lon_0=-100 +x_0=0
>> > +y_0=0"
>> >
>> > wrfLccExt<-extent(xmin,xmax,ymin,ymax)
>> >
>> > extent(r) <- extent(wrfLccExt)
>> > projection(r) <- pr
>> >
>> > # get and prepare  urosaurus locations
>> >
>> > x=read.csv("C:/Users/Agus/Dropbox/Functional Biogeography -
>> > Copy/scripts/class 4 maxent model/clean urosaurus records.csv")
>> > x=x[,1:3]
>> > colnames(x)
>> > coordinates(x)=~lon+lat
>> > proj4string(x)<- CRS("+proj=longlat +datum=WGS84")
>> > x=spTransform(x, pr)
>> >
>> >
>> > # get prepare and plot wrld_simpl
>> > require(maptools)
>> > require(sp)
>> > data(wrld_simpl)
>> > namer <- spTransform(subset(wrld_simpl, NAME %in% c("United States",
>> > "Canada", "Mexico")), prj)
>> >
>> > #plot with raster
>> > plot(r, cex.main=.7,legend=F)
>> > points(x)
>> > plot(namer, add = TRUE)
>> >
>> > # plot with mapview (cool!)
>> > m=mapview(r)
>> > u=mapview(x)
>> > m+u
>> >
>> >
>> > Thanks 

Re: [R-sig-Geo] Error running Mann-Kendall trend test on a raster stack

2016-02-26 Thread Tim Appelhans

Dear Thiago,
have a look at the gimms package, available on CRAN. The function of 
interest is significantTau().


https://github.com/environmentalinformatics-marburg/gimms/blob/master/R/significantTau.R

Hope this helps,
Tim

On 27.02.2016 09:49, Thiago V. dos Santos wrote:

Dear colleagues,

I have a raster stack with 89 layers, each layer representing yearly 
precipitation.

I am trying to use the function rkt (from package "rkt") to detect a possible 
trend in my precipitation time series.

Its usage is pretty simple, and this is how it runs on a data frame:


library(rkt)

# generate some random data
myDF <- data.frame(Date = seq.Date(as.Date("2011-01-01"), 
as.Date("2099-01-01"), by='year'),
value = runif(89,0,1200))

# extract year from date, as numeric
myDF$year <- as.numeric(format(myDF$Date, "%Y"))

# run mann-kendall test
kenn <- rkt(myDF$year, myDF$value)
kenn



Now, this is my attempt to apply this test on a raster stack, using a calc 
function:


library(raster)

# Create the date sequence
idx <- myDF$Date

# Create raster stack and apply the date
r <- raster(ncol=50, nrow=50)
s <- stack(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r), 0, 
1200
s <- setZ(s, idx)
s

# Now I define a function
tsfun <- function(x) {
 year <- as.numeric(substr(getZ(x), 1, 4))
 r.kenn <- rkt(year, x)
 return(r.kenn)
}

# I apply the function
raster.kenn <- calc(s, fun=tsfun)




And the error message I get is:

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
cannot use this function

What am I doing wrong here?
Greetings,
  -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota

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


--
#
Tim Appelhans
Department of Geography
Environmental Informatics
Philipps Universität Marburg
Deutschhausstraße 12
Raum 00A08
35032 Marburg (Paketpost: 35037 Marburg)
Germany

Tel +49 (0) 6421 28-25957

http://environmentalinformatics-marburg.de/

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


[R-sig-Geo] Error running Mann-Kendall trend test on a raster stack

2016-02-26 Thread Thiago V. dos Santos
Dear colleagues,

I have a raster stack with 89 layers, each layer representing yearly 
precipitation.

I am trying to use the function rkt (from package "rkt") to detect a possible 
trend in my precipitation time series.

Its usage is pretty simple, and this is how it runs on a data frame:


library(rkt)

# generate some random data
myDF <- data.frame(Date = seq.Date(as.Date("2011-01-01"), 
as.Date("2099-01-01"), by='year'),
value = runif(89,0,1200))

# extract year from date, as numeric
myDF$year <- as.numeric(format(myDF$Date, "%Y"))

# run mann-kendall test
kenn <- rkt(myDF$year, myDF$value)
kenn



Now, this is my attempt to apply this test on a raster stack, using a calc 
function:


library(raster)

# Create the date sequence
idx <- myDF$Date

# Create raster stack and apply the date
r <- raster(ncol=50, nrow=50)
s <- stack(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r), 0, 
1200
s <- setZ(s, idx)
s

# Now I define a function
tsfun <- function(x) {
year <- as.numeric(substr(getZ(x), 1, 4))
r.kenn <- rkt(year, x)
return(r.kenn)
}

# I apply the function
raster.kenn <- calc(s, fun=tsfun)




And the error message I get is:

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
cannot use this function

What am I doing wrong here?
Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota

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


Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

2016-02-26 Thread Jonathan Greenberg
Folks:

Just noticed this thread -- I see I didn't include a "--config" option with
any of the gdalUtils functions (it isn't one of the documented parameters
on the individual utility website, but it seems it would have allowed you
to run the GDAL_NETCDF_BOTTOMUP without setting a system environment
variable -- see
http://lists.osgeo.org/pipermail/gdal-dev/2014-July/039452.html).  If I
make this tweak and pushed it to r-forge, would one of you be willing to
see if it solves the problem?

I assume this will be something that would be needed for any gdal utility
that allows an "of" to be set?

--j

On Wed, Feb 24, 2016 at 2:37 AM Chris Reudenbach 
wrote:

> Agus,
>
> sorry for the addon but I think I have to provide a correction of the
> corner coordinates (e.g. the extent values):
>
> In the example that you have posted below I did calculate the extent
> using the domain center coordinate and the WRF grid resolution in meter
> and the number of rows and cols.
>
> Since Dominik provides a link to the file description of the netcdf file
> I think it is more accurate to reproject the  corner coordinates  as
> given by the netcdf header variables (NC_GLOBAL#corner_lats,
> NC_GLOBAL#corner_lons). Assuming that your variable "dctmx" (which I can
> not identfy it in the nc file)  is of type "Mass" (staggered = M) the
> correct corner coordinates are stored as the first 4 entrys of the dump
> snippet below:
>
> lats<-"
>
> NC_GLOBAL#corner_lats={12.355667,50.26619,50.26619,12.355667,12.308136,50.18787,50.18787,12.308136,12.210785,50.403816,50.403816,12.210785,12.163345,50.325382,50.325382,12.163345}"
>
> lons<-"
>
> NC_GLOBAL#corner_lons={-131.43678,-151.29639,-48.703613,-68.563232,-131.5851,-151.51157,-48.488434,-68.414917,-131.38828,-151.41891,-48.581085,-68.611725,-131.53641,-151.63441,-48.36557,-68.463593}"
>
> after cleaning and converting the strings you may calculate the corner
> coordinates:
>
>
> library(proj4)
> ## project mass corner coordinates to lambertian
> llMass <- ptransform(cbind(clon[1],clat[1])/180*pi,'+proj=longlat
> +datum=WGS84 +no_defs',proj)
> ulMass <- ptransform(cbind(clon[2],clat[2])/180*pi,'+proj=longlat
> +datum=WGS84 +no_defs',proj)
> lrMass <- ptransform(cbind(clon[3],clat[3])/180*pi,'+proj=longlat
> +datum=WGS84 +no_defs',proj)
> urMass <- ptransform(cbind(clon[4],clat[4])/180*pi,'+proj=longlat
> +datum=WGS84 +no_defs',proj)
> wrfLccExtMass<-extent(ulMass[1],lrMass[1],llMass[2],ulMass[2])
>
>
> According to this the correct extent for mass variables should be:
>
> extent(wrfLccExtMass)
> class   : Extent
> xmin: -3575343
> xmax: 3575342
> ymin: -2293058
> ymax: 2306330
>
>
> hope this is correct now
>
> cheers Chris
>
> Am 24.02.2016 um 05:12 schrieb Agus Camacho:
> > With the help of everybody i finally got this to work, so here is a
> script
> > that does the job of reprojecting both, a raster layer obtained from a
> .nc
> > and some locations in order to overplot them using plot.raster or
> mapview.
> > Im using a combination of the advices of Dominik, Michael and Chris.
> >
> >
> > require(ncdf4)
> > require(raster)
> >
> >
> > setwd("C:/~")
> >
> >
> > r=raster("geo_em.d01.nc",
> >varname="dctmx")# days of ctmax events
> >
> > # Set extent and projections of rasters for plotting
> > # chris gave me the orig data fom the nc file because i could not install
> > gdal
> > xmin= -3545999
> > xmax= 3546000
> > ymin = -2286000
> > ymax=2286000
> > pr<- "+proj=lcc +lat_1=25 +lat_2=45 +lat_0=38.08 +lon_0=-100 +x_0=0
> > +y_0=0"
> >
> > wrfLccExt<-extent(xmin,xmax,ymin,ymax)
> >
> > extent(r) <- extent(wrfLccExt)
> > projection(r) <- pr
> >
> > # get and prepare  urosaurus locations
> >
> > x=read.csv("C:/Users/Agus/Dropbox/Functional Biogeography -
> > Copy/scripts/class 4 maxent model/clean urosaurus records.csv")
> > x=x[,1:3]
> > colnames(x)
> > coordinates(x)=~lon+lat
> > proj4string(x)<- CRS("+proj=longlat +datum=WGS84")
> > x=spTransform(x, pr)
> >
> >
> > # get prepare and plot wrld_simpl
> > require(maptools)
> > require(sp)
> > data(wrld_simpl)
> > namer <- spTransform(subset(wrld_simpl, NAME %in% c("United States",
> > "Canada", "Mexico")), prj)
> >
> > #plot with raster
> > plot(r, cex.main=.7,legend=F)
> > points(x)
> > plot(namer, add = TRUE)
> >
> > # plot with mapview (cool!)
> > m=mapview(r)
> > u=mapview(x)
> > m+u
> >
> >
> > Thanks to all!
> > Agus
> >
> > 2016-02-23 13:11 GMT-07:00 Agus Camacho :
> >
> >> Thanks for that Dominik,
> >>
> >> Giving that projection to either the locations, the raster layer
> generated
> >> from the .nc file, or both, still did not work. I keep having locations
> >> that should be on land falling far on the sea. Might this be a problem
> >> derived from using raster with a file whose original grid distances are
> not
> >> constant?
> >>
> >> Here is a link with the original file which has the original coordinate
> >> data.
> >>
> 

Re: [R-sig-Geo] resolution of openmap() raster layers

2016-02-26 Thread Chris Reudenbach

Agus,

Mapview is using leaflet as engine. Due to this you will have the 
control icons on the map because first of all it is designed for 
interactive mapping within RStudio/R.


I think there are two different approaches to save your maps:

If you want to have a dump of the mapviewobject (but including the 
graphical buttons) you'll find a descriptat stackoverflow

http://stackoverflow.com/questions/31336898/how-to-save-leaflet-in-rstudio-map-as-png-or-jpg-file


You may also use  the spplot function from mapview which is designed for 
basic static mapping and to make usable the adavantages of spplot.


Note even if you are dealing with the mapview map object the spplot 
function  uses the Openstreetmap package for retrieving  the background 
maps  (e.g. 
http://www.inside-r.org/packages/cran/OpenStreetMap/docs/openmap). You 
can use the spplot syntax for designing your maps. Up to now this static 
plotting function is still pretty  basic but you may have a try:
spplot(mapView(nica["POP2000"]),colorkey=FALSE,  lwd= 15, alpha.regions 
= 0.9,

map.type="stamen-watercolor" )

I think for using the spplot it is better to install the current stable 
from github:

library(devtools)
install_github("environmentalinformatics-marburg/mapview", ref = "master")


cheers chris


Am 26.02.2016 um 12:27 schrieb Agustin Lobo:

Stunning!
Can I remove the buttons for saving to a bmp file?
What attribution should be used for publishing?
Agus

On Thu, Feb 25, 2016 at 7:42 PM, Chris Reudenbach
 wrote:

Hi,

if you just want to map the data, mapview could be an option that among
others avoid the pixel stretching.

require(mapview)
require(raster)
nica <- getData("GADM", country="NIC", level=0)

mapview(nica)

mapview(nica,zcol = "POP2000", color = "#FFA500", lwd= 5, alpha.regions =
0.4)


cheers Chris




Am 25.02.2016 um 18:49 schrieb Barry Rowlingson:

On Thu, Feb 25, 2016 at 5:11 PM, Agustin Lobo 
wrote:

Is there any way to download the raster layers
of openmap() with an increased resolution?
I find the quality of the labels very low,
or am I doing something wrong? i.e.

require(raster)
require(mapmisc)
nica <- getData("GADM", country="NIC", level=0)
nicabg <- openmap(nica, path="landscape")
plot(nicabg)

   Map tiles from OpenStreetMap and other map tile providers are images
designed to be shown at a fixed resolution. When you plot them in an R
graphics window you could be stretching them so that each pixel in the
original maps to 1.273 pixels on your screen. So some kind of
interpolation or nearest neighbour replacement has to be done, and
this makes text labels look bad. Other line work will look bad too.

   If you try and download more map tiles at a higher resolution then
you'll find the labels are now way too small, because what you've
downloaded are map tiles designed for a higher zoom level on a web
browser. Web map browsers have a fixed set of zoom values that
correspond to the resolution of the map tiles. With an R window, you
are free to choose odd zoom factors that give the ugly behaviour.

   If you can resize your R window exactly right then you might get
something that looks good!

   The alternative is to build a background map yourself from
OpenStreetMap *vector* data and some code and some styling. Or use a
map tile provider that doesn't have text labels and add them to
selected places with R graphics commands. Lines and polygons will
still be stretched and a bit "jaggy" but our eyes don't notice this as
much as badly scaled text.

Barry

___
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


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


Re: [R-sig-Geo] Problems with mapmisc::insetMap()

2016-02-26 Thread Patrick Brown
Hi.

1- `insetMap` returns a matrix of x and y coordinates, you needn't
plot all of them.

```
 nica <- getData("GADM", country="NIC", level=0)
 nicabg <- openmap(nica, path="landscape")
 map.new(nicabg,1,mar=c(2, 2, 2, 0) + 0.1)
 plot(nicabg,axes=TRUE)
 plot(nica,add=TRUE)
 loc <- insetMap(nica,pos=c(-89.9,8.5), width=0.3, col=NA, lty=0,pch=".")
 points(loc[1,1], loc[1,2], pch=3)
```

2- `dev.size()` gives you the size of the plotting window.  the best
way to control the plot size is to save it to a file.

```
png("myfile.png", height=400, width=300)
plot(nicabg)
dev.off()
```

p

On Fri, Feb 26, 2016 at 5:02 AM, Agustin Lobo  wrote:
> Much better now!
>
> nica <- getData("GADM", country="NIC", level=0)
> nicabg <- openmap(nica, path="landscape")
> map.new(nicabg,1,mar=c(2, 2, 2, 0) + 0.1)
> plot(nicabg,axes=TRUE)
> plot(nica,add=TRUE)
> loc <- insetMap(nica,pos=c(-89.9,8.5), width=0.3, col=NA, lty=0,pch=".")
> points(loc)
>
> Just 2 questions:
> 1. In the example you provide, the resulting mark is too thick. I've
> made a thinner version with pch=".", but
> it would nice having the option of just drawing the central point and
> not the window.
> 2. Is there a way to calculate the size of the plotting window so that
> there is no wasted white space
> between the map and the axes? I can do this by hand, but then cannot
> read the size.
>
> Thanks
>
> On Thu, Feb 25, 2016 at 6:46 PM, Patrick Brown
>  wrote:
>> Hello all, my first post here!
>>
>>
>> mapmisc (1.4.8 from r-forge) now allows
>>
>>
>> require(rgdal)
>> require(raster)
>> require(mapmisc)
>> nica <- getData("GADM", country="NIC", level=0)
>> nicabg <- openmap(nica, path="landscape")
>> map.new(nicabg, 0.7,mar=c(2, 2, 2, 0) + 0.1)
>> plot(nicabg,axes=TRUE)
>> plot(nica,add=TRUE)
>> loc = insetMap(nica,pos=c(-85,12), width=0.3, col=NA, lty=0)
>> points(loc)
>>
>> p
>>
>> ___
>> 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] odd behaviour of llgridlines() regarding axes labels

2016-02-26 Thread Michael Sumner
1) is ok with xpd clipping turned off

op <- par(xpd = NA)
plot(a,axes=FALSE)
llgridlines(a, plotLines=TRUE, plotLabels=TRUE)
par(op)

2)  same
op <- par(xpd = NA)
plot(a,axes=FALSE,mar=c(4, 3, 2, 3)+0.1)
llgridlines(a, plotLines=TRUE, plotLabels=TRUE)
par(op)

3)  Consider using graticule package instead? You can control the extents
independently of the meridian/parallel lines (and optionally provided a
target projection. )

You have to give it the longs/lats you want, there's no automatic detection
from an object or plot - but for me that's the price of control over what I
get.

install.packages("graticule")
library(graticule)
plot(a,axes=TRUE)
library(graticule)
g <- graticule(seq(-90, -82, by = 2), seq(10, 16, by = 2), xlim =
par("usr")[1:2], ylim = par("usr")[3:4])
plot(g, add = TRUE, lty = 2)

Cheers, Mike.


On Fri, 26 Feb 2016 at 22:21 Agustin Lobo  wrote:

> I see an odd behaviour of llgridlines() regarding axes labels:
>
> require(rgdal)
> require(raster)
> a <- extent(c(-90,-81.5625,8.407168,16.63619))
> a <- as(a, 'SpatialPolygons')
> projection(a) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
> +towgs84=0,0,0"
>
> No axis labels on x axis:
> plot(a,axes=FALSE)
> llgridlines(a, plotLines=TRUE, plotLabels=TRUE)
>
> #same with:
> plot(a,axes=FALSE,mar=c(4, 3, 2, 3)+0.1)
> llgridlines(a, plotLines=TRUE, plotLabels=TRUE)
>
> #and 2 vertical axes with:
> plot(a,axes=TRUE)
> llgridlines(a, plotLines=TRUE, plotLabels=TRUE)
>
> Any hint?
>
> Thanks
> Agus
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> 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
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo


Re: [R-sig-Geo] resolution of openmap() raster layers

2016-02-26 Thread Agustin Lobo
Stunning!
Can I remove the buttons for saving to a bmp file?
What attribution should be used for publishing?
Agus

On Thu, Feb 25, 2016 at 7:42 PM, Chris Reudenbach
 wrote:
> Hi,
>
> if you just want to map the data, mapview could be an option that among
> others avoid the pixel stretching.
>
> require(mapview)
> require(raster)
> nica <- getData("GADM", country="NIC", level=0)
>
> mapview(nica)
>
> mapview(nica,zcol = "POP2000", color = "#FFA500", lwd= 5, alpha.regions =
> 0.4)
>
>
> cheers Chris
>
>
>
>
> Am 25.02.2016 um 18:49 schrieb Barry Rowlingson:
>>
>> On Thu, Feb 25, 2016 at 5:11 PM, Agustin Lobo 
>> wrote:
>>>
>>> Is there any way to download the raster layers
>>> of openmap() with an increased resolution?
>>> I find the quality of the labels very low,
>>> or am I doing something wrong? i.e.
>>>
>>> require(raster)
>>> require(mapmisc)
>>> nica <- getData("GADM", country="NIC", level=0)
>>> nicabg <- openmap(nica, path="landscape")
>>> plot(nicabg)
>>
>>   Map tiles from OpenStreetMap and other map tile providers are images
>> designed to be shown at a fixed resolution. When you plot them in an R
>> graphics window you could be stretching them so that each pixel in the
>> original maps to 1.273 pixels on your screen. So some kind of
>> interpolation or nearest neighbour replacement has to be done, and
>> this makes text labels look bad. Other line work will look bad too.
>>
>>   If you try and download more map tiles at a higher resolution then
>> you'll find the labels are now way too small, because what you've
>> downloaded are map tiles designed for a higher zoom level on a web
>> browser. Web map browsers have a fixed set of zoom values that
>> correspond to the resolution of the map tiles. With an R window, you
>> are free to choose odd zoom factors that give the ugly behaviour.
>>
>>   If you can resize your R window exactly right then you might get
>> something that looks good!
>>
>>   The alternative is to build a background map yourself from
>> OpenStreetMap *vector* data and some code and some styling. Or use a
>> map tile provider that doesn't have text labels and add them to
>> selected places with R graphics commands. Lines and polygons will
>> still be stretched and a bit "jaggy" but our eyes don't notice this as
>> much as badly scaled text.
>>
>> Barry
>>
>> ___
>> 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

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


[R-sig-Geo] odd behaviour of llgridlines() regarding axes labels

2016-02-26 Thread Agustin Lobo
I see an odd behaviour of llgridlines() regarding axes labels:

require(rgdal)
require(raster)
a <- extent(c(-90,-81.5625,8.407168,16.63619))
a <- as(a, 'SpatialPolygons')
projection(a) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0"

No axis labels on x axis:
plot(a,axes=FALSE)
llgridlines(a, plotLines=TRUE, plotLabels=TRUE)

#same with:
plot(a,axes=FALSE,mar=c(4, 3, 2, 3)+0.1)
llgridlines(a, plotLines=TRUE, plotLabels=TRUE)

#and 2 vertical axes with:
plot(a,axes=TRUE)
llgridlines(a, plotLines=TRUE, plotLabels=TRUE)

Any hint?

Thanks
Agus

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


Re: [R-sig-Geo] Problems with mapmisc::insetMap()

2016-02-26 Thread Agustin Lobo
Much better now!

nica <- getData("GADM", country="NIC", level=0)
nicabg <- openmap(nica, path="landscape")
map.new(nicabg,1,mar=c(2, 2, 2, 0) + 0.1)
plot(nicabg,axes=TRUE)
plot(nica,add=TRUE)
loc <- insetMap(nica,pos=c(-89.9,8.5), width=0.3, col=NA, lty=0,pch=".")
points(loc)

Just 2 questions:
1. In the example you provide, the resulting mark is too thick. I've
made a thinner version with pch=".", but
it would nice having the option of just drawing the central point and
not the window.
2. Is there a way to calculate the size of the plotting window so that
there is no wasted white space
between the map and the axes? I can do this by hand, but then cannot
read the size.

Thanks

On Thu, Feb 25, 2016 at 6:46 PM, Patrick Brown
 wrote:
> Hello all, my first post here!
>
>
> mapmisc (1.4.8 from r-forge) now allows
>
>
> require(rgdal)
> require(raster)
> require(mapmisc)
> nica <- getData("GADM", country="NIC", level=0)
> nicabg <- openmap(nica, path="landscape")
> map.new(nicabg, 0.7,mar=c(2, 2, 2, 0) + 0.1)
> plot(nicabg,axes=TRUE)
> plot(nica,add=TRUE)
> loc = insetMap(nica,pos=c(-85,12), width=0.3, col=NA, lty=0)
> points(loc)
>
> p
>
> ___
> 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] Multiplying a raster object

2016-02-26 Thread Tadaishi Yatabe-Rodriguez
Thanks Loic,

That helps a lot, actually. I used the wrong factor for transforming.
Still adding the values of my raster don't add up to the points of the
original ppp object. So I suppose there are better approaches than
multiplying the raster by factor. If someone else wants to chip in on this
it'll be great.

Thanks again,

Tada

On Wed, Feb 24, 2016 at 4:56 AM, Loïc Dutrieux  wrote:

>
>
> On 02/23/2016 09:22 PM, Tadaishi Yatabe-Rodriguez wrote:
>
>> Hi community,
>>
>> I have a raster object of a kernel density, where the linear unit is
>> meter,
>> and when I plot it it gives me very small density values (units per sq
>> meter, I suppose). I figure that if I multiply it, say by a 1000, it will
>> give density values in units per sq km.
>>
>
> I don't know much about kernel densities, but if you want to convert
> units/m^2 to units/km^2 you'd have to multiply by 1 000 000 instead of
> 1000. Hope that helps explaining the incoherence in your results.
>
> It looks good on the legend, but I
>> wonder: is this a correct thing to do? I wonder this because when I add
>> the
>> values of the raster it does not add up to the number of points in my
>> original ppp object from which I created the density.
>>
>> Thanks,
>>
>> Tada
>>
>>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>



-- 
*Tadaishi Yatabe*
DVM, MPVM, PhD (C)
Center for Animal Disease Modeling and Surveillance (CADMS)
Department of Medicine and Epidemiology
School of Veterinary Medicine
University of California Davis

http://tadaishi.wix.com/tada

[[alternative HTML version deleted]]

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