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 32 min ago

Re: How to objectively subset cities by population

Thu, 07/27/2017 - 09:27
Dear Dr. Gräler,
Thanks for your contribution. I very much enjoyed the clustering suggestion, and it seems to be available in R's leaflet through the "markerClusterOptions" command. It could solve my problem, so I will take a closer look at that.
Regarding your first suggestion, can you point me out some example that uses the overlaid grid approach?
Thanks, -- Thiago V. dos Santos
PhD studentLand and Atmospheric ScienceUniversity of Minnesota


On Thursday, July 27, 2017, 3:00:55 AM CDT, Dr. Benedikt Gräler <[hidden email]> wrote:

Dear Thiago,

if you want them spatially evenly distributed, you could overlay a grid
and select the largest per grid box - or maybe more intuitive, select
the largest per predefined administrative areas (counties/postal
codes/...). This could also change based on zoom-level. An alternative
is to group sensors and expand and zoom in by clicking on the group (see
e.g. [1]).

HTH,

  Ben

[1] http://sensorweb.demo.52north.org/client/#/map


On 27/07/2017 06:09, Thiago V. dos Santos via R-sig-Geo wrote:
> Dear all,
>
> I have temperature records of nearly 1200 locations in southern Brazil.
>
> I am writing a shiny app that will show an interactive map with the locations plotted as circles, where the user can click a location to see its temperature time series.
>
> However, if I show all the locations in the map, it will look really bad, too cramped.
>
> Therefore, in an attempt to make the map look a bit cleaner, I am trying to think of an objective way to subset the locations. My initial approach would be to show only the "largest" locations, i.e. the ones with a population above a certain threshold.
>
> The problem is: the distribution of the population is so positively skewed that I am having a hard time determining the optimal cutoff point.
>
> Does anybody here know any tool or method, possibly spatial, that can assist me with this analysis?
>
> These are the locations I am working with:
>
> #-------------------------------
> # Download and summarize
> locs <- read.csv("https://www.dropbox.com/s/ykdd8x1mlc76klt/locations.csv?raw=1")
> hist(locs$Population)
> summary(locs$Population)
>
> # Convert to spatial points and plot
> require(sp)
> coordinates(locs) <- cbind(locs$Lon , locs$Lat)
> plot(locs)
> bubble(locs,"Population")
> #-------------------------------
>
> Thanks in advance,
>  -- Thiago V. dos Santos
>
> PhD student
> Land and Atmospheric Science
> University of Minnesota
>     [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Dr. Benedikt Gräler
52°North Initiative for Geospatial Open Source Software GmbH
Martin-Luther-King-Weg 24
48155 Muenster, Germany

E-Mail: [hidden email]
Fon: +49-(0)-251/396371-39
Fax: +49-(0)-251/396371-11

http://52north.org/
Twitter: @FiveTwoN

General Managers: Dr. Albert Remke, Dr. Andreas Wytzisk
Local Court Muenster HRB 10849
_______________________________________________
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 objectively subset cities by population

Thu, 07/27/2017 - 09:11
Hi Kent,
Thank you very much for the response.
This is exactly the same approach I am using now, with a different radius formula though.
It is not at all terrible, but do you see how overlapped the largest cities are, and how difficult it is to click on some smaller cities?
This is why I am looking for an objective way to filter out the smallest cities, while at the same time keeping as much information on the map as possible (i.e. not leaving too many "blank" areas). 
Best, -- Thiago V. dos Santos
PhD studentLand and Atmospheric ScienceUniversity of Minnesota


On Thursday, July 27, 2017, 7:29:44 AM CDT, Kent Johnson <[hidden email]> wrote:


Date: Thu, 27 Jul 2017 04:09:53 +0000 (UTC)
From: "Thiago V. dos Santos" <[hidden email]>
To: R-sig-geo Mailing List <[hidden email]>
Subject: [R-sig-Geo] How to objectively subset cities by population
Message-ID: <1634627903.292923. [hidden email]>
Content-Type: text/plain; charset="UTF-8"

Dear all,

I have temperature records of nearly 1200 locations in southern Brazil.

I am writing a shiny app that will show an interactive map with the locations plotted as circles, where the user can click a location to see its temperature time series.

However, if I show all the locations in the map, it will look really bad, too cramped.


Have you considered using leaflet to make an interactive map? This is a good start:
library(leaflet)locs <- read.csv("https://www.dropbox.com/s/ykdd8x1mlc76klt/locations.csv?raw=1")
leaflet(locs) %>% addTiles() %>%   addCircles(radius=~20*sqrt(Population), label=~as.character(Geocode),              stroke=FALSE, fillOpacity=0.5)
You can configure popups to show HTML or get Shiny events on click, for example clicking on a city could display the time series in a separate panel.
Docs and many examples on the leaflet for R web site:
https://rstudio.github.io/leaflet/

Kent
        [[alternative HTML version deleted]]

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

Re: How to objectively subset cities by population

Thu, 07/27/2017 - 07:29
>
> Date: Thu, 27 Jul 2017 04:09:53 +0000 (UTC)
> From: "Thiago V. dos Santos" <[hidden email]>
> To: R-sig-geo Mailing List <[hidden email]>
> Subject: [R-sig-Geo] How to objectively subset cities by population
> Message-ID: <[hidden email]>
> Content-Type: text/plain; charset="UTF-8"
>
> Dear all,
>
> I have temperature records of nearly 1200 locations in southern Brazil.
>
> I am writing a shiny app that will show an interactive map with the
> locations plotted as circles, where the user can click a location to see
> its temperature time series.
>
> However, if I show all the locations in the map, it will look really bad,
> too cramped.
>
Have you considered using leaflet to make an interactive map? This is a
good start:

library(leaflet)
locs <- read.csv("
https://www.dropbox.com/s/ykdd8x1mlc76klt/locations.csv?raw=1")

leaflet(locs) %>% addTiles() %>%
  addCircles(radius=~20*sqrt(Population), label=~as.character(Geocode),
             stroke=FALSE, fillOpacity=0.5)

You can configure popups to show HTML or get Shiny events on click, for
example clicking on a city could display the time series in a separate
panel.

Docs and many examples on the leaflet for R web site:
https://rstudio.github.io/leaflet/

Kent

        [[alternative HTML version deleted]]

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

Re: How to objectively subset cities by population

Thu, 07/27/2017 - 03:00
Dear Thiago,

if you want them spatially evenly distributed, you could overlay a grid
and select the largest per grid box - or maybe more intuitive, select
the largest per predefined administrative areas (counties/postal
codes/...). This could also change based on zoom-level. An alternative
is to group sensors and expand and zoom in by clicking on the group (see
e.g. [1]).

HTH,

  Ben

[1] http://sensorweb.demo.52north.org/client/#/map


On 27/07/2017 06:09, Thiago V. dos Santos via R-sig-Geo wrote:
> Dear all,
>
> I have temperature records of nearly 1200 locations in southern Brazil.
>
> I am writing a shiny app that will show an interactive map with the locations plotted as circles, where the user can click a location to see its temperature time series.
>
> However, if I show all the locations in the map, it will look really bad, too cramped.
>
> Therefore, in an attempt to make the map look a bit cleaner, I am trying to think of an objective way to subset the locations. My initial approach would be to show only the "largest" locations, i.e. the ones with a population above a certain threshold.
>
> The problem is: the distribution of the population is so positively skewed that I am having a hard time determining the optimal cutoff point.
>
> Does anybody here know any tool or method, possibly spatial, that can assist me with this analysis?
>
> These are the locations I am working with:
>
> #-------------------------------
> # Download and summarize
> locs <- read.csv("https://www.dropbox.com/s/ykdd8x1mlc76klt/locations.csv?raw=1")
> hist(locs$Population)
> summary(locs$Population)
>
> # Convert to spatial points and plot
> require(sp)
> coordinates(locs) <- cbind(locs$Lon , locs$Lat)
> plot(locs)
> bubble(locs,"Population")
> #-------------------------------
>
> Thanks in advance,
>   -- Thiago V. dos Santos
>
> PhD student
> Land and Atmospheric Science
> University of Minnesota
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Dr. Benedikt Gräler
52°North Initiative for Geospatial Open Source Software GmbH
Martin-Luther-King-Weg 24
48155 Muenster, Germany

E-Mail: [hidden email]
Fon: +49-(0)-251/396371-39
Fax: +49-(0)-251/396371-11

http://52north.org/
Twitter: @FiveTwoN

General Managers: Dr. Albert Remke, Dr. Andreas Wytzisk
Local Court Muenster HRB 10849

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
b_graeler.vcf (458 bytes) Download Attachment

How to objectively subset cities by population

Wed, 07/26/2017 - 23:09
Dear all,

I have temperature records of nearly 1200 locations in southern Brazil.

I am writing a shiny app that will show an interactive map with the locations plotted as circles, where the user can click a location to see its temperature time series.

However, if I show all the locations in the map, it will look really bad, too cramped.

Therefore, in an attempt to make the map look a bit cleaner, I am trying to think of an objective way to subset the locations. My initial approach would be to show only the "largest" locations, i.e. the ones with a population above a certain threshold.

The problem is: the distribution of the population is so positively skewed that I am having a hard time determining the optimal cutoff point.

Does anybody here know any tool or method, possibly spatial, that can assist me with this analysis?

These are the locations I am working with:

#-------------------------------
# Download and summarize
locs <- read.csv("https://www.dropbox.com/s/ykdd8x1mlc76klt/locations.csv?raw=1")
hist(locs$Population)
summary(locs$Population)

# Convert to spatial points and plot
require(sp)
coordinates(locs) <- cbind(locs$Lon , locs$Lat)
plot(locs)
bubble(locs,"Population")
#-------------------------------

Thanks in advance,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota
        [[alternative HTML version deleted]]

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

Finding overlapping polygons

Wed, 07/26/2017 - 08:30
Dear everyone,

I am studying the potential spatial distribution between vertebrate
invasive species and native ones (threatened) in the souther cone of
South America. The invasive species' distribution was grossly
approximated using an ecoregional approach, which is not relevant here.
The analysis is based on shapefiles for all terrestrial mammals and
birds of the world obtained from BirdLife International
(http://www.biodiversityinfo.org/spcdownload/r5h8a1/) and IUCN
(http://www.iucnredlist.org/technical-documents/spatial-data) public
domains.

Our goal: using a two-step approach, identify native species (polygons)
whose distributionoverlap -totally or partially- with that of the invasive.

Our problem: Some species (with small distribution) which are known to
occur within the invasive distribution range are being left out of the
final selection (checkedin QGIS). I suspect the function "isin" in the
fastshp package (2nd loop) fails to select some of the overlapping
polygons in Step 2, or some other problem. If a reproducible example is
absolutely required, I can upload to dropbox.

For the sake of simplicity, I will present present a case example
including one invasive species (e.g. wild boar) and all terrestrial
mammals.

# Analysis of Wild Boar distributional overlap with IUCN mammal species
install.packages("rgdal") install.packages("sp")
install.packages("ggplot2") install.packages("rgeos")
install.packages("raster") install.packages("fastshp", repos =
"http://rforge.net", type = "source")install.packages("cleangeo")

mypath <- ("/Users/ ... /Terrestrial mammals") files <-
list.files(path=mypath, pattern = "\\.shp$", full.names=T)
length(files)  # 5,465 mammal species # Load study area

study_area <- read.shp("C:/Users/.../Study_area.shp", format="table",
close=FALSE) salon <- range(study_area$x, na.rm=T) salat <-
range(study_area$y, na.rm=T) bounding <- c(salon, salat) box <-
data.frame(salon, salat) colnames(box) <- c("Longitud","Latitude") box

# Step 1. Run all native IUCN species' distribution through a loop and
select those which overlap with the invasive's bounding box distribution

keep <- rep(NA, length(files)) length(keep)  # 5,465 for(i in
1:length(files)){   r <- read.shp(files[i], format="table")   rlon <-
range(r$x)   rlat <- range(r$y) ## Is out of the bounding box? test.lat
<- (min(rlat) > max(salat)) |     (max(rlat) < min(salat))   test.lon <-
(min(rlon) > max(salon)) |     (max(rlon) < min(salon)) keep[i] <-
!(test.lon | test.lat) } keep <- which(keep==1) # Index of shapefiles
that passed the first filter length(keep) # 773 species passed 1st
filter # Step 2. Load invasive species distribution area

study_area <- readOGR("/Users/.../Study area","Study_area") study_area
<- spTransform(study_area, CRS("+init=epsg:32721")) # To UTM # Get a
report of geometry validity & issues for a sp spatial object report <-
clgeo_CollectionReport(study_area) summary <-
clgeo_SummaryReport(report) summary # Remove holes and simplify geometry
of study area for computational efficiency outer <-
Filter(function(x){x@ringDir==1}, study_area@polygons[[1]]@Polygons) #
PREGUNTAR study_area <- SpatialPolygons(list(Polygons(outer, ID=1)))
study_area <- disaggregate(study_area) areas <- gArea(study_area,
byid=TRUE) #Remove polygons smaller than 3000 km2 study_area <-
study_area[which(areas > 3E9),] study_area <- gSimplify(study_area,
20000, topologyPreserve=TRUE) # Simplify geometry
proj4string(study_area) <- CRS("+init=epsg:32721") study_area <-
spTransform(study_area, CRS("+proj=longlat +datum=WGS84 +no_defs
+ellps=WGS84 +towgs84=0,0,0")) xy <- fortify(study_area) #
SpatialPolygons to data frame keep2 <- rep(NA,length(keep)) # 773 sp.
that passed first filter for(i in 1:length(keep)){   r <-
read.shp(files[keep[i]], format="polygon")   isin <- any(inside(r,
x=xy$long, y=xy$lat), na.rm=T)   keep2[i] <- isin } keep2 <-
which(keep2) # Index of shapefiles that passed the second sort
length(keep2)  # 532 species that passed the second filter species <-
list.files(path=mypath, pattern = "\\.shp$",full.names=F) species <-
sub(".shp", "", species) species <- species[keep[keep2]]



---
This email has been checked for viruses by Avast antivirus software.
https://www.avast.com/antivirus

        [[alternative HTML version deleted]]


.

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

Finding overlapping polygons

Wed, 07/26/2017 - 08:28
Dear everyone,

I am studying the potential spatal distribution between vertebrate
invasive species and native, threatened ones in the souther cone of
South America. The invasive species' distribution was grossly
approximated using an ecoregion approach, which is not relevant here.
Then, I downloaded shapefiles for all terrestrial mammals and birds of
the world from BirdLife International
(http://www.biodiversityinfo.org/spcdownload/r5h8a1/) and IUCN
(http://www.iucnredlist.org/technical-documents/spatial-data) public
domains.

Our goal: using a two-step approach, identify native species (polygons)
whose distributionoverlap -totally or partially- with that of the invasive.

Our problem: Some species (with small distribution) which are known to
occur within the invasive distribution range are being left out of the
final selection (checkedin QGIS). I suspect the function "isin" in the
fastshp package (2nd loop) fails to select some of the overlapping
polygons in Step 2, or some other problem. If a reproducible example is
absolutely required, I can upload to dropbox.

For the sake of simplicity, I will present present a case example
including one invasive species (e.g. wild boar) and all terrestrial
mammals.

# Analysis of Wild Boar distributional overlap with IUCN mammal species
install.packages("rgdal") install.packages("sp")
install.packages("ggplot2") install.packages("rgeos")
install.packages("raster") install.packages("fastshp", repos =
"http://rforge.net", type = "source")install.packages("cleangeo")

mypath <- ("/Users/ ... /Terrestrial mammals") files <-
list.files(path=mypath, pattern = "\\.shp$", full.names=T)
length(files)  # 5,465 mammal species # Load study area

study_area <- read.shp("C:/Users/.../Study_area.shp", format="table",
close=FALSE) salon <- range(study_area$x, na.rm=T) salat <-
range(study_area$y, na.rm=T) bounding <- c(salon, salat) box <-
data.frame(salon, salat) colnames(box) <- c("Longitud","Latitude") box

# Step 1. Run all native IUCN species' distribution through a loop and
select those which overlap with the invasive's bounding box distribution

keep <- rep(NA, length(files)) length(keep)  # 5,465 for(i in
1:length(files)){   r <- read.shp(files[i], format="table")   rlon <-
range(r$x)   rlat <- range(r$y) ## Is out of the bounding box? test.lat
<- (min(rlat) > max(salat)) |     (max(rlat) < min(salat))   test.lon <-
(min(rlon) > max(salon)) |     (max(rlon) < min(salon)) keep[i] <-
!(test.lon | test.lat) } keep <- which(keep==1) # Index of shapefiles
that passed the first filter length(keep) # 773 species passed 1st
filter # Step 2. Load invasive species distribution area

study_area <- readOGR("/Users/.../Study area","Study_area") study_area
<- spTransform(study_area, CRS("+init=epsg:32721")) # To UTM # Get a
report of geometry validity & issues for a sp spatial object report <-
clgeo_CollectionReport(study_area) summary <-
clgeo_SummaryReport(report) summary # Remove holes and simplify geometry
of study area for computational efficiency outer <-
Filter(function(x){x@ringDir==1}, study_area@polygons[[1]]@Polygons) #
PREGUNTAR study_area <- SpatialPolygons(list(Polygons(outer, ID=1)))
study_area <- disaggregate(study_area) areas <- gArea(study_area,
byid=TRUE) #Remove polygons smaller than 3000 km2 study_area <-
study_area[which(areas > 3E9),] study_area <- gSimplify(study_area,
20000, topologyPreserve=TRUE) # Simplify geometry
proj4string(study_area) <- CRS("+init=epsg:32721") study_area <-
spTransform(study_area, CRS("+proj=longlat +datum=WGS84 +no_defs
+ellps=WGS84 +towgs84=0,0,0")) xy <- fortify(study_area) #
SpatialPolygons to data frame keep2 <- rep(NA,length(keep)) # 773 sp.
that passed first filter for(i in 1:length(keep)){   r <-
read.shp(files[keep[i]], format="polygon")   isin <- any(inside(r,
x=xy$long, y=xy$lat), na.rm=T)   keep2[i] <- isin } keep2 <-
which(keep2) # Index of shapefiles that passed the second sort
length(keep2)  # 532 species that passed the second filter species <-
list.files(path=mypath, pattern = "\\.shp$",full.names=F) species <-
sub(".shp", "", species) species <- species[keep[keep2]]



---
This email has been checked for viruses by Avast antivirus software.
https://www.avast.com/antivirus

        [[alternative HTML version deleted]]

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

Re: bivariate spatial correlation in R

Wed, 07/26/2017 - 05:07
On Wed, 26 Jul 2017, Rafael Pereira wrote:

> Roger,
>
> This example was provided only for the sake or making the code easily
> reproducible for others and I'm more interested in how the bi-variate Moran
> could be implemented in R, but your comments are very much welcomed and
> I've made changes to the question.
>
> My actual case study looks at bi-variate spatial correlation between (a)
> average household income per capita and (b) proportion of jobs in the city
> that are accessible under 60 minutes by transit. I don't think I could use
> rates in this case but I will normalize the variables using
> scale(data$variable). Please provide a reproducible example, either with a link to a data
subset, or using a builtin data set. My guess is that you do not need
bi-variate spatial correlation at all, but rather a spatial regression.

The "causal" variable would then the the proportion of jobs accessible
within 60 minutes by transit, though this is extremely blunt, and lots of
other covariates (demography, etc.) impact average household income per
capita (per block/tract?). Since there are many missing variables in your
specification, any spatial correlation would be most closely associated
with them (demography, housing costs, education, etc.), and the choice of
units of measurement would dominate the outcome.

This is also why bi-variate spatial correlation is seldom a good idea, I
believe. It can be done, but most likely shouldn't, unless it can be
motivated properly.

By the way, the weighted and FDR-corrected SAD local Moran's I p-values of
the black/white ratio for Oregon (your toy example) did deliver the goods
- if you zoom in in mapview::mapview, you can see that it detects a rate
hotspot between the rivers.

Roger

>
> best,
>
> Rafael H M Pereira
>
> On Mon, Jul 24, 2017 at 7:56 PM, Roger Bivand <[hidden email]> wrote:
>
>> On Mon, 24 Jul 2017, Rafael Pereira wrote:
>>
>> Hi all,
>>>
>>> I would like to ask whether some you conducted bi-variate spatial
>>> correlation in R.
>>>
>>> I know the bi-variate Moran's I is not implemented in the spdep library.
>>> I left a question on SO but also wanted to hear if anyone if the mainlist
>>> have come across this.
>>> https://stackoverflow.com/questions/45177590/map-of-bivariat
>>> e-spatial-correlation-in-r-bivariate-lisa
>>>
>>> I also know Roger Bivand has implemented the L index proposed by Lee
>>> (2001)
>>> in spdep, but I'm not I'm not sure whether the L local correlation
>>> coefficients can be interpreted the same way as the local Moran's I
>>> coefficients. I couldn't find any reference commenting on this issue. I
>>> would very much appreciate your thoughts this.
>>>
>>
>> In the SO question, and in the follow-up, your presumably throw-away
>> example makes fundamental mistakes. The code in spdep by Virgilio
>> Gómez-Rubio is for uni- and bivariate L, and produces point values of local
>> L. This isn't the main problem, which is rather that you are not taking
>> account of the underlying population counts, nor shrinking any estimates of
>> significance to accommodate population sizes. Population sizes vary from 0
>> to 11858, with the lower quartile at 3164 and upper 5698:
>> plot(ecdf(oregon.tract$pop2000)). Should you be comparing rates in stead?
>> These are also compositional variables (sum to pop2000, or 1 if rates) with
>> the other missing components. You would probably be better served by tools
>> examining spatial segregation, such as for example the seg package.
>>
>> The 0 count populations cause problems for an unofficial alternative, the
>> black/white ratio:
>>
>> oregon.tract1 <- oregon.tract[oregon.tract$white > 0,]
>> oregon.tract1$rat <- oregon.tract1$black/oregon.tract1$white
>> nb <- poly2nb(oregon.tract1)
>> lw <- nb2listw(nb)
>>
>> which should still be adjusted by weighting:
>>
>> lm0 <- lm(rat ~ 1, weights=pop2000, data=oregon.tract1)
>>
>> I'm not advising this, but running localmoran.sad on this model output
>> yields SAD p-values < 0.05 after FDR correction only in contiguous tracts
>> on the Washington state line in Portland between the Columbia and
>> Williamette rivers. So do look at the variables you are using before
>> rushing into things.
>>
>> Hope this clarifies,
>>
>> Roger
>>
>>
>>> best,
>>>
>>> Rafael HM Pereira
>>> http://urbandemographics.blogspot.com
>>>
>>>         [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> 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.
>> voice: +47 55 95 93 55; e-mail: [hidden email]
>> Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
>> http://orcid.org/0000-0003-2392-6140
>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> 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.
voice: +47 55 95 93 55; e-mail: [hidden email]
Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
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

Re: bivariate spatial correlation in R

Wed, 07/26/2017 - 03:18
Roger,

This example was provided only for the sake or making the code easily
reproducible for others and I'm more interested in how the bi-variate Moran
could be implemented in R, but your comments are very much welcomed and
I've made changes to the question.

My actual case study looks at bi-variate spatial correlation between (a)
average household income per capita and (b) proportion of jobs in the city
that are accessible under 60 minutes by transit. I don't think I could use
rates in this case but I will normalize the variables using
scale(data$variable).

best,

Rafael H M Pereira

On Mon, Jul 24, 2017 at 7:56 PM, Roger Bivand <[hidden email]> wrote:

> On Mon, 24 Jul 2017, Rafael Pereira wrote:
>
> Hi all,
>>
>> I would like to ask whether some you conducted bi-variate spatial
>> correlation in R.
>>
>> I know the bi-variate Moran's I is not implemented in the spdep library.
>> I left a question on SO but also wanted to hear if anyone if the mainlist
>> have come across this.
>> https://stackoverflow.com/questions/45177590/map-of-bivariat
>> e-spatial-correlation-in-r-bivariate-lisa
>>
>> I also know Roger Bivand has implemented the L index proposed by Lee
>> (2001)
>> in spdep, but I'm not I'm not sure whether the L local correlation
>> coefficients can be interpreted the same way as the local Moran's I
>> coefficients. I couldn't find any reference commenting on this issue. I
>> would very much appreciate your thoughts this.
>>
>
> In the SO question, and in the follow-up, your presumably throw-away
> example makes fundamental mistakes. The code in spdep by Virgilio
> Gómez-Rubio is for uni- and bivariate L, and produces point values of local
> L. This isn't the main problem, which is rather that you are not taking
> account of the underlying population counts, nor shrinking any estimates of
> significance to accommodate population sizes. Population sizes vary from 0
> to 11858, with the lower quartile at 3164 and upper 5698:
> plot(ecdf(oregon.tract$pop2000)). Should you be comparing rates in stead?
> These are also compositional variables (sum to pop2000, or 1 if rates) with
> the other missing components. You would probably be better served by tools
> examining spatial segregation, such as for example the seg package.
>
> The 0 count populations cause problems for an unofficial alternative, the
> black/white ratio:
>
> oregon.tract1 <- oregon.tract[oregon.tract$white > 0,]
> oregon.tract1$rat <- oregon.tract1$black/oregon.tract1$white
> nb <- poly2nb(oregon.tract1)
> lw <- nb2listw(nb)
>
> which should still be adjusted by weighting:
>
> lm0 <- lm(rat ~ 1, weights=pop2000, data=oregon.tract1)
>
> I'm not advising this, but running localmoran.sad on this model output
> yields SAD p-values < 0.05 after FDR correction only in contiguous tracts
> on the Washington state line in Portland between the Columbia and
> Williamette rivers. So do look at the variables you are using before
> rushing into things.
>
> Hope this clarifies,
>
> Roger
>
>
>> best,
>>
>> Rafael HM Pereira
>> http://urbandemographics.blogspot.com
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> 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.
> voice: +47 55 95 93 55; e-mail: [hidden email]
> Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
> http://orcid.org/0000-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
        [[alternative HTML version deleted]]

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

Error-Mapping

Tue, 07/25/2017 - 11:32
Dear all,

I am a R- starter and not really experienced, but I try my best. For  
my PhD thesis I want to create an error map of a palaeo-relief and I  
found the instruction from Tomislav Hengl 2009 in "A practical Guide  
to Geostatistical mapping" (page 221 ff) really interesting. He also  
recommended this mailing list.

Actually I have two questions.

1. Question: Is it possible to create an empty grid with a individual  
shape of the research area? I found several instruction at the  
Internet but nothing really useful. Useful in this case means, that I  
can further follow the script from Hengl.

2. Question: The gstat package developed further and the function  
"overlay", which he used, is not available anymore. Instead, I think,  
the function "over" should be used, but then I don´t get the result  
which I need for further steps. Because the message reveals: Error in  
plot(cross.ov@coords[, 1], cross.ov@data[[1]], type = "l", xlab = "X",  
  : trying to get slot "coords" from an object of a basic class  
("numeric") with no slots. But of course, the initial objects had  
slots. If I got it right the function "over" just produces numeric  
values. What is the equivalent function for "overlay" from 2009?

I would be really happy about any hint or maybe newer/better solutions  
for error mapping of DEM´s. In my case paloae-DEM´s.

Best regards,
Ulrike Grimm


--
Dipl. Geographin

Universität Leipzig
Institut für Geographie
Johannisallee 19a
D-04103 Leipzig

http://geographie.physgeo.uni-leipzig.de/landschaft/mitarbeiter/ulrike-grimm/

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

Re: bivariate spatial correlation in R

Mon, 07/24/2017 - 13:56
On Mon, 24 Jul 2017, Rafael Pereira wrote:

> Hi all,
>
> I would like to ask whether some you conducted bi-variate spatial
> correlation in R.
>
> I know the bi-variate Moran's I is not implemented in the spdep library.
> I left a question on SO but also wanted to hear if anyone if the mainlist
> have come across this.
> https://stackoverflow.com/questions/45177590/map-of-bivariate-spatial-correlation-in-r-bivariate-lisa
>
> I also know Roger Bivand has implemented the L index proposed by Lee (2001)
> in spdep, but I'm not I'm not sure whether the L local correlation
> coefficients can be interpreted the same way as the local Moran's I
> coefficients. I couldn't find any reference commenting on this issue. I
> would very much appreciate your thoughts this. In the SO question, and in the follow-up, your presumably throw-away
example makes fundamental mistakes. The code in spdep by Virgilio
Gómez-Rubio is for uni- and bivariate L, and produces point values of
local L. This isn't the main problem, which is rather that you are not
taking account of the underlying population counts, nor shrinking any
estimates of significance to accommodate population sizes. Population
sizes vary from 0 to 11858, with the lower quartile at 3164 and upper
5698: plot(ecdf(oregon.tract$pop2000)). Should you be comparing rates in
stead? These are also compositional variables (sum to pop2000, or 1 if
rates) with the other missing components. You would probably be better
served by tools examining spatial segregation, such as for example the seg
package.

The 0 count populations cause problems for an unofficial alternative, the
black/white ratio:

oregon.tract1 <- oregon.tract[oregon.tract$white > 0,]
oregon.tract1$rat <- oregon.tract1$black/oregon.tract1$white
nb <- poly2nb(oregon.tract1)
lw <- nb2listw(nb)

which should still be adjusted by weighting:

lm0 <- lm(rat ~ 1, weights=pop2000, data=oregon.tract1)

I'm not advising this, but running localmoran.sad on this model output
yields SAD p-values < 0.05 after FDR correction only in contiguous tracts
on the Washington state line in Portland between the Columbia and
Williamette rivers. So do look at the variables you are using before
rushing into things.

Hope this clarifies,

Roger

>
> best,
>
> Rafael HM Pereira
> http://urbandemographics.blogspot.com
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> 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.
voice: +47 55 95 93 55; e-mail: [hidden email]
Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
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

bivariate spatial correlation in R

Mon, 07/24/2017 - 09:42
Hi all,

I would like to ask whether some you conducted bi-variate spatial
correlation in R.

I know the bi-variate Moran's I is not implemented in the spdep library. I
left a question on SO but also wanted to hear if anyone if the mainlist
have come across this.
https://stackoverflow.com/questions/45177590/map-of-bivariate-spatial-correlation-in-r-bivariate-lisa

I also know Roger Bivand has implemented the L index proposed by Lee (2001)
in spdep, but I'm not I'm not sure whether the L local correlation
coefficients can be interpreted the same way as the local Moran's I
coefficients. I couldn't find any reference commenting on this issue. I
would very much appreciate your thoughts this.

best,

Rafael HM Pereira
http://urbandemographics.blogspot.com

        [[alternative HTML version deleted]]

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

Re: Fwd: Re: Adding text to rasterVis levelplot

Fri, 07/21/2017 - 19:10
Thank you very much you both Thiago and Florian for your help.

Just for future searches on Google, the minimal example that allows
creating the plot i was looking for was:

------ START -----
library(rasterVis)

f <- system.file("external/test.grd", package="raster")
r <- raster(f)
s <- stack(r, r+500, r-500, r+200)

col.titles = c('col1','col2','','')
row.titles = c('row1','row2')

levelplot(s, layout=c(2,2),
          names.attr=col.titles,
          ylab=row.titles,
          scales = list(y = list(rot = 90) )
          )

------ END -----

However, the example provided by Florian opened up new possibilities
with the grid package for more customisation.

Kind regards,

Mauricio

Mauricio Zambrano-Bigiarini, PhD

=====================================
Department of Civil Engineering
Faculty of Engineering and Sciences
Universidad de La Frontera, Temuco, Chile
=====================================
"The ultimate inspiration is the deadline"
(Nolan Bushnell)
=====================================
Linux user #454569 -- Linux Mint user


On 21 July 2017 at 04:21, Florian Detsch
<[hidden email]> wrote:
> Ah, I forgot to include the definition of text vectors which must be
> inserted before the for() loop - sorry for that.
>
> ---
> ## labels
> cls <- c("col1", "col2")
> rws <- c("row1", "row2")
> ---
>
>
> -------- Forwarded Message --------
> Subject:        Re: [R-sig-Geo] Adding text to rasterVis levelplot
> Date:   Fri, 21 Jul 2017 10:05:08 +0200
> From:   Florian Detsch <[hidden email]>
> To:     [hidden email]
>
>
>
> Hi Mauricio,
>
> Complementing Thiago's answer using 'names.att', here is a more manual
> approach which might come in handy when a finer control over text
> annotations is required. It is built upon lattice::trellis.focus() which
> lets you select relevant panels of your trellis graph in an iterative
> manner.
>
> ---
> p <- levelplot(s, layout = c(2, 2), names.att = rep("", 4),
>                 scales = list(y = list(rot = 90)))
>
> library(grid)
>
> png("~/rasterVis.png", width = 14, height = 16, units = "cm", res = 300L)
> grid.newpage()
> print(p, newpage = FALSE)
>
> ## loop over panels to be labelled (ie 1:3)
> panels <- trellis.currentLayout()
> for (i in 1:3) {
>
>    # focus on current panel of interest and disable clipping
>    ids <- which(panels == i, arr.ind = TRUE)
>    trellis.focus("panel", ids[2], ids[1], clip.off = TRUE)
>
>    # add labels
>    if (i %in% c(1, 3)) {
>      if (i == 1) {
>        grid.text(cls[1], x = .5, y = 1.1) # add 'col1'
>        grid.text(rws[1], x = -.35, y = .5, rot = 90) # add 'row1'
>      } else {
>        grid.text(rws[2], x = -.35, y = .5, rot = 90) # add 'row2'
>      }
>    } else {
>      grid.text(cls[2], x = .5, y = 1.1) # add 'col2'
>    }
>
>    trellis.unfocus()
> }
>
> dev.off()
> ---
>
> A very similar use case, from which the above approach definitely took
> some inspiration, can be found here:
> https://stackoverflow.com/questions/7649737/is-it-possible-to-update-a-lattice-panel-in-r.
> Have a look at ?grid.text() and ?gpar() to learn more about
> customization options.
>
> Cheers,
> Florian
>
> --
> Dr. Florian Detsch
> Environmental Informatics
> Department of Geography
> Philipps-Universität Marburg
> Deutschhausstraße 12
> 35032 (parcel post: 35037) Marburg, Germany
>
> Phone: +49 (0) 6421 28-25323
> Web: http://www.uni-marburg.de/fb19/fachgebiete/umweltinformatik/detschf/index.html
>
>
> On 20.07.2017 05:07, Thiago V. dos Santos via R-sig-Geo wrote:
>> You can manually define the names that you want and pass them as arguments to 'names.attr' (the name of each panel) and 'ylab' (the title of the y axis).
>> Following your example:
>> col.titles = c('col1','col2','','')
>> row.titles = c('row1','row2')
>>
>> levelplot(s, layout=c(2,2),
>>            names.attr=col.titles,
>>            ylab=row.titles)
>>
>> HTH,
>> Greetings, -- Thiago V. dos Santos
>> PhD studentLand and Atmospheric ScienceUniversity of Minnesota
>>
>>
>> On Wednesday, July 19, 2017, 11:13:31 AM CDT, Mauricio Zambrano Bigiarini <[hidden email]> wrote:
>>
>> Dear all,
>>
>> I would like to add custom text labels for the columns and rows of a
>> levelplot object created with he rasterVis package, similar to the
>> c("Fall", "Winter", "Spring", Summer") used to label the rows and the
>> c("a",..., "h") used to label the columns in the plot shown int he
>> following link:
>>
>> https://stackoverflow.com/questions/43007241/how-to-add-text-to-a-specific-fixed-location-in-rastervis-levelplot
>>
>> In particular, how can I add c("col1", "col2") and c("row1", "row2")
>> as column and row headers, respectively, to the following example
>> plot:
>>
>> --------- START--------
>> library(rasterVis)
>> f <- system.file("external/test.grd", package="raster")
>> r <- raster(f)
>> s <- stack(r, r+500, r-500, r+200)
>> levelplot(s, layout=c(2,2), names.att=rep("",4))
>> --------- END--------
>>
>>
>> Thanks in advance for your help,
>>
>>
>> Mauricio Zambrano-Bigiarini, PhD
>>
>> =====================================
>> Department of Civil Engineering
>> Faculty of Engineering and Sciences
>> Universidad de La Frontera, Temuco, Chile
>> =====================================
>> "The ultimate inspiration is the deadline"
>> (Nolan Bushnell)
>> =====================================
>> Linux user #454569 -- Linux Mint user
>>
>
> _______________________________________________
> 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
--
La información contenida en este correo electrónico y cualquier anexo o
respuesta relacionada, puede contener datos e información confidencial y no
puede ser usada o difundida por personas distintas a su(s) destinatario(s).
Si usted no es el destinatario de esta comunicación, le informamos que
cualquier divulgación, distribución o copia de esta información constituye
un delito conforme a la ley chilena. Si lo ha recibido por error, por favor
borre el mensaje y todos sus anexos y notifique al remitente.

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

Re: RQGIS : looking for a geoalgorithm

Fri, 07/21/2017 - 04:19
Hi Tristan,

Perhaps qtm() in package tmap? It can automatically generate a list of maps
based on different attributes of a unique shapefile of polygons.

If you work with rasters you can create a raster stack and plot all rasters
at once.

HTH
Jerome


2017-07-20 15:02 GMT+02:00 Tristan Bourgeois <[hidden email]>:
>
> Dear all,
>
> Currently working on a new project I would like to create a dynamic report
> with R (using probably XML package for writting)
>
> In this report I'll plot some graphs (I'm ok with ggplots2) but also map.
>
>
> I just discovered the RQGIS package and I think I understood its functions.
> My question is the following one :
>
> Is there any geoalgorithms which allows to create an atlas ?
>
> I searched via the "find_algorithms()" but found nothing.
>
> Thanks in advance !
>
> Tristan.
>
> --
> Tristan Bourgeois
>
>         [[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

Fwd: Re: Adding text to rasterVis levelplot

Fri, 07/21/2017 - 03:21
Ah, I forgot to include the definition of text vectors which must be
inserted before the for() loop - sorry for that.

---
## labels
cls <- c("col1", "col2")
rws <- c("row1", "row2")
---


-------- Forwarded Message --------
Subject: Re: [R-sig-Geo] Adding text to rasterVis levelplot
Date: Fri, 21 Jul 2017 10:05:08 +0200
From: Florian Detsch <[hidden email]>
To: [hidden email]



Hi Mauricio,

Complementing Thiago's answer using 'names.att', here is a more manual
approach which might come in handy when a finer control over text
annotations is required. It is built upon lattice::trellis.focus() which
lets you select relevant panels of your trellis graph in an iterative
manner.

---
p <- levelplot(s, layout = c(2, 2), names.att = rep("", 4),
                scales = list(y = list(rot = 90)))

library(grid)

png("~/rasterVis.png", width = 14, height = 16, units = "cm", res = 300L)
grid.newpage()
print(p, newpage = FALSE)

## loop over panels to be labelled (ie 1:3)
panels <- trellis.currentLayout()
for (i in 1:3) {

   # focus on current panel of interest and disable clipping
   ids <- which(panels == i, arr.ind = TRUE)
   trellis.focus("panel", ids[2], ids[1], clip.off = TRUE)

   # add labels
   if (i %in% c(1, 3)) {
     if (i == 1) {
       grid.text(cls[1], x = .5, y = 1.1) # add 'col1'
       grid.text(rws[1], x = -.35, y = .5, rot = 90) # add 'row1'
     } else {
       grid.text(rws[2], x = -.35, y = .5, rot = 90) # add 'row2'
     }
   } else {
     grid.text(cls[2], x = .5, y = 1.1) # add 'col2'
   }

   trellis.unfocus()
}

dev.off()
---

A very similar use case, from which the above approach definitely took
some inspiration, can be found here:
https://stackoverflow.com/questions/7649737/is-it-possible-to-update-a-lattice-panel-in-r.
Have a look at ?grid.text() and ?gpar() to learn more about
customization options.

Cheers,
Florian

--
Dr. Florian Detsch
Environmental Informatics
Department of Geography
Philipps-Universität Marburg
Deutschhausstraße 12
35032 (parcel post: 35037) Marburg, Germany

Phone: +49 (0) 6421 28-25323
Web: http://www.uni-marburg.de/fb19/fachgebiete/umweltinformatik/detschf/index.html


On 20.07.2017 05:07, Thiago V. dos Santos via R-sig-Geo wrote:
> You can manually define the names that you want and pass them as arguments to 'names.attr' (the name of each panel) and 'ylab' (the title of the y axis).
> Following your example:
> col.titles = c('col1','col2','','')
> row.titles = c('row1','row2')
>
> levelplot(s, layout=c(2,2),
>            names.attr=col.titles,
>            ylab=row.titles)
>
> HTH,
> Greetings, -- Thiago V. dos Santos
> PhD studentLand and Atmospheric ScienceUniversity of Minnesota
>
>
> On Wednesday, July 19, 2017, 11:13:31 AM CDT, Mauricio Zambrano Bigiarini <[hidden email]> wrote:
>
> Dear all,
>
> I would like to add custom text labels for the columns and rows of a
> levelplot object created with he rasterVis package, similar to the
> c("Fall", "Winter", "Spring", Summer") used to label the rows and the
> c("a",..., "h") used to label the columns in the plot shown int he
> following link:
>
> https://stackoverflow.com/questions/43007241/how-to-add-text-to-a-specific-fixed-location-in-rastervis-levelplot
>
> In particular, how can I add c("col1", "col2") and c("row1", "row2")
> as column and row headers, respectively, to the following example
> plot:
>
> --------- START--------
> library(rasterVis)
> f <- system.file("external/test.grd", package="raster")
> r <- raster(f)
> s <- stack(r, r+500, r-500, r+200)
> levelplot(s, layout=c(2,2), names.att=rep("",4))
> --------- END--------
>
>
> Thanks in advance for your help,
>
>
> Mauricio Zambrano-Bigiarini, PhD
>
> =====================================
> Department of Civil Engineering
> Faculty of Engineering and Sciences
> Universidad de La Frontera, Temuco, Chile
> =====================================
> "The ultimate inspiration is the deadline"
> (Nolan Bushnell)
> =====================================
> Linux user #454569 -- Linux Mint user
>
_______________________________________________
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: Adding text to rasterVis levelplot

Fri, 07/21/2017 - 03:05
Hi Mauricio,

Complementing Thiago's answer using 'names.att', here is a more manual
approach which might come in handy when a finer control over text
annotations is required. It is built upon lattice::trellis.focus() which
lets you select relevant panels of your trellis graph in an iterative
manner.

---
p <- levelplot(s, layout = c(2, 2), names.att = rep("", 4),
                scales = list(y = list(rot = 90)))

library(grid)

png("~/rasterVis.png", width = 14, height = 16, units = "cm", res = 300L)
grid.newpage()
print(p, newpage = FALSE)

## loop over panels to be labelled (ie 1:3)
panels <- trellis.currentLayout()
for (i in 1:3) {

   # focus on current panel of interest and disable clipping
   ids <- which(panels == i, arr.ind = TRUE)
   trellis.focus("panel", ids[2], ids[1], clip.off = TRUE)

   # add labels
   if (i %in% c(1, 3)) {
     if (i == 1) {
       grid.text(cls[1], x = .5, y = 1.1) # add 'col1'
       grid.text(rws[1], x = -.35, y = .5, rot = 90) # add 'row1'
     } else {
       grid.text(rws[2], x = -.35, y = .5, rot = 90) # add 'row2'
     }
   } else {
     grid.text(cls[2], x = .5, y = 1.1) # add 'col2'
   }

   trellis.unfocus()
}

dev.off()
---

A very similar use case, from which the above approach definitely took
some inspiration, can be found here:
https://stackoverflow.com/questions/7649737/is-it-possible-to-update-a-lattice-panel-in-r.
Have a look at ?grid.text() and ?gpar() to learn more about
customization options.

Cheers,
Florian

--
Dr. Florian Detsch
Environmental Informatics
Department of Geography
Philipps-Universität Marburg
Deutschhausstraße 12
35032 (parcel post: 35037) Marburg, Germany

Phone: +49 (0) 6421 28-25323
Web: http://www.uni-marburg.de/fb19/fachgebiete/umweltinformatik/detschf/index.html


On 20.07.2017 05:07, Thiago V. dos Santos via R-sig-Geo wrote:
> You can manually define the names that you want and pass them as arguments to 'names.attr' (the name of each panel) and 'ylab' (the title of the y axis).
> Following your example:
> col.titles = c('col1','col2','','')
> row.titles = c('row1','row2')
>
> levelplot(s, layout=c(2,2),
>            names.attr=col.titles,
>            ylab=row.titles)
>
> HTH,
> Greetings, -- Thiago V. dos Santos
> PhD studentLand and Atmospheric ScienceUniversity of Minnesota
>
>
> On Wednesday, July 19, 2017, 11:13:31 AM CDT, Mauricio Zambrano Bigiarini <[hidden email]> wrote:
>
> Dear all,
>
> I would like to add custom text labels for the columns and rows of a
> levelplot object created with he rasterVis package, similar to the
> c("Fall", "Winter", "Spring", Summer") used to label the rows and the
> c("a",..., "h") used to label the columns in the plot shown int he
> following link:
>
> https://stackoverflow.com/questions/43007241/how-to-add-text-to-a-specific-fixed-location-in-rastervis-levelplot
>
> In particular, how can I add c("col1", "col2") and c("row1", "row2")
> as column and row headers, respectively, to the following example
> plot:
>
> --------- START--------
> library(rasterVis)
> f <- system.file("external/test.grd", package="raster")
> r <- raster(f)
> s <- stack(r, r+500, r-500, r+200)
> levelplot(s, layout=c(2,2), names.att=rep("",4))
> --------- END--------
>
>
> Thanks in advance for your help,
>
>
> Mauricio Zambrano-Bigiarini, PhD
>
> =====================================
> Department of Civil Engineering
> Faculty of Engineering and Sciences
> Universidad de La Frontera, Temuco, Chile
> =====================================
> "The ultimate inspiration is the deadline"
> (Nolan Bushnell)
> =====================================
> Linux user #454569 -- Linux Mint user
>
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

RQGIS : looking for a geoalgorithm

Thu, 07/20/2017 - 08:02
Dear all,

Currently working on a new project I would like to create a dynamic report
with R (using probably XML package for writting)

In this report I'll plot some graphs (I'm ok with ggplots2) but also map.


I just discovered the RQGIS package and I think I understood its functions.
My question is the following one :

Is there any geoalgorithms which allows to create an atlas ?

I searched via the "find_algorithms()" but found nothing.

Thanks in advance !

Tristan.

--
Tristan Bourgeois

        [[alternative HTML version deleted]]

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

Re: Adding text to rasterVis levelplot

Wed, 07/19/2017 - 22:07
You can manually define the names that you want and pass them as arguments to 'names.attr' (the name of each panel) and 'ylab' (the title of the y axis).
Following your example:
col.titles = c('col1','col2','','')
row.titles = c('row1','row2')

levelplot(s, layout=c(2,2),
          names.attr=col.titles,
          ylab=row.titles)

HTH,
Greetings, -- Thiago V. dos Santos
PhD studentLand and Atmospheric ScienceUniversity of Minnesota


On Wednesday, July 19, 2017, 11:13:31 AM CDT, Mauricio Zambrano Bigiarini <[hidden email]> wrote:

Dear all,

I would like to add custom text labels for the columns and rows of a
levelplot object created with he rasterVis package, similar to the
c("Fall", "Winter", "Spring", Summer") used to label the rows and the
c("a",..., "h") used to label the columns in the plot shown int he
following link:

https://stackoverflow.com/questions/43007241/how-to-add-text-to-a-specific-fixed-location-in-rastervis-levelplot

In particular, how can I add c("col1", "col2") and c("row1", "row2")
as column and row headers, respectively, to the following example
plot:

--------- START--------
library(rasterVis)
f <- system.file("external/test.grd", package="raster")
r <- raster(f)
s <- stack(r, r+500, r-500, r+200)
levelplot(s, layout=c(2,2), names.att=rep("",4))
--------- END--------


Thanks in advance for your help,


Mauricio Zambrano-Bigiarini, PhD

=====================================
Department of Civil Engineering
Faculty of Engineering and Sciences
Universidad de La Frontera, Temuco, Chile
=====================================
"The ultimate inspiration is the deadline"
(Nolan Bushnell)
=====================================
Linux user #454569 -- Linux Mint user

--
La información contenida en este correo electrónico y cualquier anexo o
respuesta relacionada, puede contener datos e información confidencial y no
puede ser usada o difundida por personas distintas a su(s) destinatario(s).
Si usted no es el destinatario de esta comunicación, le informamos que
cualquier divulgación, distribución o copia de esta información constituye
un delito conforme a la ley chilena. Si lo ha recibido por error, por favor
borre el mensaje y todos sus anexos y notifique al remitente.

_______________________________________________
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: raster - unrotate?

Wed, 07/19/2017 - 19:47
That's great, thanks for the update! I didn't actually realize all that
oceandata stuff was on Thredds. You might try ROpenSci's rerddap package
too.

I will try that source with the in-dev https://Github.com/hypertidy/tidync
. Please check it out if you can, and are feeling brave!

(tidync aims to provide more general outputs than raster can, more easily
than the API wrappers.)

Cheers, Mike

On Thu, 20 Jul 2017, 10:39 Ben Tupper, <[hidden email]> wrote:

> Ahoy!
>
> My (current) solution is to user raster::merge().  Works a treat and fits
> my workflow nicely.  Since I am obtaining the raster data using OPeNDAP it
> is easy to get the two regions, rotate the extent of the eastern region,
> and then merge.  Below is a pseudo-example (that works!) as it crops a
> local file twice, once for each region, rather than using an OPeNDaP calls
> via `ncdf4::nccvar_get(start = blah, count = something)` to get the two
> regions.
>
>
> ### start
>
> library(raster)
> uri <- '
> https://oceandata.sci.gsfc.nasa.gov:443/opendap/MODISA/L3SMI/2016/001/A20160012016032.L3m_R32_SST_sst_9km.nc
> '
> R <- raster::raster(uri, varname = 'sst')
>
> R
> #class       : RasterLayer
> #dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
> #resolution  : 1, 1  (x, y)
> #extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
> #coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> #data source : in memory
> #names       : layer
> #values      : 1.482572e-05, 0.9999928  (min, max)
>
> par(mfrow = c(2,1))
> plot(R)
> lines(c(180, 100, 100, 180), c(80,80,0,0))
> lines(c(-180,-90,-90,-180), c(80,80,0,0))
>
> # eastern section
> eR <- raster::crop(R, c(100, 180, 0, 80))
>
> # rotate the extent
> eExtent <- raster::extent(eR)
> raster::extent(eR) <- c(-260, -180, eExtent[3:4])
>
> # western section
> wR <- raster::crop(R, c(-180, -90, 0, 80))
>
> # merge
> newR <- raster::merge(eR, wR)
> newR
>
> plot(newR)
>
> #### end
>
> Cheers,
> Ben
>
> > On Jun 22, 2017, at 7:26 PM, Michael Sumner <[hidden email]> wrote:
> >
> >
> >
> >
> > On Thu, 22 Jun 2017 at 21:28 Ben Tupper <[hidden email]> wrote:
> > Hi,
> >
> > Wow!  A treasure from the past!
> >
> > Hmmm.   I see what you mean; it might be best to snip-clip-and-smush
> from the native presentation.  We are testing out the performance of the
> dineof function in the sinkr package (
> https://github.com/marchtaylor/sinkr/blob/master/R/dineof.R) to
> interpolate patchy chlorophyl data.
> >
> > Thanks for the tip!
> >
> >
> >
> > For what it's worth, I'm happy to help, there's no one size fits all
> here. The dateline is one of those "easy problems" that does not have a
> general fix and will bite in many different contexts. :)
> >
> > The best distillation of my tricks is in this project, but it does
> depend on your workflow and the actual data in use.
> >
> > https://github.com/hypertidy/tabularaster
> >
> > If you have performance issues raster's cell abstraction  will help, but
> they are a bit old-school when it comes to multi-part geometries. (It's
> identical to standard database indexing concepts as are all "spatial"
> optimization tricks)
> >
> > (I see a desperate need for a general API for this kind of problem so
> that we can build tools from a general framework, which I'm working on -
> that's some kind of excuse for why this and related projects are quite raw
> and unfinished.)
> >
> > Cheers, Mike.
> >
> > Cheers,
> > Ben
> >
> >> On Jun 22, 2017, at 4:46 AM, Michael Sumner <[hidden email]> wrote:
> >>
> >>
> >> It used to do the inverse, and I preserved a copy of the old behaviour
> here:
> >>
> >>
> https://github.com/AustralianAntarcticDivision/raadtools/blob/master/R/utils.R#L9
> >>
> >> You're probably best to learn from how it approaches it rather than
> just use it, since it's not a general capability, just works for those very
> specific cases.
> >>
> >> I just tested it and it still seems to work fine, I used
> >>
> >>
> oceandata.sci.gsfc.nasa.gov/MODISA/Mapped/Monthly/9km/chlor/A20021822002212.L3m_MO_CHL_chlor_a_9km.nc
> >>
> >> obtainable with
> >>
> >>
> https://oceandata.sci.gsfc.nasa.gov/cgi/getfile/A20021822002212.L3m_MO_CHL_chlor_a_9km.nc
> >>
> >> If you're extracting points you might as well just convert them into
> the "Atlantic" range, but it's painful to get that right for polygons or
> lines (and very hard to generalize once you get it working anyway).
> >>
> >> For simple regions like the one you show I'd probably split it into two
> extents() and use those - but depends on your workflow, all these options
> are useful here and there.
> >>
> >> Cheers, Mike.
> >>
> >>
> >> On Thu, 22 Jun 2017 at 18:30 Ben Tupper <[hidden email]> wrote:
> >> Hello,
> >>
> >> We have rasters that span [-180, 180] from NASA's Ocean Color (
> https://oceancolor.gsfc.nasa.gov/)  datasets.  We are trying to extract a
> region that spans 100E to 90W, that is 100 to -90.  The region 'wraps'
> across the edges as shown by the plot at the address below.
> >>
> >> https://dl.dropboxusercontent.com/u/8433654/sst.png
> >>
> >>
> >> library(raster)
> >> uri <- '
> https://oceandata.sci.gsfc.nasa.gov:443/opendap/MODISA/L3SMI/2016/001/A20160012016032.L3m_R32_SST_sst_9km.nc
> '
> >> R <- raster::raster(uri, varname = 'sst')
> >>
> >> R
> >> #class       : RasterLayer
> >> #dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
> >> #resolution  : 1, 1  (x, y)
> >> #extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
> >> #coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> >> #data source : in memory
> >> #names       : layer
> >> #values      : 1.482572e-05, 0.9999928  (min, max)
> >>
> >> plot(R)
> >> lines(c(180, 100, 100, 180), c(80,80,0,0))
> >> lines(c(-180,-90,-90,-180), c(80,80,0,0))
> >>
> >> I see that there is the nice raster::rotate() function to rotate raster
> coordinates from [0,360] to [-180,180]. That would make extracting the
> region super easy if the inverse were available.  Is there an equivalent
> way to implement the inverse or raster::rotate()?  I could dig into the
> source of raster::rotate() to see how to code one up, but I hate like heck
> to reinvent the wheel.
> >>
> >> Thanks!
> >> Ben
> >>
> >>
> >> Ben Tupper
> >> Bigelow Laboratory for Ocean Sciences
> >> 60 Bigelow Drive, P.O. Box 380
> >> East Boothbay, Maine 04544
> >> http://www.bigelow.org
> >>
> >> _______________________________________________
> >> 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
> >>
> >
> > Ben Tupper
> > Bigelow Laboratory for Ocean Sciences
> > 60 Bigelow Drive, P.O. Box 380
> > East Boothbay, Maine 04544
> > http://www.bigelow.org
> >
> >
> >
> > --
> > Dr. Michael Sumner
> > Software and Database Engineer
> > Australian Antarctic Division
> > 203 Channel Highway
> > Kingston Tasmania 7050 Australia
>
> Ben Tupper
> Bigelow Laboratory for Ocean Sciences
> 60 Bigelow Drive, P.O. Box 380
> East Boothbay, Maine 04544
> http://www.bigelow.org
>
> Ecocast Reports: http://seascapemodeling.org/ecocast.html
> Tick Reports: https://report.bigelow.org/tick/
> Jellyfish Reports: https://jellyfish.bigelow.org/jellyfish/
>
>
>
> -- 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

Re: raster - unrotate?

Wed, 07/19/2017 - 19:39
Ahoy!

My (current) solution is to user raster::merge().  Works a treat and fits my workflow nicely.  Since I am obtaining the raster data using OPeNDAP it is easy to get the two regions, rotate the extent of the eastern region, and then merge.  Below is a pseudo-example (that works!) as it crops a local file twice, once for each region, rather than using an OPeNDaP calls via `ncdf4::nccvar_get(start = blah, count = something)` to get the two regions.


### start

library(raster)
uri <- 'https://oceandata.sci.gsfc.nasa.gov:443/opendap/MODISA/L3SMI/2016/001/A20160012016032.L3m_R32_SST_sst_9km.nc'
R <- raster::raster(uri, varname = 'sst')

R
#class       : RasterLayer
#dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
#resolution  : 1, 1  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
#data source : in memory
#names       : layer
#values      : 1.482572e-05, 0.9999928  (min, max)

par(mfrow = c(2,1))
plot(R)
lines(c(180, 100, 100, 180), c(80,80,0,0))
lines(c(-180,-90,-90,-180), c(80,80,0,0))

# eastern section
eR <- raster::crop(R, c(100, 180, 0, 80))

# rotate the extent
eExtent <- raster::extent(eR)
raster::extent(eR) <- c(-260, -180, eExtent[3:4])

# western section
wR <- raster::crop(R, c(-180, -90, 0, 80))

# merge
newR <- raster::merge(eR, wR)
newR

plot(newR)

#### end

Cheers,
Ben

> On Jun 22, 2017, at 7:26 PM, Michael Sumner <[hidden email]> wrote:
>
>
>
>
> On Thu, 22 Jun 2017 at 21:28 Ben Tupper <[hidden email]> wrote:
> Hi,
>
> Wow!  A treasure from the past!
>
> Hmmm.   I see what you mean; it might be best to snip-clip-and-smush from the native presentation.  We are testing out the performance of the dineof function in the sinkr package (https://github.com/marchtaylor/sinkr/blob/master/R/dineof.R) to interpolate patchy chlorophyl data.  
>
> Thanks for the tip!
>
>
>
> For what it's worth, I'm happy to help, there's no one size fits all here. The dateline is one of those "easy problems" that does not have a general fix and will bite in many different contexts. :)
>
> The best distillation of my tricks is in this project, but it does depend on your workflow and the actual data in use.
>
> https://github.com/hypertidy/tabularaster
>
> If you have performance issues raster's cell abstraction  will help, but they are a bit old-school when it comes to multi-part geometries. (It's identical to standard database indexing concepts as are all "spatial" optimization tricks)
>
> (I see a desperate need for a general API for this kind of problem so that we can build tools from a general framework, which I'm working on - that's some kind of excuse for why this and related projects are quite raw and unfinished.)
>
> Cheers, Mike.  
>
> Cheers,
> Ben
>
>> On Jun 22, 2017, at 4:46 AM, Michael Sumner <[hidden email]> wrote:
>>
>>
>> It used to do the inverse, and I preserved a copy of the old behaviour here:
>>
>> https://github.com/AustralianAntarcticDivision/raadtools/blob/master/R/utils.R#L9
>>
>> You're probably best to learn from how it approaches it rather than just use it, since it's not a general capability, just works for those very specific cases.
>>
>> I just tested it and it still seems to work fine, I used
>>
>> oceandata.sci.gsfc.nasa.gov/MODISA/Mapped/Monthly/9km/chlor/A20021822002212.L3m_MO_CHL_chlor_a_9km.nc
>>
>> obtainable with
>>
>> https://oceandata.sci.gsfc.nasa.gov/cgi/getfile/A20021822002212.L3m_MO_CHL_chlor_a_9km.nc 
>>
>> If you're extracting points you might as well just convert them into the "Atlantic" range, but it's painful to get that right for polygons or lines (and very hard to generalize once you get it working anyway).
>>
>> For simple regions like the one you show I'd probably split it into two extents() and use those - but depends on your workflow, all these options are useful here and there.
>>
>> Cheers, Mike.
>>
>>
>> On Thu, 22 Jun 2017 at 18:30 Ben Tupper <[hidden email]> wrote:
>> Hello,
>>
>> We have rasters that span [-180, 180] from NASA's Ocean Color (https://oceancolor.gsfc.nasa.gov/)  datasets.  We are trying to extract a region that spans 100E to 90W, that is 100 to -90.  The region 'wraps' across the edges as shown by the plot at the address below.
>>
>> https://dl.dropboxusercontent.com/u/8433654/sst.png
>>
>>
>> library(raster)
>> uri <- 'https://oceandata.sci.gsfc.nasa.gov:443/opendap/MODISA/L3SMI/2016/001/A20160012016032.L3m_R32_SST_sst_9km.nc'
>> R <- raster::raster(uri, varname = 'sst')
>>
>> R
>> #class       : RasterLayer
>> #dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
>> #resolution  : 1, 1  (x, y)
>> #extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
>> #coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
>> #data source : in memory
>> #names       : layer
>> #values      : 1.482572e-05, 0.9999928  (min, max)
>>
>> plot(R)
>> lines(c(180, 100, 100, 180), c(80,80,0,0))
>> lines(c(-180,-90,-90,-180), c(80,80,0,0))
>>
>> I see that there is the nice raster::rotate() function to rotate raster coordinates from [0,360] to [-180,180]. That would make extracting the region super easy if the inverse were available.  Is there an equivalent way to implement the inverse or raster::rotate()?  I could dig into the source of raster::rotate() to see how to code one up, but I hate like heck to reinvent the wheel.
>>
>> Thanks!
>> Ben
>>
>>
>> Ben Tupper
>> Bigelow Laboratory for Ocean Sciences
>> 60 Bigelow Drive, P.O. Box 380
>> East Boothbay, Maine 04544
>> http://www.bigelow.org
>>
>> _______________________________________________
>> 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
>>
>
> Ben Tupper
> Bigelow Laboratory for Ocean Sciences
> 60 Bigelow Drive, P.O. Box 380
> East Boothbay, Maine 04544
> http://www.bigelow.org
>
>
>
> --
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> Kingston Tasmania 7050 Australia
Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org

Ecocast Reports: http://seascapemodeling.org/ecocast.html
Tick Reports: https://report.bigelow.org/tick/
Jellyfish Reports: https://jellyfish.bigelow.org/jellyfish/

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

Pages