Re: [R-sig-Geo] Multiplying a raster object
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
Re: [R-sig-Geo] Problems with mapmisc::insetMap()
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
[R-sig-Geo] odd behaviour of llgridlines() regarding axes labels
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] resolution of openmap() raster layers
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] odd behaviour of llgridlines() regarding axes labels
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] odd behaviour of llgridlines() regarding axes labels
Thanks, that fixes it. I will try graticules as well. Agus On Fri, Feb 26, 2016 at 12:45 PM, Michael Sumner wrote: > 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 ___ 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()
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] resolution of openmap() raster layers
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] adapting spatial points and wrld_smpl to a reference system implicit in a .nc file
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. > >> > >> https://www.dropbox.com/s/qpt5twtunhy3x3x/geo_em.
[R-sig-Geo] Error running Mann-Kendall trend test on a raster stack
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] Error running Mann-Kendall trend test on a raster stack
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
Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a reference system implicit in a .nc file
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 to all! >> > Agus >>