Re: [R-sig-Geo] Using CRS() method on SpatialPointDataFrame to project coordinates

[MacQueen, Don via R-sig-Geo]

First, it looks like prcp_sites_ll does not have CRS information, because 
@projargs is NA. You'll need something like

  proj4string(prcp_sites_ll) <- CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 
+no_defs +towgs84=0,0,0')

That string may have more than is necessary; you could also try

   CRS('+init=epsg:4326')

Otherwise, your spTransform() command looks ok to me.

(have you come across spatialreference.org? I find it useful. Example
   http://spatialreference.org/ref/epsg/wgs-84/
)

-Don

--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509
 
 

On 9/20/18, 7:29 AM, "R-sig-Geo on behalf of Rich Shepard" 
 wrote:

   I have a SpatialPointsDataFrame with geographic coordinates:

#> str(prcp_sites_ll)
Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
   ..@ data   :'data.frame':58 obs. of  3 variables:
   .. ..$ name : Factor w/ 58 levels "Blazed Alder",..: 1 2 3 4 5 6 7 8 
9 10 ...
   .. ..$ elev : num [1:58] 1112.5 159.1 162.2 45.4 57.3 ...
   .. ..$ mean_prcp: num [1:58] 0.362 0.155 0.16 0.12 0.128 ...
   ..@ coords.nrs : int [1:2] 2 3
   ..@ coords : num [1:58, 1:2] -122 -122 -122 -123 -123 ...
   .. ..- attr(*, "dimnames")=List of 2
   .. .. ..$ : NULL
   .. .. ..$ : chr [1:2] "easting" "northing"
   ..@ bbox   : num [1:2, 1:2] -122.8 45 -121.7 45.5
   .. ..- attr(*, "dimnames")=List of 2
   .. .. ..$ : chr [1:2] "easting" "northing"
   .. .. ..$ : chr [1:2] "min" "max"
   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
   .. .. ..@ projargs: chr NA

and I want to project this using CRS('+init=epsg:2838') to a new SPDF
called prcp_sites_lcc.

   I think the proper syntax would be:

prcp_sites_lcc <- spTransform(prcp_sites_ll,CRS('+init=epsg:2838'))

   Is this correct?

Regards,

Rich

___
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] Help with simple Map of US states to predefined regions

[MacQueen, Don via R-sig-Geo]

I know this is not a complete solution -- and it's a very different approach -- 
but it should at least show you one way to reliably get states colored by 
region.
I also left out Alaska and Hawaii.

require(sp)
require(rgdal)

## US Census Bureau Tiger file -- polygons of each US State
## try this URL for download
##   
https://www.census.gov/cgi-bin/geo/shapefiles/index.php?year=2017=States+%28and+equivalent%29

ustf <- readOGR('.', 'tl_2017_us_state', stringsAsFactors=FALSE)

## note, the Tiger file includes 6 additional territories
dim(ustf)
## [1] 56 14

## get rid of the extra six territories  (state.name comes with R)
cus <- subset(ustf, NAME %in% state.name)

## cheap rename
cus$state <- cus$NAME
cus$abb <- cus$STUSPS

## invent ridiculous groupings of states
cus$grp <- 'a'
cus$grp[11:20] <- 'b'
cus$grp[21:30] <- 'c'
cus$grp[31:40] <- 'd'
cus$grp[41:50] <- 'e'

## assign colors to the groups
cus$color <- 'red'
cus$color[cus$grp=='b'] <- 'green'
cus$color[cus$grp=='c'] <- 'blue'
cus$color[cus$grp=='d'] <- 'brown'
cus$color[cus$grp=='e'] <- 'cyan'

## exclude Alaska, Hawaii
cus <- subset(cus, !(state %in% c('Alaska','Hawaii')))

## get rid of extraneous variables (optional)
cus <- cus[ , c('state','REGION','abb', 'grp') ]

## plot colored by regions as defined in the Census Bureau Tiger file
plot(cus, col=cus$REGION, usePolypath=FALSE)
## color "1" is black, looks bad, do this instead
plot(cus, col=as.numeric(cus$REGION)+1, usePolypath=FALSE)
text(coordinates(cus), cus$abb, col='white', cex=0.75)

## colors specified by a color variable in the data
plot(cus, col=cus$color, usePolypath=FALSE)
text(coordinates(cus), cus$abb, col='white', cex=0.75)


--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509
 
 

On 9/12/18, 11:27 AM, "R-sig-Geo on behalf of Bill Poling" 
 wrote:

Hi

I have this df with three columns ProviderState, ProviderStateCode, 
ProviderRegion I wanted to use to create a simple 5 color map

I have reviewed fiftystater pkg and map pkg but not sure how to simply take 
these three columns and plot a simple 5 color map based on the Region the state 
is in?

After looking at these and trying to apply these ideas to my data

https://cran.r-project.org/web/packages/fiftystater/vignettes/fiftystater.html
https://cran.r-project.org/web/packages/maps/maps.pdf

I found tutorial at: 
https://uchicagoconsulting.wordpress.com/tag/r-ggplot2-maps-visualization/

I used the tutorial data and subset in my regions

So now I have come up with the 5 segmented maps and my question becomes how 
to put this all into one map of the US?


install.packages("maps")
library(maps)
library(ggplot2)

#load us map data
all_states <- map_data("state") View(all_states)
#plot all states with ggplot
p <- ggplot()
p <- p + geom_polygon( data=all_states, aes(x=long, y=lat, group = 
group),colour="white", fill="blue" )
p

#http://sape.inf.usi.ch/quick-reference/ggplot2/colour

#Pacificstates
Pacificstates <- subset(all_states, region %in% c( "alaska", "arizona", 
"california", "hawaii", "nevada", "oregon","washington") )
p <- ggplot()
p <- p + geom_polygon( data=Pacificstates, aes(x=long, y=lat, group = 
group),colour="white", fill="deepskyblue4" ) +
  labs(title = "Pacificstates")
p

#Frontierstates
Frontierstates <- subset(all_states, region %in% c( "colorado", "idaho", 
"kansas", "montana", "new mexico", "oklahoma","texas", "utah", "wyoming") )
p <- ggplot()
p <- p + geom_polygon( data=Frontierstates, aes(x=long, y=lat, group = 
group),colour="white", fill="dodgerblue1" ) +
  labs(title = "FrontierStates")
p

#Midweststates
Midweststates <- subset(all_states, region %in% c( "iowa", "illinois", 
"indiana", "michigan", "minnesota", "missouri","north dakota", "nebraska", 
"ohio","south dakota","wisconsin") )
p <- ggplot()
p <- p + geom_polygon( data=Midweststates, aes(x=long, y=lat, group = 
group),colour="white", fill="dodgerblue1" ) +
  labs(title = "MidwestStates")
p

#Southernstates
Southernstates <- subset(all_states, region %in% c( "alabama", "arkansas", 
"florida", "georgia", "kentucky", "louisiana","mississippi"
,"north carolina", 
"south carolina","tennessee","virginia","west virginia") )
p <- ggplot()
p <- p + geom_polygon( data=Southernstates, aes(x=long, y=lat, group = 
group),colour="white", fill="royalblue2" ) +
  labs(title = "Southernstates")
p

# Northeaststates
Northeaststates <- subset(all_states, region %in% c( "connecticut", 
"district of columbia", "delaware", "massachusetts", "maryland", "maine","new 
hampshire"
 , "new jersey", "new 
york","pennsylvania","rhode 

Re: [R-sig-Geo] Help with simple Map of US states with predefined regions Version 2

[MacQueen, Don via R-sig-Geo]

I know this is not a complete solution -- and it's a very different approach -- 
but it should at least show you a way to reliably get states colored by region.
(I also left out Alaska and Hawaii, since the point here is how to color the 
regions)

require(sp)
require(rgdal)

## US Census Bureau Tiger file -- polygons of each US State
## try this URL for download
##  
https://www.census.gov/cgi-bin/geo/shapefiles/index.php?year=2017=States+%28and+equivalent%29

## unzip to working directory ( '.' )
ustf <- readOGR('.', 'tl_2017_us_state', stringsAsFactors=FALSE)

## note, the Tiger file includes 6 additional territories
dim(ustf)
## [1] 56 14

## get rid of the extra six territories  (state.name comes with R)
cus <- subset(ustf, NAME %in% state.name)

## cheap rename
cus$state <- cus$NAME
cus$abb <- cus$STUSPS

## invent ridiculous groupings of states
cus$grp <- 'a'
cus$grp[11:20] <- 'b'
cus$grp[21:30] <- 'c'
cus$grp[31:40] <- 'd'
cus$grp[41:50] <- 'e'

## assign colors to the groups
cus$color <- 'red'
cus$color[cus$grp=='b'] <- 'green'
cus$color[cus$grp=='c'] <- 'blue'
cus$color[cus$grp=='d'] <- 'brown'
cus$color[cus$grp=='e'] <- 'cyan'

## exclude Alaska, Hawaii
cus <- subset(cus, !(state %in% c('Alaska','Hawaii')))

## get rid of extraneous variables (optional)
cus <- cus[ , c('state','REGION','abb', 'grp') ]

## plot colored by regions as defined in the Census Bureau Tiger file
plot(cus, col=cus$REGION, usePolypath=FALSE)

## color "1" is black, looks bad, do this instead
plot(cus, col=as.numeric(cus$REGION)+1, usePolypath=FALSE)
text(coordinates(cus), cus$abb, col='white', cex=0.75)

## colors specified by a color variable in the data
plot(cus, col=cus$color, usePolypath=FALSE)
text(coordinates(cus), cus$abb, col='white', cex=0.75)

(my preferred graphics device does not support Polypath, but probably most 
others do, so one can omit usePolypath=FALSE)

-Don

--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509
 
 

On 9/13/18, 5:17 AM, "R-sig-Geo on behalf of Bill Poling" 
 wrote:

Hi,

I hope someone can help me finalize this please.

I am coming close to what I need using variations from two ggplot2 
tutorials.

This first gives me the map of the US with AK & HI but I cannot figure out 
how to get my 5 regions colored


#https://stackoverflow.com/questions/38021188/how-to-draw-u-s-state-map-with-hi-and-ak-with-state-abbreviations-centered-us?rq=1


library(ggplot2)
install.packages("ggalt")
library(ggalt) # coord_proj
library(albersusa) # devtools::install_github("hrbrmstr/albersusa")
install.packages("ggthemes")
library(ggthemes)  # theme_map
install.packages("rgeos")
library(rgeos) # centroids
library(dplyr)

# composite map with AK & HI
usa_map <- usa_composite()

# calculate the centroids for each state gCentroid(usa_map, byid=TRUE) %>%
  as.data.frame() %>%
  mutate(state=usa_map@data$iso_3166_2) -> centroids

# make it usable in ggplot2
usa_map <- fortify(usa_map)

View(usa_map)
t1 <- head(usa_map,n=5)
knitr::kable(t1, row.names=FALSE, align=c("l", "l", "r", "r", "r"))

#

  # |long  |lat  | group| order|  region|subregion |
  # |:-|:|-:|-:|---:|:-|
  # |-87.46201 |30.38968 | 1| 1| alabama|NA|
  # |-87.48493 |30.37249 | 1| 2| alabama|NA|
  # |-87.52503 |30.37249 | 1| 3| alabama|NA|
  # |-87.53076 |30.33239 | 1| 4| alabama|NA|
  # |-87.57087 |30.32665 | 1| 5| alabama|NA|

usa_map <- fortify(usa_map)
gg <- ggplot()
gg <- gg + geom_map(data=usa_map, map=usa_map,
aes(long, lat, map_id=id),
color="#2b2b2b", size=0.1, fill=NA)

gg <- gg + geom_text(data=centroids, aes(x, y, label=state), size=2) gg <- 
gg + coord_proj(us_laea_proj) gg <- gg + theme_map() gg





#/

This second is an alternative (however not liking AK, not coming into 
the map like scenario one above) but also ignoring new Mexico (because 
recognizing a seventh field value) and I suspect it will do the same for new 
York and new jersey etc.. when I add them to the list.

Error in scan(file = file, what = what, sep = sep, quote = quote, dec = 
dec,  :  line 12 did not have 6 elements

When I use newmexico (all one word) it appears white in the map like the 
other states not in the table statement


#https://stackoverflow.com/questions/3832/r-code-to-generating-map-of-us-states-with-specific-colors

library(ggplot2)

Re: [R-sig-Geo] help: Problem getting centroids inside lists

[MacQueen, Don via R-sig-Geo]

If all of your data frames had enough points then this should work:

myfun1 <- function(le) {
  list(one=geosphere::centroid(coordinates(le$one)),
   two=geosphere::centroid(coordinates(le$two))
   )
}

lapply(ct, myfun1)

To handle the one point case, try the following, which I think works with your 
example data, but is untested if there just two points.
It also assumes the data frames are named 'one' and 'two', and will ignore any 
others. To handle other names, I think you
could replace  $one  with  [[1]]  and  $two  with  [[2]]  .

myfun <- function(le) {
  list(one=if (length(le$one)>1) geosphere::centroid(coordinates(le$one)) else 
coodinates(le$one),
   two=if (length(le$two)>1) geosphere::centroid(coordinates(le$two)) else 
coordinates(le$two)
   )
}

lapply(ct, myfun)


To see a little bit of what's going on, try

  myfun(ct[[1]])
  myfun1(ct[[1]])
  myfun1(ct[[2]])
  myfun(ct[[2]])



I would also verify that the length method for objects of class SpatialPoints 
does return the number of points, as it appears to:

> class(ct$a$two)
[1] "SpatialPoints"
attr(,"package")
[1] "sp"

> length(ct$a$one)
[1] 3
> length(ct$a$two)
[1] 3
> length(ct$a$two)
[1] 3
> length(ct$b$two)
[1] 1

--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509
 
 

On 9/10/18, 10:56 AM, "R-sig-Geo on behalf of Ariel Fuentesdi" 
 wrote:

Hi everyone,

I have a list of coordinates called "ct" and I want to extract the
centroids of each sublist, but it only works when it has only 3 or more
points. In ct$b$two only has one point; in that case, I would rescue that
point as If it were a centroid.

This is what I did:

library(dplyr)
library(geosphere)

ct <- list(a = list(one = data.frame(lon = c(-180, -160, -60), lat = c(-20,
5, 0)),
two = data.frame(lon = c(-18, -16, -6), lat = c(-2, 50,
10))),
   b = list(one = data.frame(lon = c(-9, -8, -3), lat = c(-1, 25,
5)),
two = data.frame(lon = c(-90), lat = c(-1

coordinates(ct$a$one) <- ~lon+lat
coordinates(ct$a$two) <- ~lon+lat
coordinates(ct$b$one) <- ~lon+lat
coordinates(ct$b$two) <- ~lon+lat

s <- 1:length(ct)
ctply <- list()
for (i in s){
 for (j in 1:length(ct[[i]])) {
  ctply[[i]][j] <- ifelse(test = lengths(ct[[i]][1]) > 2, yes = sapply(X =
ct[[i]][j],
FUN = function(y) geosphere::centroid(slot(y, "coords"))), no =
ct[[i]][j] )
 }
}

Thanks in advance

Regards,
Ariel

[[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] CRAN releases of sp, rgdal and rgeos

[MacQueen, Don via R-sig-Geo]

Success!

This morning I upgraded a Mac to OS 10.13.5 (the so-called High Sierra 
version), then

 - Installed R 3.5.0 from CRAN (installed ever 3.3.x)
 - Installed sp 1.3-1, rgdal 1.3-2, rgeos 0.3-28, also sf 0.6-3, from the CRAN 
binaries that are now available
 - Ran my personal test suite, which exercises the kinds of tasks I perform 
with those packages

All tests succeeded.

I forgot to control the order of installation, so I don't know if I installed 
sp first, as advised. But I expect that with the binary versions it doesn't 
matter.

(I don't think it matters for the above, but I also installed the clang and 
gfortran version provided on CRAN's Mac "tools" page, and successfully compiled 
some source packages that require fortran, and others that require C).

Thank again for all your work
-Don

--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509
 
 

On 6/8/18, 11:15 AM, "R-sig-Geo on behalf of Roger Bivand" 
 wrote:

There are new releases of sp, rgdal and rgeos on CRAN. Please install sp 
first, then the other two, which link to the installed sp. They all 
address so-called rchk issues, which have not so far been a problem, but 
might have become more fragile as R's internal memory management is made 
even more efficient. This involves compiled code using memory allocated by 
R to be freed by R's garbage collector, which has to know if an object is 
still being used. Tomas Kalibera, the author of rchk, helped resolve and 
explain the issues encountered - what was good coding practice fifteen 
years ago isn't always still good practice.

In addition, the earliest versions of GDAL and PROJ with which rgdal will 
work have been updated, and set to PROJ 4.8.0 and GDAL 1.11.4. The current 
released versions of PROJ and GDAL are to be prefered, as bugs have been 
fixed and new features and drivers introduced. A check has been put 
in place to trap attempts to install rgdal without a C++11-capable 
compiler when the GDAL version is >=2.3.0 - which requires C++11. rgeos is 
ready for the forthcoming version of GEOS.

The CRAN team has also been very supportive of our efforts to bring 
compiled code in these packages into rchk compliance.

Please get in touch if you see any loose ends in these releases.

Roger

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: roger.biv...@nhh.no
http://orcid.org/-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0J=en

___
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] Error using spsample on SpatialLines in Linux OS

[MacQueen, Don]

The code also works fine on my Mac, and an RHEL (linux) machine to which I have 
access (though both are still at R 3.4.2 and sp_1.2-5).

You could focus in on the problem a bit by trying the following (I expect you 
will see the same error)

> sp::LineLength(pts)
[1] 123.0096

Starting with spsample(), and drilling down, eventually LineLength() calls the 
C function sp_lengths.  LineLength() is pretty simple, and I don't see any 
evidence that it depends on other packages. I'm not the expert that Edzer is, 
but it pretty much looks like an installation problem with sp, sorry to say. I 
can't think of anything else...

Sorry I'm not more helpful,
-Don


> sp::LineLength(pts, sum=FALSE)
[1] 15.81139  7.81025 42.37924 57.00877

> SpatialLinesLengths(L)
[1] 123.0096

> LinesLength(L@lines[[1]])
[1] 123.0096


> sessionInfo()
R version 3.4.2 (2017-09-28)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: OS X El Capitan 10.11.6

Matrix products: default
BLAS: 
/Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: 
/Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base 

other attached packages:
[1] sp_1.2-5

loaded via a namespace (and not attached):
 [1] compiler_3.4.2 colorspace_1.3-2   scales_0.5.0.9000  lazyeval_0.2.1
 [5] plyr_1.8.4 rgdal_1.2-13   tools_3.4.2gtable_0.2.0  
 [9] tibble_1.3.4   Rcpp_0.12.14   ggplot2_2.2.1.9000 grid_3.4.2
[13] rlang_0.1.2munsell_0.4.3  lattice_0.20-35   

--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509
 
 

On 4/26/18, 9:52 AM, "R-sig-Geo on behalf of Glenn Stauffer" 
 wrote:

Do you mean the sp package? I should have mentioned that I already tried
that.

-Original Message-
From: R-sig-Geo [mailto:r-sig-geo-boun...@r-project.org] On Behalf Of Edzer
Pebesma
Sent: Thursday, April 26, 2018 11:35 AM
To: r-sig-geo@r-project.org
Subject: Re: [R-sig-Geo] Error using spsample on SpatialLines in Linux OS

I would try reinstalling the package.

On 04/26/2018 06:31 PM, Glenn Stauffer wrote:
> I ran across a problem when trying to use spsample to sample points 
> along a line, while in a Linux OS (Mint 18.3). The following code 
> (which works fine in Win 10):
> 
>  
> 
> pts <- matrix(c(c(6000,6015,6021,6035,6050),
> c(6000,5995,6000,6040,6095)),5,2,byrow=FALSE)
> L = 
> SpatialLines(list(Lines(list(Line(coordinates(pts))),"X")),proj4string 
> =
> CRS("+init=epsg:3071"))
> plot(L)
> Lsamp <- spsample(L,10,type="regular",offset=0.5) # this is the line 
> that generates the error
> 
>  
> 
> produces the following error:
> 
>  
> 
> Error in .C("sp_lengths", x, y, n, lengths, lonlat, PACKAGE = "sp") : 
>   "sp_lengths" not available for .C() for package "sp"
> 
> Could this be related to some wrong versions of certain dependencies 
> (e.g., GDAL)? I don't know much about that, but I did run into that 
> issue with trying to install other packages (rgeos, sf). In any case, 
> how can I prevent this error?
> 
>  
> 
> Thanks,
> 
> Glenn
> 
>  
> 
>  
> 
> 
>   [[alternative HTML version deleted]]
> 
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

___
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] Change in projection results relative to 3 years ago

[MacQueen, Don]

That solved the problem, thanks!

-Don

--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509



From: Roger Bivand <roger.biv...@nhh.no>
Date: Monday, October 30, 2017 at 8:54 AM
To: MacQueen_Don <macque...@llnl.gov>
Cc: "r-sig-geo@r-project.org" <r-sig-geo@r-project.org>
Subject: Re: [R-sig-Geo] Change in projection results relative to 3 years ago

Fill in from http://download.osgeo.org/proj/proj-datumgrid-1.6.zip, this was 
temporarily missing in the CRAN OSX binary.
Roger
Roger Bivand
Norwegian School of Economics
Bergen, Norway



On Mon, Oct 30, 2017 at 4:46 PM +0100, "MacQueen, Don" 
<macque...@llnl.gov<mailto:macque...@llnl.gov>> wrote:

I should have checked this before my last email:



Since I have the KyngChaos frameworks:



[69]% ls /Library/Frameworks/PROJ.framework/unix/share/proj/

./  WI  nad.lst 
proj_def.dat

../ WO  nad27   prvi

CH  alaska  nad83   stgeorge

FL  conus   ntf_r93.gsb stlrnc

GL27epsgntv1_can.datstpaul

IGNFesrinullworld

MD  esri.extra  nzgd2kgrid0005.gsb

TN  hawaii  other.extra



If my rgdal was installed from the CRAN binary, which I believe it was, then 
perhaps if I install from source it will pick up those files.



-Don



--

Don MacQueen

Lawrence Livermore National Laboratory

7000 East Ave., L-627

Livermore, CA 94550

925-423-1062

Lab cell 925-724-7509







On 10/30/17, 7:55 AM, "Roger Bivand"  wrote:



On Mon, 30 Oct 2017, MacQueen, Don wrote:



> Three years ago Roger and Hermann Peifer graciously helped me solve a 
projection issue in R. Something has since changed; results are different now.

>

> I'd once again appreciate assistance... and am more or less hoping that 
someone will have some idea about what has changed, hence what I can do to 
return to my previously blissful state.

>

> Thanks in advance,

> -Don

>



Briefly, please check the current contents of:



list.files("/Users/macqueen1/Library/R/3.4/library/rgdal/proj")



I get:



> locs.26743 <- spTransform(locs.ll, CRS("+init=epsg:26743"))

> coordinates(locs.ref)-coordinates(locs.26743)

  coords.x1 coords.x2

[1,]  3.746539 -1.876668

[2,]  3.746607 -1.876466



with



> library(rgdal)

Loading required package: sp

rgdal: version: 1.2-15, (SVN revision 691)

  Geospatial Data Abstraction Library extensions to R successfully loaded

  Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15

  Path to GDAL shared files: /usr/local/share/gdal

  GDAL binary built with GEOS: TRUE

  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]

  Path to PROJ.4 shared files: (autodetected)

  Linking to sp version: 1.2-5



(1.2-15 is the development version, but nothing changed in this respect

vis. 1.2-13)



and



> list.files("/usr/local/share/proj")

  [1] "alaska" "CH" "conus"

  [4] "epsg"   "esri"   "esri.extra"

  [7] "FL" "GL27"   "hawaii"

[10] "IGNF"   "MD" "nad.lst"

[13] "nad27"  "nad83"  "ntf_r93.gsb"

[16] "ntv1_can.dat"   "null"   "nzgd2kgrid0005.gsb"

[19] "other.extra""proj_def.dat"   "prvi"

[22] "stgeorge"   "stlrnc" "stpaul"

[25] "TN" "WI" "WO"

[28] "world"



Maybe you are being thrown back just to the ellipsoid if the NAD grids are

absent?



Hope this helps,



Roger



>

>

> My original request is here

>   https://stat.ethz.ch/pipermail/r-sig-geo/2014-August/021611.html

> and Roger's solution is here

>   https://stat.ethz.ch/pipermail/r-sig-geo/2014-August/021623.html

> (actually, they did more than help, the pretty much did all the work!)

>

>

> I recently discovered the coordinates produced by that solution are now

> roughly 100m different from what they were then.

> (I've done quite a bit of ch

Re: [R-sig-Geo] Change in projection results relative to 3 years ago

[MacQueen, Don]

I should have checked this before my last email:

Since I have the KyngChaos frameworks:

[69]% ls /Library/Frameworks/PROJ.framework/unix/share/proj/
./  WI  nad.lst 
proj_def.dat
../ WO  nad27   prvi
CH  alaska  nad83   stgeorge
FL  conus   ntf_r93.gsb stlrnc
GL27epsgntv1_can.datstpaul
IGNFesrinullworld
MD  esri.extra  nzgd2kgrid0005.gsb
TN  hawaii  other.extra

If my rgdal was installed from the CRAN binary, which I believe it was, then 
perhaps if I install from source it will pick up those files.

-Don

--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509
 
 

On 10/30/17, 7:55 AM, "Roger Bivand" <roger.biv...@nhh.no> wrote:

On Mon, 30 Oct 2017, MacQueen, Don wrote:

> Three years ago Roger and Hermann Peifer graciously helped me solve a 
projection issue in R. Something has since changed; results are different now.
>
> I'd once again appreciate assistance... and am more or less hoping that 
someone will have some idea about what has changed, hence what I can do to 
return to my previously blissful state.
>
> Thanks in advance,
> -Don
>

Briefly, please check the current contents of:

list.files("/Users/macqueen1/Library/R/3.4/library/rgdal/proj")

I get:

> locs.26743 <- spTransform(locs.ll, CRS("+init=epsg:26743"))
> coordinates(locs.ref)-coordinates(locs.26743)
  coords.x1 coords.x2
[1,]  3.746539 -1.876668
[2,]  3.746607 -1.876466

with

> library(rgdal)
Loading required package: sp
rgdal: version: 1.2-15, (SVN revision 691)
  Geospatial Data Abstraction Library extensions to R successfully loaded
  Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15
  Path to GDAL shared files: /usr/local/share/gdal
  GDAL binary built with GEOS: TRUE
  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
  Path to PROJ.4 shared files: (autodetected)
  Linking to sp version: 1.2-5

(1.2-15 is the development version, but nothing changed in this respect 
vis. 1.2-13)

and

> list.files("/usr/local/share/proj")
  [1] "alaska" "CH" "conus"
  [4] "epsg"   "esri"   "esri.extra"
  [7] "FL" "GL27"   "hawaii"
[10] "IGNF"   "MD" "nad.lst"
[13] "nad27"  "nad83"  "ntf_r93.gsb"
[16] "ntv1_can.dat"   "null"   "nzgd2kgrid0005.gsb"
[19] "other.extra""proj_def.dat"   "prvi"
[22] "stgeorge"   "stlrnc" "stpaul"
[25] "TN" "WI" "WO"
[28] "world"

Maybe you are being thrown back just to the ellipsoid if the NAD grids are 
absent?

Hope this helps,

Roger

>
>
> My original request is here
>   https://stat.ethz.ch/pipermail/r-sig-geo/2014-August/021611.html
> and Roger's solution is here
>   https://stat.ethz.ch/pipermail/r-sig-geo/2014-August/021623.html
> (actually, they did more than help, the pretty much did all the work!)
>
>
> I recently discovered the coordinates produced by that solution are now
> roughly 100m different from what they were then.
> (I've done quite a bit of checking to make sure I'm not making some
> dumb mistake, and feel confident I haven't. Hopefully I'm right!)
>
> A small set of example locations is transformed from long/lat to a local 
coordinate system using
>  (A) a reference method deemed correct
>  (B) an R method using sp::spTransform()
> The goal was to have B reproduce A. Three years ago, it did. Now it does 
not.
>
> (The accuracy of method A is not the issue. It uses standard methods 
implemented in ESRI software.)
>
> The reproducible example below (same one as 3 years ago) does not use the 
full version of method B, but it does illustrate how something has changed.
>
>
> ##
> ## reproducible example begins
> ##
>
>

Re: [R-sig-Geo] Change in projection results relative to 3 years ago

[MacQueen, Don]

I get

> list.files("/Users/macqueen1/Library/R/3.4/library/rgdal/proj")
 [1] "CH"   "GL27" "IGNF" "epsg" "esri"
 [6] "esri.extra"   "nad.lst"  "nad27""nad83""other.extra" 
[11] "proj_def.dat" "world"   

Considerably different indeed, and conus in particular is missing.

I have a somewhat vague memory of some r-sig-geo traffic some months ago that 
might be relevant.

Thanks
-Don 

--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509
 
 

On 10/30/17, 7:55 AM, "Roger Bivand" <roger.biv...@nhh.no> wrote:

On Mon, 30 Oct 2017, MacQueen, Don wrote:

> Three years ago Roger and Hermann Peifer graciously helped me solve a 
projection issue in R. Something has since changed; results are different now.
>
> I'd once again appreciate assistance... and am more or less hoping that 
someone will have some idea about what has changed, hence what I can do to 
return to my previously blissful state.
>
> Thanks in advance,
> -Don
>

Briefly, please check the current contents of:

list.files("/Users/macqueen1/Library/R/3.4/library/rgdal/proj")

I get:

> locs.26743 <- spTransform(locs.ll, CRS("+init=epsg:26743"))
> coordinates(locs.ref)-coordinates(locs.26743)
  coords.x1 coords.x2
[1,]  3.746539 -1.876668
[2,]  3.746607 -1.876466

with

> library(rgdal)
Loading required package: sp
rgdal: version: 1.2-15, (SVN revision 691)
  Geospatial Data Abstraction Library extensions to R successfully loaded
  Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15
  Path to GDAL shared files: /usr/local/share/gdal
  GDAL binary built with GEOS: TRUE
  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
  Path to PROJ.4 shared files: (autodetected)
  Linking to sp version: 1.2-5

(1.2-15 is the development version, but nothing changed in this respect 
vis. 1.2-13)

and

> list.files("/usr/local/share/proj")
  [1] "alaska" "CH" "conus"
  [4] "epsg"   "esri"   "esri.extra"
  [7] "FL" "GL27"   "hawaii"
[10] "IGNF"   "MD" "nad.lst"
[13] "nad27"  "nad83"  "ntf_r93.gsb"
[16] "ntv1_can.dat"   "null"   "nzgd2kgrid0005.gsb"
[19] "other.extra""proj_def.dat"   "prvi"
[22] "stgeorge"   "stlrnc" "stpaul"
[25] "TN" "WI" "WO"
[28] "world"

Maybe you are being thrown back just to the ellipsoid if the NAD grids are 
absent?

Hope this helps,

Roger

>
>
> My original request is here
>   https://stat.ethz.ch/pipermail/r-sig-geo/2014-August/021611.html
> and Roger's solution is here
>   https://stat.ethz.ch/pipermail/r-sig-geo/2014-August/021623.html
> (actually, they did more than help, the pretty much did all the work!)
>
>
> I recently discovered the coordinates produced by that solution are now
> roughly 100m different from what they were then.
> (I've done quite a bit of checking to make sure I'm not making some
> dumb mistake, and feel confident I haven't. Hopefully I'm right!)
>
> A small set of example locations is transformed from long/lat to a local 
coordinate system using
>  (A) a reference method deemed correct
>  (B) an R method using sp::spTransform()
> The goal was to have B reproduce A. Three years ago, it did. Now it does 
not.
>
> (The accuracy of method A is not the issue. It uses standard methods 
implemented in ESRI software.)
>
> The reproducible example below (same one as 3 years ago) does not use the 
full version of method B, but it does illustrate how something has changed.
>
>
> ##
> ## reproducible example begins
> ##
>
> ## define two example points in WGS84 long/lat (same as 3 years ago)
> locs.xy <- cbind(
>  c(-121.524764291826,-121.523480804667),
>  c(37.6600366036405,37.6543604613483)
> )
>
> locs.ll <- SpatialPoints(locs.xy, proj4string=CRS("+proj=longlat 
+datum=W

[R-sig-Geo] Change in projection results relative to 3 years ago

[MacQueen, Don]

Three years ago Roger and Hermann Peifer graciously helped me solve a 
projection issue in R. Something has since changed; results are different now.

I'd once again appreciate assistance... and am more or less hoping that someone 
will have some idea about what has changed, hence what I can do to return to my 
previously blissful state.

Thanks in advance,
-Don



My original request is here
   https://stat.ethz.ch/pipermail/r-sig-geo/2014-August/021611.html
and Roger's solution is here
   https://stat.ethz.ch/pipermail/r-sig-geo/2014-August/021623.html
(actually, they did more than help, the pretty much did all the work!)


I recently discovered the coordinates produced by that solution are now
roughly 100m different from what they were then.
(I've done quite a bit of checking to make sure I'm not making some
dumb mistake, and feel confident I haven't. Hopefully I'm right!)

A small set of example locations is transformed from long/lat to a local 
coordinate system using
  (A) a reference method deemed correct
  (B) an R method using sp::spTransform()
The goal was to have B reproduce A. Three years ago, it did. Now it does not.

(The accuracy of method A is not the issue. It uses standard methods 
implemented in ESRI software.)

The reproducible example below (same one as 3 years ago) does not use the full 
version of method B, but it does illustrate how something has changed.


##
## reproducible example begins
##

## define two example points in WGS84 long/lat (same as 3 years ago)
locs.xy <- cbind(
  c(-121.524764291826,-121.523480804667),
  c(37.6600366036405,37.6543604613483)
)

locs.ll <- SpatialPoints(locs.xy, proj4string=CRS("+proj=longlat +datum=WGS84") 
)


## Outside of R, reference method A converts from long/lat
## to a local system, EPSG 26743, which is:
##   California State Plane Zone III, NAD27, feet;
##   http://spatialreference.org/ref/epsg/26743/
## using an ESRI "two-step process" (some detail at the end),
## saved as a shapefile, loaded into R using readOGR().

## Here is the "dput" of those coordinates from three years ago (the example 
location transformed by method A)
locs.ref <- new(
  "SpatialPoints",
  coords = structure(c(1703671.30566227, 1704020.20113366,
   424014.398045834, 421943.708664294),
 .Dim = c(2L, 2L),
 .Dimnames = list(NULL, c("coords.x1", "coords.x2")))
, bbox = structure(
c(1703671.30566227, 421943.708664294,
  1704020.20113366, 424014.398045834),
.Dim = c(2L, 2L),
.Dimnames = list(c("coords.x1",  "coords.x2"), c("min", "max")))
, proj4string =
new("CRS",
projargs = "+proj=lcc +lat_1=37.07 +lat_2=38.43 
+lat_0=36.5 +lon_0=-120.5 +x_0=609601.2192024384 +y_0=0 +datum=NAD27 
+units=us-ft +no_defs +ellps=clrk66 +nadgrids=@conus, at alaska, at ntv2_0.gsb, 
at ntv1_can.dat"
)
)

## locs.ref and locs.ll are created as above, the same as they were 3 years ago
## (same reproducible example data as then)

## use spTransform to go from WGS84 directly to the local system
## using spTransform()
locs.26743 <- spTransform(locs.ll, CRS("+init=epsg:26743"))

## distances relative to the reference transformation in August 2014
coordinates(locs.ref)-coordinates(locs.26743)
##  coords.x1 coords.x2
## [1,]  3.746539 -1.876668
## [2,]  3.746607 -1.876466

## distances relative to the reference transformation now
coordinates(locs.ref)-coordinates(locs.26743)
##  coords.x1 coords.x2
## [1,]  309.9325  20.05890
## [2,]  309.9136  20.03086

##
## a transformation that formerly was within a few feet of the example location 
is now some 300 feet away
##



## the ESRI "two step" starts with the lon/lat,
## Step 1 converts transforms them using what ESRI calls
## NAD_1983_To_WGS_1984_5   (wkid 1515)
## Step 2 transforms the resulting coordinates using
## NAD_1927_To_NAD1983_NADCON   (wkid 1241)


  Session info now:

library(rgdal)
Loading required package: rgdal
Loading required package: sp
rgdal: version: 1.2-13, (SVN revision 686)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
Path to GDAL shared files: /Users/macqueen1/Library/R/3.4/library/rgdal/gdal
Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
Path to PROJ.4 shared files: /Users/macqueen1/Library/R/3.4/library/rgdal/proj
Linking to sp version: 1.2-5

sessionInfo()
R version 3.4.2 (2017-09-28)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: OS X El Capitan 10.11.6

Matrix products: default
BLAS: 
/Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: 
/Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] openxlsx_4.0.17 Rkml_1.6-5  mymaps_1.4-2rmacq_1.3-5 
rgdal_1.2-13

[R-sig-Geo] Help with plotKML::kml() ?

[MacQueen, Don]

Can I get some help with the
   plotKML::kml()
function? I'm looking at various documentation and having trouble finding what 
I need.

Suppose "foo" is a SpatialPolygonsDataFrame object with a lat/long coordinate 
reference system. Then in R I can do, for example,

  plot(foo, col='yellow', border='red')

I want to write a KML file that will display the polygon in Google Earth with 
the same colors.

  plotKML::kml(foo, colour='yellow')

gets the fill, but I can't find how to set the border color.

I'd also like to be able to specify no fill, border only.

Thanks
-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062


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

[R-sig-Geo] writeOGR layer options for the KML driver?

[MacQueen, Don]

Hi,

Where can I find documentation for available layer_options (if any) for the 
writeOGR KML driver?

Specifically, I want to be able to specify the line color when writing a 
SpatialLines or SpatialLinesDataFrame object to a KML file for viewing in 
Google Earth.
(also specifyng a fill color and fill transparency would be nice, but I can 
certainly do without it)

I see that the plotKML package has a kml() function. I'll send a separate email 
asking for assistance with it.

Thanks
-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062


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

Re: [R-sig-Geo] Creating points along polylines

[MacQueen, Don]

sp::spsample can create equally spaced points along lines

I don't know if there's a function that will directly find the nearest points 
between two layers, but if necessary you can program it yourself, possibly 
starting with spatstat::crossdist to calculate the distances.

More generally, I would look into the rgeos and spatstat packages. And visit 
the CRAN spatial task view; that might lead you to additional packages.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062


On 2/10/17, 5:21 AM, "R-sig-Geo on behalf of Tristan Bourgeois" 
 
wrote:

Dear all,

I'm looking for a way to create points along polylines in 2 different cases
:

- First case :

I want to create points along  2 polyline layers (hydrographic network)
with a step between each point of *x *meters (as the QGIs pluggin "Locate
points along lines" is able to do)

Then, for these two "Points layers" I want to join the created points from
one layer to the nearest points of the other one (as the  QGIs pluggin
"MMGIS / HubDistance" does so)


- Second case :

I have a polyline layer (hydrographic network) which is sectionned into
several parts. I want to create points at each sections' extremities and
get teh coordinates of each of these points.

Hope my explanations are enough clear ?

Thanks in advance and have a sweet weekend !

Cheers.



-- 
Tristan Bourgeois

[[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] R - QGIS integration

[MacQueen, Don]

And if Kent is correct, then there is also:


subplot  package:Hmisc  R Documentation

Embed a new plot within an existing plot

Description:

 Subplot will embed a new plot within an existing plot at the
 coordinates specified (in user units of the existing plot).



(Which would not require learning ggplot2)
-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 11/22/16, 6:05 AM, "R-sig-Geo on behalf of Kent Johnson"
 wrote:

>If I understand correctly your question is not about QGIS but about using
>R
>to make maps with inset charts. There are several examples of charts on
>maps here:
>http://stackoverflow.com/questions/10368180/plotting-pie-graphs-on-map-in-
>ggplot
>
>Also ggmap + ggmap::inset looks like it will do what you want, see the
>ggplot2::annotation_custom example for a little help.
>
>Kent
>
>On Mon, Nov 21, 2016 at 6:00 AM,  wrote:
>
>
>> -
>>
>> Message: 7
>> Date: Mon, 21 Nov 2016 10:40:26 +0100
>> From: Domagoj Culinovic 
>> To: R-sig-Geo@r-project.org
>> Subject: [R-sig-Geo] R - QGIS integration
>> Message-ID:
>> > gmail.com>
>> Content-Type: text/plain; charset="UTF-8"
>>
>> Hi,
>> I have Polygon and point data in SHP format which i have use in QGIS.
>>also
>> i have several csv files where i have different statistical data about
>> points and polygons in SHP files.
>> This csv data are generated every week.
>> How to make R script which can make georeferenced charts, diagrams or
>> something like heathmap i R, and return picture or charts at theit
>>position
>> on the map?
>>
>> Regards,
>> Domagoj
>>
>> [[alternative HTML version deleted]]
>>
>>
>>
>> --
>>
>> Subject: Digest Footer
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>> --
>>
>> End of R-sig-Geo Digest, Vol 159, Issue 15
>> **
>>
>
>   [[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] DEM plot3D and overlay an aerial image

[MacQueen, Don]

Have you tried
  plotRGB(gruissan)
?

plotRGB is in the raster package.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 10/3/16, 2:38 AM, "R-sig-Geo on behalf of Mathieu Rajerison"

wrote:

>Hi R-List,
>
>
>I have a DEM on one hand, and on the other hand, I have an RGB aerial
>image
>
>I tried rasterVis and plot3D function, but I didn't find how to use the
>colors of my RGB aerial image.
>
>For the moment, I used rgl instead although it doesn't seem very
>appropriate for georeferenced data..
>
>Here's the RGL code :
>
>gruissan =
>stack("../MAISON/DATA/geo/raster/bdortho_marseille/marseille.tif")
>
>r = values(gruissan[[1]])/255
>g = values(gruissan[[2]])/255
>b = values(gruissan[[3]])/255
>col = rgb(r, g, b)
>
>library(rgl)
>persp3d(1:nrow(mnt), 1:ncol(mnt), values(mnt), col=col)
>
>Best,
>
>Mathieu
>
>
>ᐧ
>
>   [[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] create a triangular grid ?

[MacQueen, Don]

 spsample( type='hexagonal') will generate points in what some might call
a triangular grid.

They'll just be independent points, with no underlying structure declaring
them to be a grid, but at least you'll have them.

sp::over() may help you trim a rectangular grid; ?sp::over has an example
applying it to a grid.

-Don


-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 9/30/16, 4:26 AM, "R-sig-Geo on behalf of Fernando Paim"

wrote:

> 
>
>Hello everyone,
>
>I need to create a triangular grid for kriging.
>I've tried with makegrid, expand.grid but with just a rectangular grid.
>I have not been able to find a method to trim the rectangular grid to
>the newdata same shape of the original data, is this possible?
>
>The
>data (x,y,z) are available in
>http://www.paim.pro.br/downloads/dados.csv
>
>Can anyone help?
>
>Thanks in
>advance,
>
>Fernando Paim 
>
> 
>   [[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] offsetting overlapping points

[MacQueen, Don]

I don't know of a function for precisely what you're looking for, but in
the "roll your own" department, have you looked at the jitter() function?

You would probably have to run it separately on x and y, and tinker with
its arguments that control the offset. But it certainly should be doable.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 12/23/15, 8:44 AM, "R-sig-Geo on behalf of Paul Lantos"
 wrote:

>Hi,
>
> I have a point dataset in which some points overlap. I'd like to find a
>way to randomly offset these points within a certain tolerance so that
>none perfectly overlap. I can't find a way to do this within ArcGIS. I'd
>like to find a way to either impose a grid and randomize all points that
>fall in each grid cell, or specify a tolerance (i.e. offset each point
>randomly by a certain number of map distance units and not allow both x
>and y to be identical for any one point).
>
>
>Thanks!
>
>Paul
>
>
>__
>Paul M. Lantos, M.D., FIDSA, FAAP, FACP
>Pediatric Infectious Diseases
>Hospital Medicine
>Depts of Internal Medicine and Pediatrics
>Duke University Medical Center
>
>__
>
>This message and any included attachments are confidential and are
>intended only for the addressee(s).  The information contained herein may
>be confidential under the attorney/client privilege and/or the quality
>assurance and peer review privilege.  Unauthorized review, forwarding,
>printing, copying or otherwise distributing or using such information is
>prohibited and may be unlawful.  If you have received this message in
>error, please notify the sender promptly by email or telephone and delete
>the message.
>
>   [[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] spsample

[MacQueen, Don]

I'm not so deeply familiar with spsample that I know of a way to do this
directly.

However, I think it would be pretty easy to use spsample to first obtain
more points than needed, then calculate all the distances from the center,
then use the base::sample function to select a weighted sample, with
weights based on the distance from the center (see the prob argument to
sample).

-Don


-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 7/10/15, 2:59 AM, R-sig-Geo on behalf of Luca Candeloro
r-sig-geo-boun...@r-project.org on behalf of luca.candel...@gmail.com
wrote:

Hi,
following the example found in SpeciesDistributionModelling, a given
number
of points is drawn randomly within a circle:

x - circles(pts, d=5, lonlat=T, col='light gray')
# sample randomly from all circles
 samp1 - spsample(x@polygons, 250, type='random', iter=25)

Is there a way to specify weigths decreasing from the circle center, so
that random points are most likely near it?

thanks ,
Luca

   [[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] (no subject)

[MacQueen, Don]

This is a question for R-help, not R-sig-Geo.
Please also provide a subject, and do not send html.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 12/9/14, 7:17 AM, R-sig-Geo on behalf of rpaolo1967 RP
r-sig-geo-boun...@r-project.org on behalf of rpaolo1...@gmail.com wrote:

Hi everybody,



I apologize for this simple question. I¹m running a code which, at each
step of a loop, generates a data frame similar to this one:



d - rep(1:10, each = 6, len = 60)

a - rep(0:5, each = 1, len = 60)

v - rep(sort(runif(5, 30, 50), decreasing = T), each = 1, len = 60)

x - data.frame(t(cbind(d,a,v)))



I want to create and save plots of the results, once data frames are
created. I¹m using this code within the loop, after the creation of the
data frame:



# library(lattice)

svg(eraseMe.svg) # Format doesn¹t matter, and name is iteratively
created

xyplot(v ~ a | d, x, aspect = 5, type = o, xlab = a, ylab = v)

update(trellis.last.object(),

strip = strip.custom(strip.names = FALSE, strip.levels = i),

par.strip.text = list(cex = 1.1))

dev.off()



where ³i² is something as simple as



i - as.numeric(c(1:10))



It works troublelessly on my laptop (under 32 bits Windows Vista), but I
get the following message when running in computers under ubuntu or 64
bits
Windows 7:



Error in !lattice.getStatus(current.plot.saved, prefix = prefix) :

  invalid argument type



Notwithstanding, when I simply copy and paste again the following code,
after getting the message, it works fine (but loop is interrupted).



svg(eraseMe.svg) # No matter the format

xyplot(v ~ a | d, x, aspect = 5, type = o, xlab = a, ylab = v)

update(trellis.last.object(),

   strip = strip.custom(strip.names = FALSE, strip.levels = i),

   par.strip.text = list(cex = 1.1))

dev.off()



Any suggestion? Thanks a lot


Paolo

   [[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] error projecting with rgdal::project()

[MacQueen, Don]

If it helps (and since I could do it quickly and easily), I have success
with
  R 3.2.0
  gdal 1.11.2
  rgdal 0.9-3
  sp_1.1-0
  proj.4 4.9.1
(on a Mac; see below)

If I read your email correctly, that is a combination you haven't tried
(you have rgdal 0.9-2 with R 3.2.0).


 require(rgdal)
Loading required package: rgdal
Loading required package: sp
rgdal: version: 0.9-3, (SVN revision 530)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
 Path to GDAL shared files:
/Users/macqueen1/Library/R/3.2/library/rgdal/gdal
 Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
 Path to PROJ.4 shared files:
/Users/macqueen1/Library/R/3.2/library/rgdal/proj
 Linking to sp version: 1.1-0
 
 
 tst - cbind(c(-75.13574,-75.02346, -74.98747, -74.87270, -74.86581),
+ c(38.89733, 38.74623, 38.56410, 38.56395, 38.45108))
 
 xy - project(tst, +proj=tpeqd +lon_1=-80 +lat_1=33 +lon_2=-75
+lat_2=39)
 
 
 sessionInfo()
R version 3.2.0 (2015-04-16)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.9.5 (Mavericks)

locale:
[1] C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] rgdal_0.9-3 sp_1.1-0

loaded via a namespace (and not attached):
[1] grid_3.2.0  lattice_0.20-31




I don't know what autodetected in the path to proj.4 means, but that's
another difference between mine and yours.

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 6/8/15, 9:17 AM, Danielle Haulsee dhaul...@udel.edu wrote:

Hello,

I recently updated R and a bunch of packages and am now finding bits of
old
code that no longer run (of course!). Specifically I am having issues with
project() and spTransform(), likely both due to the same issue. It seems
to
perhaps be an issue with rgdal or perhaps one of its dependencies because
the problem does not appear related to the r-version.

I ran the following example code on two different machines, both running

R version 3.0.3 (2014-03-06) -- Warm Puppy, and either got an error or
no
error depending on the version of rgdal installed.

*This version was successful. *
require(rgdal)

Loading required package: rgdal
Loading required package: sp
rgdal: version: 0.9-1, (SVN revision 518)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08
Path to GDAL shared files:
/Library/Frameworks/R.framework/Versions/3.0/Resources/library/rgdal/gdal
Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
Path to PROJ.4 shared files:
/Library/Frameworks/R.framework/Versions/3.0/Resources/library/rgdal/proj
tst - cbind(c(-75.13574,-75.02346, -74.98747, -74.87270, -74.86581),
c(38.89733, 38.74623, 38.56410, 38.56395, 38.45108))
xy - project(tst, +proj=tpeqd +lon_1=-80 +lat_1=33 +lon_2=-75
+lat_2=39)

*This version gives me an error. *
require(rgdal)
Loading required package: rgdal
Loading required package: sp
rgdal: version: 0.8-16, (SVN revision 498)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08
Path to GDAL shared files: /usr/share/gdal
Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
Path to PROJ.4 shared files: (autodetected)

tst - cbind(c(-75.13574,-75.02346, -74.98747, -74.87270, -74.86581),
c(38.89733, 38.74623, 38.56410, 38.56395, 38.45108))
xy - project(tst, +proj=tpeqd +lon_1=-80 +lat_1=33 +lon_2=-75
+lat_2=39)
*Error in project(tst, +proj=tpeqd +lon_1=-80 +lat_1=33 +lon_2=-75
+lat_2=39) : *
*  major axis or radius = 0 or not given*

I also get the same error when running R version 3.2.0 (2015-04-16) --
Full of Ingredients and rgdal version:
Loading required package: rgdal
Loading required package: sp
rgdal: version: 0.9-2, (SVN revision 526)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
Path to GDAL shared files:
Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
Path to PROJ.4 shared files: (autodetected)

Does anyone have any ideas as to what is going on? I'm also getting the
same error (major axis or radius = 0 or not given), when using the
spTransform function. I've tried searching for this error and haven't come
up with much. Let me know if more information is needed. Thanks in advance
for any insight you might have.

  [image: photo]
*Danielle Haulsee*
PhD Candidate, University of Delaware
 m:(717)451-7636 | e:dhaul...@udel.edu | w:www.orb.ceoe.udel.edu | a:700
Pilottown Rd. Lewes, DE 19958
Get a signature like this:
https://ws-stats.appspot.com/r?rdata=eyJydXJsIjogImh0dHA6Ly93d3cud2lzZXN0
YW1wLmNvbS8/dXRtX3NvdXJjZT1leHRlbnNpb24mdXRtX21lZGl1bT1lbWFpbCZ1dG1fY2FtcG
FpZ249cHJvbW9fNDUiLCAiZSI6ICJwcm9tb180NV9jbGljayJ9
Click
here!

Re: [R-sig-Geo] Hole attribute in Polygon lost when made into Polygons object?

[MacQueen, Don]

Thanks, Roger, it does clarify.
-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 5/26/15, 10:18 PM, Roger Bivand roger.biv...@nhh.no wrote:

On Wed, 27 May 2015, MacQueen, Don wrote:

 I must be missing something basic
 (either in understanding, or possibly being blind to a typo!)

 Here is the example data:

 m1 -structure(c(2.8, 8, 8.2, 3.9, 1.6,
9.2, 6.8, 3, 1.1, 4.2),
   .Dim = c(5L, 2L))


 Pnh - Polygon(m1, hole=FALSE)
 Ph  -  Polygon(m1, hole=TRUE)

 Pnhs - Polygons(list(Pnh), ID=1)
 Phs - Polygons(list(Ph), ID=1)



 Then:
 Pnh@hole
 [1] FALSE
 Ph@hole
 [1] TRUE

 Pnhs@Polygons[[1]]@hole
 [1] FALSE
 Phs@Polygons[[1]]@hole
 [1] FALSE

 It appears that the hole attribute has been lost??


Yes, by design. A Polygons object is like an OGC/SFS MultiPolygon or
Polygon object, and must have at least one exterior ring (non-hole
sp:Polygon object). If the only Polygon object in a Polygons object has
its hole slot set to TRUE, this is treated as a misunderstanding and
silently changed to FALSE (and ring direction reversed if need be).
Please 
see the note in ?Polygons-class.

Hope this clarifies,

Roger


 sp_1.1-0
 R 3.1.2

 Thanks
 -Don



-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: roger.biv...@nhh.no


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

Re: [R-sig-Geo] overlap between shapefiles

[MacQueen, Don]

In the rgeos package there are  gIntersects() and gIntersection(), that
might be enough to get you started.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 5/8/15, 10:46 AM, Karla Shikev karlashi...@gmail.com wrote:

Dear all,

I'd like to take two shapefiles and to calculate the area of overlap based
on some world projection. Any suggestions about possible functions?

Thanks!

Karla

   [[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] Delimit a polygon for the region which is 1000 m from my raster altitude

[MacQueen, Don]

Look into the help page for the contour function in the raster package
It refers to
  RasterToContour

Follow the example given in ?RasterToContour
(copied here)
f - system.file(external/test.grd, package=raster)
r - raster(f)
x - rasterToContour(r, levels=500)
class(x)
plot(r)
plot(x, add=TRUE)


Then see:
class(x)

Then x can be written to a shapefile using writeOGR in the rgdal package.


-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 4/30/15, 3:12 PM, sadaoui sadaouimah...@outlook.com wrote:

*Hello,
I am new to this forum

My goal is to define a polygon for the region 1000 m from raster of
altitude
  then export to shapefile.

I tried with this code:

#import the raster: (attachment)

library (raster)
map - raster (Adige.tif)

#Delimit the area ( 1000 m) with the contour


contour (map, add = T, levels = c (1000))


but I can not export the contour to shapefile

thank you in advance for helping me

sadaoui*



--
View this message in context:
http://r-sig-geo.2731867.n2.nabble.com/Delimit-a-polygon-for-the-region-wh
ich-is-1000-m-from-my-raster-altitude-tp7588147.html
Sent from the R-sig-geo mailing list archive at Nabble.com.

___
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] Fitting a 3D anisotropic variogram

[MacQueen, Don]

There may be something in the RandomFields or georob packages.

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 4/13/15, 6:00 AM, Carlo Innocenti carlo.innoce...@isprambiente.it
wrote:

Hello,

anybody know a function to fit a theoretical 3D variogram model to the
experimental data?

I work with pollution data from boreholes samples of marine sediment.
Generally the model that best fits to the data is a nested variogram with
both geometric and zonal anisotropy, composed by:

- nugget effect
- spherical model with short range ( 1 m) along the z axis and very long
range (1000 km) along  x and y axis
- spherical model with very long range  (10^6 m) along the z axis and
medium range ( 1 km) along  x and y axis

In the past I used Isatis to visually adapt the range of the model to the
experimental variogram and automatically fit the sill.
Now I'm moving on open source software and I'm looking for a way to fit at
least the sill of the 3D variogram to the data.

I found that gstat allows to use the 3D variogram, but doesn't fit the
anisotropic ones, and geoR fit the anisotropic variogram but not the 3D
ones.

Anybody knows a way in R to build and fit a variogram as that described
above?

Thanks,
Carlo





-- 


Carlo Innocenti, PhD

ISPRA
Istituto Superiore per la Protezione e la Ricerca Ambientale
High Institute for Environmental Protection and Research

Via Vitaliano Brancati, 60
00144 Roma
Italy

phone: +39 06 5007 4639
email: carlo.innoce...@isprambiente.it

   [[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] inconsistent as.data.frame(SpatialPointsDF)

[MacQueen, Don]

In my experience, relying on column names to extract the coordinates is
not at all a good idea. I would strongly recommend that you take the time
to update all of your scripts to use the coordinates() function. I think
it will be worth it in the long run.

It's not a good idea because the column names of the coordinates depend on
how the SpatialPointsDataFrame was originally created, and in my own
applications that is highly variable.  Sometimes ('x','y'), sometimes
('lon','lat'), or any of several other variations of how to spell or
abbreviate latitude and longitude (with or without capitalization). Or
('easting','northing'). Or, or, or... Trying to carefully control all that
is more trouble than it's worth; I just use, for example,
coordinates(obj)[,1] and coordinates(obj)[,2] if I want to pull them out
as vectors. Ugly, but I can count on it.

That said, if
  as.data.frame(locs)
produces different names for the coordinates when used in different
contexts, then you've got something else going on that should not be going
on. This is where Frede's suggestions might help. You will need to
carefully track the construction of your locs object and see if it is
somehow different in the two situations. I don't know of any designed
circumstance that would explain this.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 3/19/15, 4:02 PM, dschneiderch dominik.schnei...@colorado.edu wrote:

I have a spatial points dF that is causing me trouble. I've figured out
what
is happening but without a clue why.

at the prompt, I do
 locs
class   : SpatialPointsDataFrame
features: 10
extent  : -112.0623, -109.0571, 33.65387, 36.32678  (xmin, xmax, ymin,
ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
variables   : 7
names   : network, State, Station_ID, Site_ID,Site_Name,
Elevation_ft, Elevation_m
min values  :SNTL,AZ, 09N05S, 308,  BAKER BUTTE,
  
7100,2164
max values  :SNTL,AZ, 12P01S,1143, HANNAGAN MEADOWS,
  
9200,2804

 head(as.data.frame(locs))
  xy network State Station_ID Site_ID   Site_Name
1 -111.4064 34.45660SNTLAZ 11R06S 308 BAKER BUTTE
2 -111.3827 34.45547SNTLAZ 11R07S1140 BAKER BUTTE SMT
3 -109.5034 33.97883SNTLAZ 09S01S 310   BALDY
4 -109.2166 33.69144SNTLAZ 09S06S 902 BEAVER HEAD
5 -109.0571 36.32678SNTLAZ 09N05S1143   BEAVER SPRING
6 -112.0623 35.26247SNTLAZ 12P01S1139   CHALENDER
  Elevation_ft Elevation_m
1 73002225
2 77002347
3 91252781
4 79902435
5 92002804
6 71002164

so as expected(?) my coordinate names get converted from Longitude,
Latitude
to x, y.

However, when I run my script, the output of head(as.data.frame(locs)) is:
Longitude   Latitude network State Station_ID Site_ID   Site_Name
1 -111.4064 34.45660SNTLAZ 11R06S 308 BAKER BUTTE
2 -111.3827 34.45547SNTLAZ 11R07S1140 BAKER BUTTE SMT
3 -109.5034 33.97883SNTLAZ 09S01S 310   BALDY
4 -109.2166 33.69144SNTLAZ 09S06S 902 BEAVER HEAD
5 -109.0571 36.32678SNTLAZ 09N05S1143   BEAVER SPRING
6 -112.0623 35.26247SNTLAZ 12P01S1139   CHALENDER
  Elevation_ft Elevation_m
1 73002225
2 77002347
3 91252781
4 79902435
5 92002804
6 71002164


I found out the hard way because i was doing
as.data.frame(locs)[,c('x','y')] to get the coordinates I switched
this
line to coordinates(locs) but I have other lines in the same code that use
as.data.frame() so I'm wondering if there are designed circumstance for
one
behavior compared to the other.  I did notice that
data.frame(locs)[,c('x','y')] seems to always maintain the original
coordinate names but I confirmed that the script uses as.data.frame()

Below is dput of a sample of my data. does anyone else get this behavior?

 dput(locs)
new(SpatialPointsDataFrame
, data = structure(list(network = c(SNTL, SNTL, SNTL, SNTL,
SNTL,
SNTL, SNTL, SNTL, SNTL, SNTL), State = c(AZ, AZ,
AZ, AZ, AZ, AZ, AZ, AZ, AZ, AZ), Station_ID = c(11R06S,
11R07S, 09S01S, 09S06S, 09N05S, 12P01S, 09S07S, 11P02S,
11P13S, 09S11S), Site_ID = c(308L, 1140L, 310L, 902L, 1143L,
1139L, 416L, 1121L, 488L, 511L), Site_Name = c(BAKER BUTTE,
BAKER BUTTE SMT, BALDY, BEAVER HEAD, BEAVER SPRING, CHALENDER,
CORONADO TRAIL, FORT VALLEY, FRY, HANNAGAN MEADOWS),
Elevation_ft = c(7300L, 7700L, 9125L, 7990L, 9200L, 7100L,
8400L, 7350L, 7200L, 9020L), Elevation_m = c(2225L, 2347L,
2781L, 2435L, 2804L, 2164L, 2560L, 2240L, 2195L, 2749L)), .Names =
c(network,
State, Station_ID, Site_ID, Site_Name, Elevation_ft,
Elevation_m), 

Re: [R-sig-Geo] inconsistent as.data.frame(SpatialPointsDF)

[MacQueen, Don]

Thanks, Edzer, this is helpful.

-Don


-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 3/20/15, 9:28 AM, Edzer Pebesma edzer.pebe...@uni-muenster.de wrote:

There is some logic: sp tries to track coordinate names when it can, but
if coordinates are set as a nameless matrix, as in the last example
below, it will choose names itself. coordnames() helps you discover how
the coordinates of a SpatialPointsDataFrame are called:

 df = data.frame(x=1:2, y=2:1, z = 3:4)
 df1 = df
 library(sp)
 coordinates(df1) = ~x+y
 as.data.frame(df1)
  x y z
1 1 2 3
2 2 1 4
 coordnames(df1)
[1] x y
 df1=df
 coordinates(df1) = df[1:2]
 coordnames(df1)
[1] x y
 df1=df
 coordinates(df1) = cbind(df$x,df$y) # nameless
 df1
  coordinates x y z
1  (1, 2) 1 2 3
2  (2, 1) 2 1 4
 coordnames(df1)
[1] coords.x1 coords.x2


a bit of a semantic trap is this: if your coordinates are longitude and
latitude and carry these names, after you project the object with
rgdal::spTransform they're still called longitude and latitude, although
they are no longer understood as such. You can then solve this yourself by

 coordnames(df1) = c(x, y)

after which

 coordnames(df1)
[1] x y


On 03/20/2015 05:01 PM, MacQueen, Don wrote:
 In my experience, relying on column names to extract the coordinates is
 not at all a good idea. I would strongly recommend that you take the
time
 to update all of your scripts to use the coordinates() function. I think
 it will be worth it in the long run.
 
 It's not a good idea because the column names of the coordinates depend
on
 how the SpatialPointsDataFrame was originally created, and in my own
 applications that is highly variable.  Sometimes ('x','y'), sometimes
 ('lon','lat'), or any of several other variations of how to spell or
 abbreviate latitude and longitude (with or without capitalization). Or
 ('easting','northing'). Or, or, or... Trying to carefully control all
that
 is more trouble than it's worth; I just use, for example,
 coordinates(obj)[,1] and coordinates(obj)[,2] if I want to pull them out
 as vectors. Ugly, but I can count on it.
 
 That said, if
   as.data.frame(locs)
 produces different names for the coordinates when used in different
 contexts, then you've got something else going on that should not be
going
 on. This is where Frede's suggestions might help. You will need to
 carefully track the construction of your locs object and see if it is
 somehow different in the two situations. I don't know of any designed
 circumstance that would explain this.
 
 -Don
 

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi),  University of Münster,
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers  Geosciences:   http://elsevier.com/locate/cageo/
Spatial Statistics Society http://www.spatialstatistics.info


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

Re: [R-sig-Geo] still having problems installing rgdal on osx mavericks

[MacQueen, Don]

You didn't say how you installed GDAL, nor did you show any error
messages, so it is difficult to say much. My guess would be that you have
some version incompatibility.

It is not the only way, but over the years, I have had success using
  http://www.kyngchaos.com/software/frameworks
and there is an rgdal there (and I do have a working rgdal on Mavericks).

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 2/28/15, 1:11 PM, Roger Bivand roger.biv...@nhh.no wrote:

On Sat, 28 Feb 2015, Cyrus Mohammadian wrote:

 Hey Everyone,

 I¹m having problems installing and loading rgal on mavericks. I know
 there is a lot written on the net about it but I have not been able to
 successfully install the package. I downloaded the .tgz file and I
 installed the package from source but every time I attempt to load the
 library my R workspace crashes. Any help would be fantastic. If

As you should know, files with the *.tgz extension are Mac Binary
packages. Those with *.tar.gz are source packages. Is this the cause of
your confusion?

Roger

 upgrading to Yosemite solves the problem I¹m willing to upgrade but
only 
 if I can be sure. Thanks everyone!

 Best,

 Cyrus Mohammadian
 ___
 R-sig-Geo mailing list
 R-sig-Geo@r-project.org
 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; fax +47 55 95 91 00
e-mail: roger.biv...@nhh.no

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

Re: [R-sig-Geo] Integrating a trajectory

[MacQueen, Don]

In two dimensions, interp() from the akima package would be a good
starting point. You might have to first generate a set of suitably spaced
lon,lat coordinates along the trajectory, depending on how the trajectory
is stored.

In the absence of something already existing for your 3d situation, I
suppose it could be done by first interpolating in 2d at altitudes above
and below the trajectory, then doing vertical interpolation between them.
I doubt that would be the most accurate interpolation, but it might be
good enough.

If you go to CRAN task view CRAN Task View: Handling and Analyzing
Spatio-Temporal Data, there is a whole section Moving objects,
Trajectories.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 2/26/15, 5:50 AM, Paul Woods p.wo...@qub.ac.uk wrote:

Hi everyone,

Perhaps this is a little removed from strict geography, but you seem to
be a clever bunch.

I¹m looking for a method in R or Python to integrate along a trajectory,
through a regular grid. So, imagine an aeroplane flying from London to
New YorkŠ I want to calculate how much material in the atmosphere the
aeroplane comes into contact with. I have a regular grid with the density
of the atmosphere at each lon, lat, altitude grid point, but of course
the plane does not travel along grid lines, so some interpolation would
be needed. Any ideas on where to start in calculating how much atmosphere
a plane would scoop up from London to New York?

Thanks,

Paul.
___
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] Problem with getGridIndex() ?

[MacQueen, Don]

I have defined an example SpatialPointsDataFrame, and then a GridTopology
following the example in Applied Spatial Data Analysis with R (2008, pg
48-49). With them, getGridIndex() fails. Suggestions would be much
appreciated.

Here's the (reproducible) example:


require(sp)

spd - data.frame(x=c(0,100), y=c(2,95), z=1:2)
coordinates(spd) - c('x','y')

cdim - 5
bb - bbox(spd)
cs - c(cdim, cdim)   ## cellsize
cc - bb[,1] + cs/2   ## the lower left cell center
cd - ceiling(diff(t(bb))/cs) ## cells to over all the data

gr - GridTopology(cellcentre.offset=cc,
   cellsize=cs,
   cells.dim=cd)

cnum - getGridIndex( coordinates(spd), gr)



Which results in:

 cnum - getGridIndex( coordinates(spd), gr)
   Min. 1st Qu.  MedianMean 3rd Qu.Max.
  0   5  10  10  15  20
Error in getGridIndex(coordinates(spd), gr) : this.idx out of range



This looks like a boundary issue, in that the upper and lower boundaries
of the grid cells in the x dimension are 0, 100, and so are the upper and
lower x coordinates.


I looked at getGridIndex, and found that changing
   outside = this.idx = g...@cells.dim[i] | this.idx  0
to
   outside = this.idx  g...@cells.dim[i] | this.idx  0

fixes the error in this example.

However, in my original example, which I will share upon request, it gets
past the out of range error, but then fails with the index outside
boundaries error.




R was started using
   R --vanilla

 sessionInfo()
R version 3.1.2 (2014-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)

locale:
[1] C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] sp_1.0-17

loaded via a namespace (and not attached):
[1] grid_3.1.2  lattice_0.20-29



-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062

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

Re: [R-sig-Geo] maptools issue

[MacQueen, Don]

The first thing to do is to find out the class of the object you are
plotting. For example, in this example:

library(maptools)
xx - readShapePoints(system.file(shapes/baltim.shp,
package=maptools)[1])
plot(xx)

Then:
  class(xx)
[1] SpatialPointsDataFrame
attr(,package)
[1] sp



This indicates that you are using the plot() function from the sp package.

 find('plot')
[1] package:sp   package:graphics


The plot function in the sp package is built on the plot function in the
graphics package, that is, the  plot() that is included with R. This
implies that the possible arguments include those of the basic plot
function. For example:

plot(xx, cex=2)
plot(xx, cex=2, col='red', pch=4)


Exactly what the additional arguments do depends on the class of the
object you're plotting.

For lines objects, see ?lines
For polygon objects, see ?polygons
For points objects, see ?points

(that list is an abbreviated version of Table 3.2 in the book Applied
Spatial Data Analysis with R, by Bivand, Pebesma, and Gomez-Rubio; get a
copy if you can!)

Final suggestion; if you are using the read functions from maptools as
in the example above, I would suggest installing the rgdal package and
using readOGR() instead. I think you will be better off in the long run
(and I'm pretty sure this is what the experts have been recommending).

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 12/25/14, 1:43 AM, kun17 kzh...@geo.ecnu.edu.cn wrote:

Hi all,

I want to use maptools to plot map, there are many plot examples in the
document, but I can not find the usage explaination for plot function as
well as the possible parameters.

Could someone give me advice?

Kun






--
View this message in context:
http://r-sig-geo.2731867.n2.nabble.com/maptools-issue-tp7587580.html
Sent from the R-sig-geo mailing list archive at Nabble.com.

___
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] readOGR cannot open .shp file

[MacQueen, Don]

What does
  list.files()
show ‹ after having set the working directory, of course.

Have you successfully loaded other shape files before, or is this the
first time?

Do any of the shape file examples in ?readOGR succeed?

Is the shape file complete? That is, are the .prj, .dbf, .shx, and so on,
files there?
(I¹m not sure what the complete required list is)

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 10/3/14, 12:15 PM, GraysonLee08 gray.vanderli...@gmail.com wrote:

I'm certain this is just user error however I can't seem to load a
specific
shape file into my current project. Any assistance would be greatly
appreciated. 

To start off, I'm trying to load the 2010 census tracts for the city of
Chicago so I can show population density. I've installed and loaded the
following packages:

 ## Install packages
 install.packages(rgdal)

 ## Load packages
 library(rgdal)
 library(ggplot2)
 library(rgeos)
 library(map tools)
 library(ggmap)


In the loading process for rgdal I notice that the version of rgdal that
was
loaded is 0.9-1 and the version of GDAL is 1.9.2. Additionally, Path to
GDAL
shared files is
/Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgdal/gdal

Loaded PROJ.4 runtime: Rel. 4.8.0
Path to PROJ.4 shared files:
/Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgdal/prom

I've set my working directory to the location of the .shp file
 getwd()
 [1] /Users/Gray/Dropbox/School/Food Desert

and use the following command to store my shape file:
 ChicagoCensusTract2010 - readOGR(dsn = ., ChicagoCensusTract)

I'm receiving the following error and am not really sure what to do next.
 Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding,
use_iconv = use_iconv,  :
   Cannot open file

I'd really appreciate any help I can get on this.



--
View this message in context:
http://r-sig-geo.2731867.n2.nabble.com/readOGR-cannot-open-shp-file-tp7587
228.html
Sent from the R-sig-geo mailing list archive at Nabble.com.

___
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] Identifying which points are in which cluster

[MacQueen, Don]

In the following reproducible example I have created a
SpatialPointsDataFrame with three clusters of points. What I¹m looking for
is a (good) way to add to the SPDF a column which identifies which cluster
each point is in.

‹‹ begin example ---

require(sp)
require(rgeos)


## construct at SpatialPointsDataFrame with three clusters
ctrs - cbind(x=c( 7000,  8000,  9000),
  y=c(12000, 13000, 14000))

pts - cbind(rep(ctrs[,1],3)+runif(9,-20,20),
 rep(ctrs[,2],3) + runif(9, -20,20))

plot(pts)

pts - SpatialPointsDataFrame(pts, data.frame(name=letters[1:9]) ,
proj4string=CRS('+init=epsg:26943'))

plot(pts)

## now pretend I don't which points are in which cluster
tmp1 - gBuffer(pts, width=30, byid=TRUE)
tmp2 - gUnaryUnion(tmp1)

## tmp2 is now a SpatialPolygons object
## with three Polygons, one for each cluster

plot(tmp2, usePolypath=FALSE)
plot(pts, add=TRUE)
points(ctrs, col='red', pch=3, cex=0.5)

‹ end example ‹

Now I need a way to identify which point is in which polygon, and add a
variable to the data frame slot of pts with that information.

I¹m sure I can work it out by digging into the structure of tmp2 and
pulling out the individual polygons using a loop, but I¹m hoping there¹s a
higher level solution, possibly using some combination of lapply() or
sapply() with over(). But I have not been able to come up with it.

Thanks
-Don


By the way, searching through old r-sig-geo emails, I found
  tmp.cc - hclust(dist(coordinates(pts)), complete)
  tmp.50 - cutree(tmp.cc, h=50)

and this works for this example (thanks to marcelino.delac...@upm.es), but
it¹s not clear to me which approach will be better in the long run for my
applications.

And I see something could be done with spdep:dnearneigh(), but the output
structure is a little complex and I don¹t understand it well enough.


So I would still appreciate suggestions for a solution based on points in
polygons and the data structures in the example.

Thanks very much
-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062

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

Re: [R-sig-Geo] Using spTransform() to reproduce another software package's transformation

[MacQueen, Don]

 +no_defs
 +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat
 +towgs=-3.746315,1.876856
 
 This corresponds to what you can find on www.epsg-registry.org
 
 Is there a corresponding EPSG code for the ESRI NAD_1983_To_WGS_1984_5
 transformation? Or can you by any other means find something similar
to the 
 CRS arguments? If so we can compare the CRS arguments from R and ESRI.
 
 
 
 
 Yours sincerely / Med venlig hilsen
 
 
 Frede Aakmann Tøgersen
 Specialist, M.Sc., Ph.D.
 Plant Performance  Modeling
 
 Technology  Service Solutions
 T +45 9730 5135
 M +45 2547 6050
 fr...@vestas.com
 http://www.vestas.com
 
 Company reg. name: Vestas Wind Systems A/S
 This e-mail is subject to our e-mail disclaimer statement.
 Please refer to www.vestas.com/legal/notice
 If you have received this e-mail in error please contact the sender.
 
 
 -Original Message-
 From: r-sig-geo-boun...@r-project.org [mailto:r-sig-geo-bounces@r-
 project.org] On Behalf Of MacQueen, Don
 Sent: 29. august 2014 01:58
 To: r-sig-geo@r-project.org
 Subject: [R-sig-Geo] Using spTransform() to reproduce another software
 package's transformation
 
 The program I work for has specified the use of a local coordinate
 reference
 system and a method for transforming and projecting from WGS84
long/lat
 to the local system. They use ESRI products to convert from long/lat
to 
 the
 local system.
 
 Since I do everything in R, naturally I wish to use spTransform() to
 replicate
 their conversion. I've been using spTransform() for a number of years
now,
 and thought I understood what I've been doing.
 
 But I have run into trouble. I would appreciate any advice.
 
 I believe I have a reproducible example. Toward the end of this email
are 
 R
 expressions (based on dput) that will create two SpatialPoints
objects 
 that
 are used in the example. They need to be created first, before
running the
 example.
 
 
 ## before adding further detail and the example, here are some
references
 
 
 (1)
 
http://downloads2.esri.com/support/TechArticles/Geographic_Transformati
 ons_10.1.zip
 (2)
 
http://resources.arcgis.com/en/help/main/10.2/index.html#/Equation_base
 d_methods/003r001200/
 (3)
 http://resources.arcgis.com/en/help/main/10.2/index.html#/Grid_based_m
 ethods/003r001300/
 
 
 
 The programs's specified CRS is epsg 26743 = California State Plane
Zone 3
 NAD27 US feet (out of my control!)
 
 The specified method for transforming and projecting from WGS84
long/lat
 to the local CRS consists of two steps:
  1) transform and project to epsg 2227 = California State Plane Zone
3 
 NAD83
 US feet
  2) transform to epsg 26743 = California State Plane Zone 3 NAD27 US
feet
 
 When doing the steps in the ESRI software's projection tool:
  step 1) use what ESRI calls NAD_1983_To_WGS_1984_5  (wkid 1515 in
 reference 1)
  step 2) use what ESRI calls NAD_1927_To_NAD_1983_NADCON  (wkid 1241
 in reference 1)
 
 According to reference 1 NAD_1983_To_WGS_1984_5 is a coordinate
 frame transformation.
 Based on reference 2, this means it is a 7 parameter Bursa-Wolf method
 Also based on reference 1, NAD_1927_To_NAD_1983_NADCON is a grid-
 based method
 
 
 As I hope you will see, a naive use of spTransform() produces
coordinates
 that differ from the ESRI two-step process by approximately 3.7 ft
(x) and 
 -
 1.9 ft (y). This is too large for our use case. I also believe as a
matter 
 of
 principle that it should be possible to do better (I'd like to
believe 
 that any
 transformation possible in ESRI is also possible using PROJ.4).
 
 
 ##
 ## reproducible example begins
 ##
 
 ## define two example points in WGS84 long/lat
 locs.xy - cbind(
  c(-121.524764291826,-121.523480804667),
  c(37.6600366036405,37.6543604613483)
  )
 
 locs.ll - SpatialPoints(locs.xy, proj4string=CRS(+proj=longlat
 +datum=WGS84) )
 
 ## source the expressions near the bottom of this email to create
 ##locs.ref
 ##locs.step1.esri
 
 ## use spTransform to go from WGS84 to the local system in one step:
 locs.26743 - spTransform(locs.ll, CRS(+init=epsg:26743))
 
 ## not close enough:
 coordinates(locs.ref)-coordinates(locs.26743)
 ##  coords.x1 coords.x2
 ## [1,]  3.746539 -1.876668
 ## [2,]  3.746607 -1.876466
 
 
 ## spTransform equivalent of ESRI step 1
 locs.step1.proj4 - spTransform(locs.ll, CRS(+init=epsg:2227))
 
 ## not close enough, essentially the same difference as above
 coordinates(locs.step1.esri)-coordinates(locs.step1.proj4)
 ##  coords.x1 coords.x2
 ## [1,]  3.746244 -1.877057
 ## [2,]  3.746315 -1.876856
 
 
 ## next, try the spTransform equivalent of ESRI step 1, but
specifying the
 seven parameters
 ## note, had to reverse the sign of the rotation args from wkid 1515
in
 reference 1;
 ## evidently the PROJ.4 default is the position vector method
(reference 
 2)
 
 crs.step1.cf - CRS('+proj=lcc +lat_1=38.43
 +lat_2=37.07 +lat_0=36.5 +lon_0=-120.5\
  +x_0=200.0 +y_0

[R-sig-Geo] Using spTransform() to reproduce another software package's transformation

[MacQueen, Don]

The program I work for has specified the use of a local coordinate reference 
system and a method for transforming and projecting from WGS84 long/lat to the 
local system. They use ESRI products to convert from long/lat to the local 
system.

Since I do everything in R, naturally I wish to use spTransform() to replicate 
their conversion. I’ve been using spTransform() for a number of years now, and 
thought I understood what I’ve been doing.

But I have run into trouble. I would appreciate any advice.

I believe I have a reproducible example. Toward the end of this email are R 
expressions (based on dput) that will create two SpatialPoints objects that are 
used in the example. They need to be created first, before running the example.


## before adding further detail and the example, here are some references


(1) 
http://downloads2.esri.com/support/TechArticles/Geographic_Transformations_10.1.zip
(2) 
http://resources.arcgis.com/en/help/main/10.2/index.html#/Equation_based_methods/003r001200/
(3) 
http://resources.arcgis.com/en/help/main/10.2/index.html#/Grid_based_methods/003r001300/



The programs’s specified CRS is epsg 26743 = California State Plane Zone 3 
NAD27 US feet (out of my control!)

The specified method for transforming and projecting from WGS84 long/lat to the 
local CRS consists of two steps:
 1) transform and project to epsg 2227 = California State Plane Zone 3 NAD83 US 
feet
 2) transform to epsg 26743 = California State Plane Zone 3 NAD27 US feet

When doing the steps in the ESRI software's projection tool:
 step 1) use what ESRI calls NAD_1983_To_WGS_1984_5  (wkid 1515 in reference 
1)
 step 2) use what ESRI calls NAD_1927_To_NAD_1983_NADCON  (wkid 1241 in 
reference 1)

According to reference 1 NAD_1983_To_WGS_1984_5 is a coordinate frame 
transformation.
Based on reference 2, this means it is a 7 parameter Bursa-Wolf method
Also based on reference 1, NAD_1927_To_NAD_1983_NADCON is a grid-based method


As I hope you will see, a naive use of spTransform() produces coordinates that 
differ from the ESRI two-step process by approximately 3.7 ft (x) and -1.9 ft 
(y). This is too large for our use case. I also believe as a matter of 
principle that it should be possible to do better (I’d like to believe that any 
transformation possible in ESRI is also possible using PROJ.4).


##
## reproducible example begins
##

## define two example points in WGS84 long/lat
locs.xy - cbind(
 c(-121.524764291826,-121.523480804667),
 c(37.6600366036405,37.6543604613483)
 )

locs.ll - SpatialPoints(locs.xy, proj4string=CRS(+proj=longlat +datum=WGS84) 
)

## source the expressions near the bottom of this email to create
##locs.ref
##locs.step1.esri

## use spTransform to go from WGS84 to the local system in one step:
locs.26743 - spTransform(locs.ll, CRS(+init=epsg:26743))

## not close enough:
coordinates(locs.ref)-coordinates(locs.26743)
##  coords.x1 coords.x2
## [1,]  3.746539 -1.876668
## [2,]  3.746607 -1.876466


## spTransform equivalent of ESRI step 1
locs.step1.proj4 - spTransform(locs.ll, CRS(+init=epsg:2227))

## not close enough, essentially the same difference as above
coordinates(locs.step1.esri)-coordinates(locs.step1.proj4)
##  coords.x1 coords.x2
## [1,]  3.746244 -1.877057
## [2,]  3.746315 -1.876856


## next, try the spTransform equivalent of ESRI step 1, but specifying the 
seven parameters
## note, had to reverse the sign of the rotation args from wkid 1515 in 
reference 1;
## evidently the PROJ.4 default is the position vector method (reference 2)

crs.step1.cf - CRS('+proj=lcc +lat_1=38.43 
+lat_2=37.07 +lat_0=36.5 +lon_0=-120.5\
 +x_0=200.0 +y_0=50.0 +ellps=GRS80 +datum=NAD83\
 +units=us-ft +no_defs\
 +towgs84=-0.991,1.9072,0.5129,0.025789908,0.0096501,0.0116599,0.0’)

## by the way, this alternative to specifying the CRS gives the same result
##  crs.step1.cf - CRS(+init=epsg:2227 
+towgs84=-0.991,1.9072,0.5129,0.025789908,0.0096501,0.0116599,0.0”)

locs.step1.cf - spTransform(locs.ll, crs.step1.cf)

## good enough (hooray!)
coordinates(locs.step1.esri)-coordinates(locs.step1.cf)
##  coords.x1 coords.x2
## [1,] -3.469177e-06 -5.122274e-08
## [2,] -3.418885e-06 -7.380731e-08



## now try for step 2 using spTranform()
locs.step2.cf - spTransform(locs.step1.cf, CRS(+init=epsg:26743))

## the original difference is back!
coordinates(locs.ref)-coordinates(locs.step2.cf)
##  coords.x1 coords.x2
## [1,]  3.746539 -1.876668
## [2,]  3.746608 -1.876466

## the implication is that in doing the transformation to epsg 26743, it 
reversed the effect of the 7-parameter method

## attempt to prevent that:

locs.tmp - locs.step1.cf
proj4string(locs.tmp) - CRS(+init=epsg:2227)
## Warning message:
## In `proj4string-`(`*tmp*`, value = S4 object of class CRS) :
##   A new CRS was assigned to an object with an existing CRS:
## +proj=lcc +lat_1=38.43 

Re: [R-sig-Geo] Randomly moving a locality (within set limits)

[MacQueen, Don]

My first thought would be to use spTransform() to convert to a (local)
cartesian coordinate system, do the shift [perhaps using
maptools::elide()], and then convert back.

Identifying appropriate (local) coordinate systems for each locality could
be a chore, but other than that it looks to me like a conceptually simple
process.

elide() will, I believe, remove the proj4string from the elided object so
you¹d have to put it back in afterwards. Maybe calculating on the
coordinates would be easier, but that would depend on what kind of Spatial
objects your localities are.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 8/22/14, 4:12 AM, Alastair Potts pott...@gmail.com wrote:

Hi all,

I was wondering if someone could help point me in the right direction here
- I can't seem to find a function or post that focuses on this.

I have localities around the world. I want to be able to randomly move a
given locality within a set radius (defined by km). So, I have a point at
xy and want it to be shifted to some other locality within, say, 15 km of
of its current locality.

This is simple enough using something like runif(1,-15,15), but it's the
lat-long conversion that is confusing me (how to work out how many decimal
degrees this might be around the point at different global localities).

Any help or pointers would be greatly appreciated.

Thanks in advance,

Cheers,
Alastair

   [[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] Dramatically slow plotting

[MacQueen, Don]

I'm not sure this will help, but try
   usePolypath=FALSE
in your call to spplot().

For further information, see
  help.search('polypath')
I get:

Help files with alias or concept or title matching 'polypath' using fuzzy
matching:
   graphics::polypath Path Drawing
   sp::SpatialPolygons-class  Class SpatialPolygons

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 3/17/14 3:23 PM, Rolando Valdez rvald...@gmail.com wrote:

Hello, 

Recently, I acquired a MacBook Pro, Core i7, 8 GB ram. I Installed the
newest R version, 3.0.3 from the web page. The problem is when I¹m
plotting maps, because is going very, very slow, about 3 or 4 minutes
just for a single map, while I¹ve done this in a few seconds in Windows
with Core i5 and 4 GB ram.

This is what I have:

R version 3.0.3 (2014-03-06) -- Warm Puppy
Copyright (C) 2014 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin10.8.0 (64-bit)

[R.app GUI 1.63 (6660) x86_64-apple-darwin10.8.0]

I found a reproducible example in web and I took time with proc.time()

ptm - proc.time()
library(sp)
library(lattice) # required for trellis.par.set():
trellis.par.set(sp.theme()) # sets color ramp to bpy.colors()

# prepare nc sids data set:
library(maptools)
nc - readShapePoly(system.file(shapes/sids.shp,
package=maptools)[1], proj4string=CRS(+proj=longlat +datum=NAD27))
arrow = list(SpatialPolygonsRescale, layout.north.arrow(),
offset = c(-76,34), scale = 0.5, which = 2)
#scale = list(SpatialPolygonsRescale, layout.scale.bar(),
#offset = c(-77.5,34), scale = 1, fill=c(transparent,black),
which = 2)
#text1 = list(sp.text, c(-77.5,34.15), 0, which = 2)
#text2 = list(sp.text, c(-76.5,34.15), 1 degree, which = 2)
## multi-panel plot with filled polygons: North Carolina SIDS
spplot(nc, c(SID74, SID79), names.attr = c(1974,1979),
colorkey=list(space=bottom), scales = list(draw = TRUE),
main = SIDS (sudden infant death syndrome) in North Carolina,
sp.layout = list(arrow), as.table = TRUE)

#sp.layout = list(arrow, scale, text1, text2), as.table = TRUE)
proc.time() - ptm

user  system elapsed
  2.408   0.064   2.616

It was quick.

Then I did a single plot with my shape:

mapa - readShapePoly(³Entidades_2013.shp²)
ptm - proc.time()
spplot(mapa[1]); proc.time() - ptm

user  system elapsed
 87.575   0.786  88.068

Why it take a lot of time? I worked with same shapes in Windows and never
took that time.

Hope you can help me,

Regards,

Rolando Valdez

___
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] build polygons from ellipses, then sample points within the polygons

[MacQueen, Don]

If you have the x,y coordinates of an ellipse, you can create a
SpatialPolygons or SpatialPolygonsDataFrame following the example in ?over
(package sp). 

For sampling, see ?spsample

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 3/10/14 12:19 PM, Anthony Fischbach afischb...@usgs.gov wrote:

I have estimated location error ellipses generated from a continuous-time
version of the correlated random walk model for animal telemetry data
(see http://cran.r-project.org/web/packages/crawl/index.html, functions:
crwMLE and crwPredict).
From these error ellipses, I wish to randomly select points.
Toward this end, I have two questions.

Question 1: Is there a handy function that readily allows creation of
polygons from ellipses?
## toy example
ellipses-data.frame(x=c(251753, 252966, 254179),
   y=c(343394, 343315, 343237),
   majorAxis=c(0.0, 2753.5, 
 2114.8),
   minorAxis=c(0.0,  2754.1, 
 2118.0),
   angle=c(0,0,0)) ## dataFrame 
 with Ellipse information
   ## in this case the majorAxis 
 is assumed to lie along the x axis
## x, y define the
centroid
of the ellipse.
coordinates(ellipses) - x~y  ## cast to a spatialPointsDataFrame
ellipses.poly -
Wonder.function.that.Makes.spatialPolygons.from.ellipses(ellipses)

Question 2:  I know that there is a handy means to randomly generate
locations from within a SpatialPolygon.  Please help me find it.



-
Tony Fischbach, Wildlife Biologist
Walrus Research Program
Alaska Science Center
U.S. Geological Survey
4210 University Drive
Anchorage, AK 99508-4650

afischb...@usgs.gov
http://alaska.usgs.gov/science/biology/walrus
--
View this message in context:
http://r-sig-geo.2731867.n2.nabble.com/build-polygons-from-ellipses-then-s
ample-points-within-the-polygons-tp7585922.html
Sent from the R-sig-geo mailing list archive at Nabble.com.

___
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] Represent data on map

[MacQueen, Don]

Cannot help much without more information, but try:

  names(map)

Then pick one of the column names. Apparently, 2000 is not the name of a
column of data in the map object.
If map is like unem, then perhaps

  spplot(map, 'X2000')

is what you need.

Note, the c() in
  c(2000)
is not necessary. Simply
  2000
would be sufficient.

I am not familiar with the readShapeSpatial function; I always use
readOGR() from the sp package.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 1/3/14 1:00 PM, Rolando Valdez rvald...@gmail.com wrote:

Hi,

I'm trying to represent data (unemployment rate) on a map.

This is what I have:

 library(sp)
 library(maptools)
Checking rgeos availability: TRUE
 map - readShapeSpatial(ESTADOS)
 read.csv(C:\\Users\\Rolando\\Documents\\Maestría en Economía\\3er
Semestre\\AEMT\\desempleo_est.csv)
 spplot(map, c(2000))
Error en `[.data.frame`(obj@data, zcol) : undefined columns selected

 unem - read.csv(C:\\Users\\Rolando\\Documents\\Maestría en
Economía\\3er Semestre\\AEMT\\desempleo_est.csv)
 names(unem)
[1] num_est est X2000   X2005   X2010   X2013

I hope you can help me, this is my first time in R.

Regards,

Rol~

___
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] making an inset in a map figure

[MacQueen, Don]

Perhaps use subplot() from the Hmisc package

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 9/10/13 9:00 AM, Waichler, Scott R scott.waich...@pnnl.gov wrote:

Hi, Does anyone have an example of making a map figure with an inset,
where a portion of the map in the large figure is shown at a larger scale
in the inset?

Scott Waichler
Hydrology Group
Pacific Northwest National Laboratory
Richland, WA, USA

___
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] Global mapping of multiple categories by colour

[MacQueen, Don]

Rhona,

I am completely unfamiliar with the package and functions you're using,
so I had supposed you would use functions from the sp and rgdal packages.
For what it's worth, I'll give an example below.

However, probably you will need to share at least a representative portion
of your data in order for people to give useful help. One way is to
include the ouput of the dput() function. In your case, probably,
  dput(data3)
or perhaps
  dput(head(data3))
if data3 is big.


Here is an example that requires that your countries be stored in a
SpatialPolygonsDataFrame. In this example, first we create some arbitrary
polygons, then plot them. It's not complete solution for your situation,
but is, I hope, analogous.

Probably, my email software will mangle the line lengths.

## example data, creating a SpatialPolygonsDataFrame
## copied from help('over',package='sp')
r1 = cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409,
 c(332349, 332057, 332342, 333250, 333558, 333676,   180162, 180114),
 332618, 332413, 332349))
r2 = cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437,
  179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683,
 331133, 331623, 332152, 332357, 332373))
r3 = cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875,
  179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
  c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004,
329783, 329665, 329720, 329933, 330478, 331062, 331086))
r4 = cbind(c(180304, 180403,179632,179420,180304),
  c(332791, 333204, 333635, 333058, 332791))
 
sr1=Polygons(list(Polygon(r1)),r1)
sr2=Polygons(list(Polygon(r2)),r2)
sr3=Polygons(list(Polygon(r3)),r3)
sr4=Polygons(list(Polygon(r4)),r4)
sr=SpatialPolygons(list(sr1,sr2,sr3,sr4))
srdf=SpatialPolygonsDataFrame(sr, data.frame(cbind(1:4,5:2),
row.names=c(r1,r2,r3,r4)))


## now plot
plot(srdf)
plot(srdf,col=1:4)




-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 8/28/13 8:58 AM, Rhona Govender rho...@gmail.com wrote:

Dear All,

Sorry for cross posting, I have been kindly directed here by the r help
group.

I have to make a last minute global map. I'm trying to do this in r,
but I am a beginner.

I have a global dataset of dependent values (eg. rate of common cold)
for 100/240 countries/offshore territories. I have a predictor
variable (eg. Human development index or HDI). What I would like to do
is colour a global map by HDI grouping (low [0-0.5], medium[0.51-0.7],
high[0.71-1.0]). So for countries that I have data for, if it is low
HDI it would be colour1, medium colour 2, low colour3. For countries
with no data they will be a fourth colour (grey).

Additionally, some small island nations will be too small to be
visualized. Is there a way to manually increase the size of the
colouration to be visible on a map at the global scale?

My data is as such (where blanks are no data):

Country HDI Colds Albania 0.42 .72 Austria 0.89 China 0 .76 .12

I have found this example to try and build off:

http://stackoverflow.com/questions/11225343/how-to-create-a-world-map-in-r
-with-specific-countries-filled-in?rq=1

To change this example, I imported my csv, picked out the columns I
wanted (country ISO3 codes [country] and the hdi predictor[hdi]) and
named them(country=data3$country) then and made it into a data frame
(data3).

I have this code I am confused with:

malMap - joinCountryData2Map(dF = data3, joinCode = ISO3,
nameJoinColumn = country)

# This will join your malDF data.frame to the country map data

mapCountryData(malMap, nameColumnToPlot=hdi, catMethod =
categorical, missingCountryCol = gray(.8))

head(data3)

country   hdi
  0
  0
  0
  ASM 0.827
  0
  0

Firstly I get this error as it won't match the country names up:

Error in joinCountryData2Map(dF = data3, joinCode = ISO3,
nameJoinColumn = country) :  your chosen nameJoinColumn :'country'
seems not to exist in your data, columns =

How can I change this sample code to colour the three levels of hdi?
Is it possible with this code? Also, is it possible to manually make
the small island nations bigger as they will not appear at this scale.

I'm sorry this is so long. Any help would be greatly appreciated.

Cheers,

Sarah

___
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 wirh R 3.0.1 and maptools

[MacQueen, Don]

What happens if you try these, derived from the examples in ?rgdal:

## from ?readOGR
dsn - system.file(vectors, package = rgdal)[1]

## then read two of the shapefiles in that directory
cit - readOGR(dsn,'cities')
plot(cit)

tmp - readOGR(dsn,'trin_inca_pl03')
plot(tmp)


They work for me, also having recently upgraded to R 3.0.1 on a Mac.

I see that error message from time to time, and if I recall correctly it
usually has something to do with which packages I have (or don't have)
attached. 

Similar to what Lee said, I would make sure you have up to date versions
of everything that maptools and rgdal depend on.

Compare with the versions I have:
 require(rgdal)
Loading required package: rgdal
Loading required package: sp
rgdal: version: 0.8-10, (SVN revision 478)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08
Path to GDAL shared files:
/Library/Frameworks/R.framework/Versions/3.0/Resources/library/rgdal/gdal
Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
Path to PROJ.4 shared files:
/Library/Frameworks/R.framework/Versions/3.0/Resources/library/rgdal/proj



I'd also suggest including the the output of sessionInfo() in your posts
(standard practice for all requests for R help)

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 8/17/13 1:20 PM, sara martino saramart...@yahoo.com wrote:

Hi
i am new to using shapefiles in R. I just want to plot some maps and I
used to do it as: 


shp=readShapeSpatial(file.shp)
plot(shp)

using maptools or 


shp -readOGR(path/to/shpfiles,IND_adm2)plot(shp)

using rgdal.

I have now upgraded to R 3.0.1 and since then everytime i try to plot the
shapefile i get the error:

Error in as.double(y) :
cannot coerce type 'S4' to vector of type 'double'

I get the same error also if i try to reproduce the example in
?readShapeSpatial.

Can anyone tell me what is the problem?
thanks
sara

   [[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] Error in plotKML

[MacQueen, Don]

Apparently, because you don't have gnome-open installed on your system.
Why it wants gnome-open is a mystery.
But this is probably how it is trying to run GoogleEarth.
The equivalent command in OS X is simply 'open'

But, if you open a Finder window to your working directory [type
   getwd()
in your R session] and look for a file named
   eberg__CLYMHT_A__.kml
and double-click on it, it should open in GoogleEarth.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 8/2/13 4:09 PM, Manuel Spínola mspinol...@gmail.com wrote:

Dear list members,

Working on X OS mountain lion and R 3.0.1

I am running the example from the plotKML tutorial (
http://gsif.isric.org/doku.php?id=wiki:tutorial_plotkml) and I get this
error:

 data(eberg)
 coordinates(eberg) - ~X+Y
 proj4string(eberg) - CRS(+init=epsg:31467)
 eberg - eberg[runif(nrow(eberg)).2,]
 bubble(eberg[CLYMHT_A])
 plotKML(eberg[CLYMHT_A])
Plotting the first variable on the list
KML file opened for writing...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Closing  eberg__CLYMHT_A__.kml
sh: gnome-open: command not found

Why I get this error?

Best,

Manuel

-- 
*Manuel Spínola, Ph.D.*
Instituto Internacional en Conservación y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
mspin...@una.ac.cr
mspinol...@gmail.com
Teléfono: (506) 2277-3598
Fax: (506) 2237-7036
Personal website: Lobito de río
https://sites.google.com/site/lobitoderio/
Institutional website: ICOMVIS http://www.icomvis.una.ac.cr/

   [[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] projection issue - can one rotate a projection?

[MacQueen, Don]

The elide function in the maptools package has a 'rotate' argument and may
be able to do what you need.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 7/31/13 3:17 AM, Tanya Rijkenberg Compton
tanyajcomp...@yahoo.com.au wrote:

Hi all

My colleague, a physical oceanographer, and I are having problems getting
his geographic layers into the correct projection (after importing as a
netcdf file), and we are hoping someone can help us. We are using the
Raster package.

My colleague obtains a Dutch Rijksdriehoek in the first instance but then
rotates this by 17 degrees to run his models. He then provides me with a
netcdf file where the layers are still projected at 17 degrees from the
normal Rijksdriehoek.

Is there a way in R to rotate the layers back from the 17 degrees?

If anyone can help we would greatly appreciate it!

Thanks

Tanya and Matias



R Code 
# Load the packages:
library(sp)
library(ncdf4)
library(raster)

# get the netcdf and read in to R
d - raster(test_layers.nc, varname = time_exp)

# rijksdriehoek projection
RD -CRS(+proj=sterea +lat_0=52.156160 +lon_0=5.387639
+k=0.079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs)
##http://spatialreference.org/ref/epsg/28992/proj4/,

# reproject to Rijksdriehoek
pr1 - projectRaster(d, crs=newproj)
   [[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] Linux or OS X on Macbook Pro Retina?

[MacQueen, Don]

I have used all of these except spgrass6 on Mac OS X for quite a few years
now.

I use the kyngchaos framework builds of the supporting libraries (proj4,
gdal, an others), and also his builds of QGIS and GRASS.
I have never used macports or homebrew for these (I do use macports to get
an X11-aware version of emacs).

In some circles, people advise against using macports, I guess mostly to
avoid confusion with their versions of things, and the versions that come
with the OS. I neither agree nor disagree with that position. I suppose
that reasoning might apply to homebrew. I'm more inclined to say, if I
don't need those systems I don't use them, and I haven't needed them for
GIS work in Mac OS.

The caveat is that I haven't updated to R 3.x.x yet, nor OS to OS X 10.8

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 6/19/13 1:13 AM, Rainer M Krug rai...@krugs.de wrote:

Hi

I finally received my new computer - a Mac Powerbook Retina.

I am a Linux / Ubuntu user, but I am thinking about possibly using OS
X. How does OS X play with these spatial applicatios, mainly R (and
associated packages), gdal, proj4, GRASS? Are there problems /
complications with using spgrass6 under OS X?

Any experiences, things to be wary about?

Thanks,

Rainer

-- 
Rainer M. Krug

email: RMKrugatgmaildotcom

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

[R-sig-Geo] writeRaster and RGB vs Gray

[MacQueen, Don]

Here's a little example script using functions from the raster package,

r1 - r2 - r3 - raster(ncol=25,nrow=25,
 xmn=0, ymn=0, xmx=1, ymx=1)
r1[] - round(runif(25^2,0,255))
r2[] - round(runif(25^2,0,255))
r3[] - round(runif(25^2,0,255))
s - stack(r1,r2,r3)

plotRGB(s)

writeRaster(s,'s.tif', format='GTiff', overwrite=TRUE,
bandorder='BSQ')


When I open the file s.tif in QGIS, it's fine; I see the same image as in
R.

But the Mac platform's default viewer, Preview.app, thinks the file's
color model is Gray, rather than RGB, and it does not render properly.
Preview does recognize, for example, a jpeg photo as having an RGB color
model.
(Adobe Reader (for Mac) wouldn't open it at all.)

Is there some better set of arguments to give writeRaster that will
somehow make the file more broadly recognizable as RGB?

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062

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

Re: [R-sig-Geo] make polygons transparent in plotKML

[MacQueen, Don]

Just a pointer to hopefully get you started:

see ?col2rgb and look at the alpha parameter.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 4/18/13 9:03 AM, Seth Bigelow s...@swbigelow.net wrote:

Greetings: 

 

Could anyone advise me how to make transparent polygons in plotKML? I am
overlaying polygons in a SpatialPolygonsDataFrame object into Google
Earth -
Seems I can make the polygons any color (e.g., yellow - see example)
except
transparent

 

plotKML(mySPDF, colour_scale =#00)

 

Appreciately, Seth Bigelow


   [[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] spplot

[MacQueen, Don]

To use spplot:
  Your data must have coordinates.
  The coordinates do not have to be longitude and latitude.

You have to tell R where to find the coordinates, it cannot find them for
you.

Telling R where to find the coordinates might be very easy, depending on
what is already in your data, or where it came from and how you loaded it
into R.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 2/15/13 6:06 AM, Milan Sharma milansharma2...@yahoo.com wrote:

Dear list/ Dr. Bivand,
is it possible to plot spplot if there is no Longitude and Latitude in my
data? (i.e. asking R to find latitude and longitude itself, is there any
way to do that)?
Milan

___
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] converting latitude and longitude to UTM

[MacQueen, Don]

To get started, try the following:

require(sp)

bbox(NYC2)
bbox(NYBorough)
range(myfile$X)
range(myfile$Y)

and see if the ranges of X,Y are even close to the coordinate ranges of
your shapefiles.

Then try

coordinates(myfile) - c('X','Y')
bbox(myfile)
plot(myfile)

If the problem is what you suspect, i.e., lat/long vs. UTM (which seems
likely), then you will have to start learning about how to specify
coordinate systems in R. It's a much bigger topic than I can go into here.
See the help for
  ?proj4string
  ?spTransform   (as suggested by Mike Sumner in another response)


You will need the rgdal package for spTransform.

For what it's worth, I use Quantum GIS (QGIS) when I need to find out
information about a particular projection.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 1/31/13 5:03 AM, Abby Rudolph arudo...@pire.org wrote:

Hello,

I am using the following code to create ppp files from csv data and map
shape files, but I am getting some errors which I have been unable to fix
by searching them online:

library(spatstat)
library(maps)
library(maptools)

NYC2-readShapePoly(nybb.shp) # this is a map of the NYC boroughs
without waterways and no census tract divisions (but it does include
lines separating the 5 boroughs)
plot(NYC2)

NYBorough-readShapePoly(NYBoroughsShapeMerge.shp) # this is a map of
the NYC boroughs with census tract divisions
plot(NYBorough)

myfile-read.csv(myfile.csv)
# my csv file contains 4 columns (ID, variable_type, X, Y)
# my goal is to determine whether individuals of variable_type=1 are
spatially distributed more than what would be expected due to chance, if
individuals of variable_type=0 are spatially distributed more than what
would be expected due to chance, and finally if the spatial distribution
of variable_type=1 differs from that of variable_type=0.

x-NYBorough@polygons[[1]]@Polygons
# ppp(x.coordinates, y.coordinates, x.range, y.range)
coords-slot(x[[1]],coords)
coords
coords-coords[19:1,]
coords-coords[-1,]
coords

u-myfile$type==0

Type_0-ppp(myfile$X[u],myfile$Y[u],
window=owin(poly=list(x=coords[,1],y=coords[,2])))
Type_1-ppp(myfile$X[!u],myfile$Y[!u],
window=owin(poly=list(x=coords[,1],y=coords[,2])))

ERROR Message:
Warning message:
In ppp(myfile$X[u], myfile$Y[u], window = owin(poly = list(x = coords[,  :
  217 points were rejected as lying outside the specified window
Warning message:
In ppp(myfile$X[!u], myfile$Y[!u], window = owin(poly = list(x = coords[,
 :
  435 points were rejected as lying outside the specified window

# I only have 652 points, so all of my points are rejected as lying
outside the specified window


Below is the sample code I was using as a template:

x-BaltCity2@polygons[[1]]@Polygons
coords-slot(x[[1]],coords)
coords-coords[167:1,]
coords-coords[-1,]
coords.mi-coords*0.000621371192
# convert meters to miles

play2-ppp(play$X_m*0.000621371192,play$Y_m*0.000621371192,
   window=owin(poly=list(x=coords.mi[,1],y=coords.mi[,2])))

There are two things I can think of that might be different between my
data set and the data that is used for this template code - and might
cause the warning:
1) My shapefile is for NYC and includes 5 different boroughs, whereas the
shapefile for Baltimore looks more like a regular polygon (like a
rectangle).
2) My data has only lat/lon coordinates and the data corresponding with
the template code I am working off of is using what I think is UTM (X_m
and Y_m).

Is there another way I can run the code above that will utilize the
latitude and longitude coordinates and will not return a warning?

Thanks,
Abby

From: Abby Rudolph
Sent: Wednesday, January 30, 2013 2:54 PM
To: r-sig-geo@r-project.org
Subject: converting latitude and longitude to UTM

Hello,
Is there a formula to convert x-/y-coordinates (latitude /longitude) into
UTM.  My data is in the form of latitude and longitude, but the code I
wish to use relies on x_m and y_m which I believe are UTM measurements (x
in meters and y in meters).
Thanks
Abby

___
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] lapply and SpatialGridDataFrame

[MacQueen, Don]

, SpatialPolygonsDataFrame))
 user system elapsed
 2.875 0.043 2.945

If the SpatialGridDataFrame objects all share the same GridTopology,
you 
should make a single SpatialGridDataFrame in a loop, or a single
RasterBrick or RasterStack in the raster package, and go from there.
But 
your proposed workflow is so opaque that helping you fix lapply isn't
going to move you forward at all.

You must think through your workflow carefully - it is very possible
that 
you can create the BayesX bnd object directly from the GridTopology of
the 
input data, without any intermediate SpatialPolygons objects or
shapefiles. The authors of BayesX would have done everyone a big favour
if 
they had supported sp and spacetime classes, for example by coercion -
spatstat is a good example of how to use both sp and package-specific
classes, and possibly used facilities in sp and spacetime for
visualisation.

Roger

 I have tried to update the previous R code to save each output file
to a
 folder:

 library(maptools)
 library(Grid2Polygons)

 readfunct - function(x)
 {
 u - readAsciiGrid(x)
 }

 modfilesmore - paste0(MaxFloodDepth_, 1:54, .txt)
 modeldepthsmore - lapply(modfilesmore, readfunct)

 myfolder - paste0(/MaxFloodDepthImages/MaxFloodDepthPolygon_,
1:54)
 maxdepth.plys - lapply(modeldepthsmore, Grid2Polygons, level=FALSE,
 cat, file = myfolder, append = TRUE)

 Error in FUN(X[[1L]], ...) :
 unused argument(s) (file = myfolder, append = TRUE)


 -Original Message-
 From: MacQueen, Don [macque...@llnl.gov]
 Sent: 1/28/2013 12:34:56 PM
 To: iruc...@mail2world.com;r-sig-geo@r-project.org
 Subject: Re: [R-sig-Geo] lapply and SpatialGridDataFrame

 I suspect you're not passing arguments correctly to lapply()

 I'd try
 maxdepth.plys - lapply(modeldepthsmore, Grid2Polygons, level=FALSE)

 -Don



 --
 Don MacQueen

 Lawrence Livermore National Laboratory
 7000 East Ave., L-627
 Livermore, CA 94550
 925-423-1062





 On 1/27/13 2:47 AM, Irucka Embry iruc...@mail2world.com wrote:

 Hi all, I have a set of 54 files that I need to convert from ASCII
 grid
 format to .shp files to .bnd files for BayesX.

 I have the following R code to operate on those files:

 library(maptools)
 library(Grid2Polygons)

 readfunct - function(x)
 {
 u - readAsciiGrid(x)
 }

 modfilesmore - paste0(MaxFloodDepth_, 1:54, .txt)
 modeldepthsmore - lapply(modfilesmore, readfunct)

 maxdepth.plys - lapply(modeldepthsmore,
 Grid2Polygons(modeldepthsmore,
 level = FALSE))

 layers - paste0(examples/floodlayers_, 1:54)
 writePolyShape(maxdepth.plys, layers)
 shpName - sub(pattern=(.*)\\.dbf, replacement=\\1,
 x=system.file(examples/Flood/layer_.dbf, package=BayesX))
 floodmaps - shp2bnd(shpname=shpName, regionnames=SP_ID)

 ## draw the map
 drawmap(map=floodmaps)


 This is the error message that I receive:
 maxdepth.plys - lapply(modeldepthsmore,
 Grid2Polygons(modeldepthsmore, level = FALSE))
 Error in Grid2Polygons(modeldepthsmore, level = FALSE) :
 Grid object not of class SpatialGridDataFrame


 Can someone assist me in modifying the R code so that I can convert
 the
 set of files to .shp files?

 Thank-you.

 Irucka Embry

 ___
 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


-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: roger.biv...@nhh.no


span id=m2wTlpfont face=Arial, Helvetica, sans-serif size=2
style=font-size:13.5px_
__BRGet the Free email that has everyone talking at a
href=http://www.mail2world.com
target=newhttp://www.mail2world.com/abr  font
color=#99Unlimited Email Storage #150; POP3 #150; Calendar #150;
SMS #150; Translator #150; Much More!/font/font/span
   [[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] lapply and SpatialGridDataFrame

[MacQueen, Don]

I suspect you're not passing arguments correctly to lapply()

I'd try
 maxdepth.plys - lapply(modeldepthsmore, Grid2Polygons, level=FALSE)

-Don



-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 1/27/13 2:47 AM, Irucka Embry iruc...@mail2world.com wrote:

Hi all, I have a set of 54 files that I need to convert from ASCII grid
format to .shp files to .bnd files for BayesX.

I have the following R code to operate on those files:

library(maptools)
library(Grid2Polygons)

readfunct - function(x)
{
u - readAsciiGrid(x)
}

modfilesmore - paste0(MaxFloodDepth_, 1:54, .txt)
modeldepthsmore - lapply(modfilesmore, readfunct)

maxdepth.plys - lapply(modeldepthsmore, Grid2Polygons(modeldepthsmore,
level = FALSE)) 

layers - paste0(examples/floodlayers_, 1:54)
writePolyShape(maxdepth.plys, layers)
shpName - sub(pattern=(.*)\\.dbf, replacement=\\1,
x=system.file(examples/Flood/layer_.dbf, package=BayesX))
floodmaps - shp2bnd(shpname=shpName, regionnames=SP_ID)

## draw the map
drawmap(map=floodmaps)


This is the error message that I receive:
 maxdepth.plys - lapply(modeldepthsmore,
Grid2Polygons(modeldepthsmore, level = FALSE))
Error in Grid2Polygons(modeldepthsmore, level = FALSE) :
Grid object not of class SpatialGridDataFrame


Can someone assist me in modifying the R code so that I can convert the
set of files to .shp files?

Thank-you.

Irucka Embry 


span id=m2wTlpfont face=Arial, Helvetica, sans-serif size=2
style=font-size:13.5px_
__BRGet the Free email that has everyone talking at a
href=http://www.mail2world.com
target=newhttp://www.mail2world.com/abr  font
color=#99Unlimited Email Storage #150; POP3 #150; Calendar #150;
SMS #150; Translator #150; Much More!/font/font/span
   [[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] readOGR or writeOGR not closing connections?

[MacQueen, Don]

Roger,

My earlier reply to your request to test a newer version of rgdal seems to
have not been sent.

I can't receive attachments that are zip files. Could you change the
suffix to, say .doc and resend?

For everyone else on r-sig-geo, I apologize for the email noise; this is
an attempt to work around another email problem.

Thanks
-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062

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

[R-sig-Geo] readOGR or writeOGR not closing connections?

[MacQueen, Don]

I am using R to help me interactively create and manage GIS layers, stored
as shapefiles, and I'm encountering a problem. From time to time I get an
error:

Error in writeOGR(tmpp, shps, fn, a.name, driver = ESRI Shapefile, :
Layer creation failed

I can work around it by quitting and restarting R, but it would be nice if
I didn't have to.


Below is a little script that reproduces the behavior (on my machine).
('shps' is a sub-directory)



The implication, as I see it, is that one or both of readOGR or writeOGR
is not closing some connections. If I'm right about this, it would be nice
to have a way to force-close them whenever necessary, but I haven't found
a way to do that. Information for how, or any other suggestions, would be
appreciated.



Note that in each iteration there are 4 files written and 4 files read.
  2*4*14=112
which is getting close to the number available according to ?open:

-- quote --
A maximum of 128 connections can be allocated (not necessarily
 open) at any one time.  Three of these are pre-allocated (see
 'stdout').  The OS will impose limits on the numbers of
 connections of various types, but these are usually larger than
 125.
-- end quote --





 begin script
require(sp)
require(rgdal)

print(sessionInfo())

data(meuse.grid)
coordinates(meuse.grid) = ~x+y
tmpp - meuse.grid[1:10,]

for(i in 1:100) {
  fn - paste('tmpp',i,sep='')
  writeOGR(tmpp, 'shps',fn,'a.name',driver='ESRI Shapefile',
overwrite_layer=TRUE)
  foo - readOGR('shps',fn)
}

print(sessionInfo())
 end script


And here is the output when I source the script in a fresh R session

-
source('try-problem.r')
Loading required package: sp
Loading required package: rgdal
rgdal: (SVN revision 360)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.9.0, released 2011/12/29
Path to GDAL shared files:
/Library/Frameworks/R.framework/Versions/2.15/Resources/library/rgdal/gdal
Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
Path to PROJ.4 shared files:
/Library/Frameworks/R.framework/Versions/2.15/Resources/library/rgdal/proj

R version 2.15.2 (2012-10-26)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] rgdal_0.7-24 sp_1.0-2

loaded via a namespace (and not attached):
[1] grid_2.15.2 lattice_0.20-10

OGR data source with driver: ESRI Shapefile
Source: shps, layer: tmpp1
with 10 features and 5 fields
Feature type: wkbPoint with 2 dimensions
OGR data source with driver: ESRI Shapefile
Source: shps, layer: tmpp2
with 10 features and 5 fields
Feature type: wkbPoint with 2 dimensions

--- output from iterations 2 through 12 omitted ---

OGR data source with driver: ESRI Shapefile
Source: shps, layer: tmpp13
with 10 features and 5 fields
Feature type: wkbPoint with 2 dimensions
OGR data source with driver: ESRI Shapefile
Source: shps, layer: tmpp14
with 10 features and 5 fields
Feature type: wkbPoint with 2 dimensions
Error in writeOGR(tmpp, shps, fn, a.name, driver = ESRI Shapefile,  :
  Layer creation failed


## then manually:

print(sessionInfo())
Error in gzfile(file, rb) : cannot open the connection
In addition: Warning message:
In gzfile(file, rb) :
  cannot open compressed file
'/Library/Frameworks/R.framework/Versions/2.15/Resources/library/rgdal/Meta
/package.rds', probable reason 'Too many open files'



I'm beginning to suspect this has to do with gzcon or gzfile
connections (not sure of correct terminology).
?open also says,

-- quote --
 'close' closes and destroys a connection.  This will happen
 automatically in due course (with a warning) if there is no longer
 an R object referring to the connection.
-- end quote --

But ?showConnections says

-- quote --
... However, if there is no R level object referring to
 the connection it will be closed automatically at the next garbage
 collection (except for 'gzcon' connections).
-- end quote --


So if read/write OGR involve zipping and unzipping, I guess the
connections don't get closed. (?)

Anyway, it would be very nice if I didn't have to quit and restart my R
sessions every so often when I'm working with my shapefiles.


Thanks
-Don

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

Re: [R-sig-Geo] R and GIS

[MacQueen, Don]

Here's a book:

The Geospatial Desktop
Open Source GIS  Mapping

by Gary Sherman

isbn 978-0-9868052-1-9

http://www.locatepress.com/

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 11/7/12 11:17 AM, marco milella vru...@gmail.com wrote:

Dear all,
I'm interested in the analysis of spatial data with R and I'm completely
new to the world of GIS. However, I would like to try to use R and GIS
data
for analyzing archaeological data. On the web I found several links, and
it
seems that a good combination could be  the use of R with GRASS. However,
I'm a little bit confused due to the big amount of possible links, guides,
manuals. Do you have any suggestion (online material or books) for a not
experienced R user and completely ignorant in GIS?
Thanks for any feedback
marco

   [[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] as() and conversion from character to factor

[MacQueen, Don]

Thank you.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 10/7/12 8:26 AM, Edzer Pebesma edzer.pebe...@uni-muenster.de wrote:

For now, you should be able to prevent this by starting your script with
an

options(stringsAsFactors = FALSE)

but in future releases, sp will not do this conversion.

Thanks,

On 10/04/2012 09:31 PM, MacQueen, Don wrote:
 I just noticed that in the following sequence the final step to
 SpatialGridDataFrame converts the character column to factor, whereas
the
 previous ones do not. Source the following little script, and it will
 hopefully show what I'm seeing. Shouldn't this be consistent?
 
 Thanks
 -Don
 
 
 require(sp)
 require(maptools)
 
 tstdf - expand.grid(x=1:3, y=1:3)
 tstdf$ch - rep('a',nrow(tstdf))
 
 tstsp - tstdf
 coordinates(tstsp) - c('x','y')
 proj4string(tstsp) - CRS('+init=epsg:32610')
 
 tstpx - as(tstsp, SpatialPixelsDataFrame)
 
 tstgrd - as(tstpx, SpatialGridDataFrame)
 
 print(class(tstdf$ch))
 print(class(tstsp$ch))
 print(class(tstpx$ch))
 print(class(tstgrd$ch))
 
 
 #
 I get:
 
 print(class(tstdf$ch))
 [1] character
 print(class(tstsp$ch))
 [1] character
 print(class(tstpx$ch))
 [1] character
 print(class(tstgrd$ch))
 [1] factor
 #
 
 
 
 
 
 
 
 sessionInfo()
 R version 2.15.1 (2012-06-22)
 Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
 
 locale:
 [1] C
 
 attached base packages:
 [1] stats utils datasets  grDevices graphics  methods   base
 
 other attached packages:
 [1] maptools_0.8-16 lattice_0.20-6  foreign_0.8-50  sp_1.0-0
 
 loaded via a namespace (and not attached):
 [1] grid_2.15.1 rmacq_1.1-8
 
 
 
 
 
 
 
 

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics  e.pebe...@wwu.de

___
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] as() and conversion from character to factor

[MacQueen, Don]

I just noticed that in the following sequence the final step to
SpatialGridDataFrame converts the character column to factor, whereas the
previous ones do not. Source the following little script, and it will
hopefully show what I'm seeing. Shouldn't this be consistent?

Thanks
-Don


require(sp)
require(maptools)

tstdf - expand.grid(x=1:3, y=1:3)
tstdf$ch - rep('a',nrow(tstdf))

tstsp - tstdf
coordinates(tstsp) - c('x','y')
proj4string(tstsp) - CRS('+init=epsg:32610')

tstpx - as(tstsp, SpatialPixelsDataFrame)

tstgrd - as(tstpx, SpatialGridDataFrame)

print(class(tstdf$ch))
print(class(tstsp$ch))
print(class(tstpx$ch))
print(class(tstgrd$ch))


#
I get:

 print(class(tstdf$ch))
[1] character
 print(class(tstsp$ch))
[1] character
 print(class(tstpx$ch))
[1] character
 print(class(tstgrd$ch))
[1] factor
#







 sessionInfo()
R version 2.15.1 (2012-06-22)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] C

attached base packages:
[1] stats utils datasets  grDevices graphics  methods   base

other attached packages:
[1] maptools_0.8-16 lattice_0.20-6  foreign_0.8-50  sp_1.0-0

loaded via a namespace (and not attached):
[1] grid_2.15.1 rmacq_1.1-8








-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062

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

Re: [R-sig-Geo] Apparent rotation when plotting a SpatialGridDataFrame?

[MacQueen, Don]

Thank you, Roger.

I see that in essence I separated my coords from my data,
in effect re-ordered the coords by making a grid topology,
and then rejoined the coords with the data in the original
order.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 10/1/12 11:45 AM, Roger Bivand roger.biv...@nhh.no wrote:

On Mon, 1 Oct 2012, MacQueen, Don wrote:


 I have simple data coming to me from an outside source.
 x,y are coordinates in a UTM projection, and z is a
 measured value at each location.

 I make a SpatialPointsDataFrame from it, plot it with spplot(),
 and it appears as expected.

 I convert to a SpatialGridDataFrame, plot using spplot(),
 it appears to have been rotated 90 degrees.

 See example code to illustrate below.

 Am I doing something wrong?


No, this is a consequence of code cleanup, revision 1140 in R-forge SVN:

https://r-forge.r-project.org/scm/viewvc.php/pkg/sp/R/Class-SpatialGrid.R?
root=rspatialr1=1140r2=1204

removing coords and grid.index slots from the definition of
the SpatialGrid class. This split SpatialPixels and SpatialGrid cleanly,
but meant that you can't smuggle the data slot across unpunished. This
works:

slot(tstsp, data)$z
tstpxdf - as(tstsp, SpatialPixelsDataFrame)
slot(tstpxdf, data)$z
tstgrdf - as(tstpxdf, SpatialGridDataFrame)
slot(tstgrdf, data)$z # showing the re-ordering
tmppl4 - spplot(tstgrdf)
print(tmppl4)


Hope this helps,

Roger


 I used the same sequence of conversions about 7 months ago on
 a different set of data and did not see this problem. Inspecting
 the two versions of the data does not suggest a problem with
 the conversion process. The bounding box is consistent before
 and after the conversion.

 Here is a simple example to illustrate, followed by sessionInfo()
 It should be possible to just source the example code into an
 R session.

 ##
 ## simple example to illustrate the question
 ##
 require(sp)

 ## example data
 tstdf - expand.grid(x=1:3, y=1:3)
 tstdf$z - 5
 tstdf$z[tstdf$x==tstdf$y] - 10

 ## move points to somewhere actually in the zone
 tstdf$x - tstdf$x +  608000
 tstdf$y - tstdf$y + 4167900

 ## convert to spatial points data frame; UTM projection
 tstsp - tstdf
 coordinates(tstsp) - c('x','y')
 proj4string(tstsp) - CRS('+init=epsg:32610')
 print(bbox(tstsp))

 ## plot
 tmppl1 - spplot(tstsp)
 print(tmppl1)


 ## convert to spatial grid
 tstgr - points2grid(tstsp)
 tstsgr - SpatialGrid(tstgr,
 proj4string=CRS('+init=epsg:32610'))
 print(bbox(tstsgr))
 tstsgrd - SpatialGridDataFrame(tstgr,
data=tstsp@data,
proj4string=CRS('+init=epsg:32610'))
 print(bbox(tstsgrd))

 ## plot -- compare with first plot
 tmppl2 - spplot(tstsgrd)
 print(tmppl2)

 ##
 ## sessionInfo()
 ##
 R version 2.15.1 (2012-06-22)
 Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

 locale:
 [1] C

 attached base packages:
 [1] stats graphics  grDevices utils datasets  methods   base

 other attached packages:
 [1] sp_1.0-0

 loaded via a namespace (and not attached):
 [1] grid_2.15.1lattice_0.20-6


 ---
 By the way, the reason for converting to a spatial grid is
 so that I can then convert it to a SpatialPolygonsDataFrame,
 spTransform() to long/lat and export to KML,
 then display in GoogleEarth with a raster-like appearance.


 There are no doubt other and likely better ways to accomplish
 the same end goal, a raster-like display in GoogleEarth,
 but this is a way I found.
 



-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: roger.biv...@nhh.no


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

Re: [R-sig-Geo] New (?) use of polypath() in package sp causing a problem?

[MacQueen, Don]

(see below)


On 8/10/12 1:04 AM, Roger Bivand roger.biv...@nhh.no wrote:

On Fri, 10 Aug 2012, MacQueen, Don wrote:

 I have encountered a problem I've never seen before.

 The problem is:

 x11(type='Xlib')
 plot(sfbound)
 Warning message:
 In polypath(x = mcrds[, 1], y = mcrds[, 2], border = border, col = col,
:
 Path drawing not available for this device

 And nothing is plotted.

 The object sfbound is a SpatialPolygonsDataFrame object with one simple
 polygon in it.
 Previously I have had no problem plotting such objects with
 x11(type='Xlib').x11(type='cairo') succeeds.


 It looks very likely that this relates to this entry in NEWS for R
2.12.0:

 -- quote --
 Added support for polygons with holes to the graphics
 engine. This is implemented for the pdf(),
 postscript(), x11(type=cairo), windows(),
 and quartz() devices (and associated raster formats),
 but not for x11(type=Xlib) or xfig() or
 pictex(). The user-level interface is the
 polypath() function in graphics and
 grid.path() in grid.
 -- end quote --

 That is, presuming a relatively recent change in sp in which polyplot()
is
 now used and previously was not, the change described in NEWS would
 explain why I'm seeing this.



 The problem is, I have found that for plots in which a very large number
 of objects are being plotted, the Xlib version of X11() is much faster
 than the cairo version. And I frequently plot spatial objects with
 thousands of line segments in them, where the performance difference is
 substantial, to the point of making use of the cairo version
intolerable.

 I guess the question is -- am I correct in this analysis?
 If so, do I have any alternative when I need to plot spatial objects
with
 huge numbers of points/lines/polygons in them?

Don,

The analysis is correct (r1248). If the hole status of Polygon objects is
known to be correct, and the par(bg=) and pbg= values are set to match,
the plots will be the same. polypath analyses the polygon rings itself
and 
provides transparency in polygon painting - it doesn't paint in holes.
The 
legacy polygon function first painted the whole exterior ring, then
overpainted holes with the pbg= colour, which is transparent by
default. 
For many purposes, polypath is preferable, but as you know, cairo is much
slower than X11. So using X11 ought still to be OK, in the plot method
set 
usePolypath=FALSE. We can add back an option and set the argument by
default to that option if that would help.

The details for the plot method in the code:

https://r-forge.r-project.org/scm/viewvc.php/pkg/sp/R/SpatialPolygons-disp
layMethods.R?root=rspatialview=log

and for spplot methods look in the same repository.

Hope this helps,

Roger

It does help, especially the pointer to usePolypath=FALSE, which lets me
plot with X11/Xlib. Thank you.

In terms of adding back an option, would that be a global option of some
sort, so that users like me don't have to supply usePolypath=FALSE every
time?  That would indeed be helpful, if it's not too much trouble.

-Don



 Thanks
 -Don


 str(sfbound)
 Formal class 'SpatialPolygonsDataFrame' [package sp] with 5 slots
  ..@ data   :'data.frame': 1 obs. of  1 variable:
  .. ..$ name: Factor w/ 1 level Boundary.kml: 1
  ..@ polygons   :List of 1
  .. ..$ :Formal class 'Polygons' [package sp] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package sp] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] -121.7 37.7
  .. .. .. .. .. .. ..@ area   : num 3.71e-06
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:8, 1:2] -122 -122 -122 -122 -122
 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt: num [1:2] -121.7 37.7
  .. .. .. ..@ ID   : chr 1
  .. .. .. ..@ area : num 3.71e-06
  ..@ plotOrder  : int 1
  ..@ bbox   : num [1:2, 1:2] -121.7 37.7 -121.7 37.7
  .. ..- attr(*, dimnames)=List of 2
  .. .. ..$ : chr [1:2] x y
  .. .. ..$ : chr [1:2] min max
  ..@ proj4string:Formal class 'CRS' [package sp] with 1 slots
  .. .. ..@ projargs: chr +proj=longlat +ellps=GRS80 +datum=NAD83
 +no_defs +towgs84=0,0,0


 sessionInfo()
 R version 2.15.1 (2012-06-22)
 Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

 locale:
 [1] C

 attached base packages:
 [1] stats graphics  grDevices utils datasets  methods   base

 other attached packages:
 [1] sp_0.9-99   rmacq_1.1-8

 loaded via a namespace (and not attached):
 [1] grid_2.15.1lattice_0.20-6 tools_2.15.1





-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: roger.biv...@nhh.no


-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062

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

Re: [R-sig-Geo] How to identify unknown Coordinate systems?

[MacQueen, Don]

I had data in an unknown projection -- but I knew exactly where it was
supposed to be located.

What I did was make a guess at its projection, transform it to lat/long,
export to KML, then open in Google Earth as an easy way to check whether
it was where it was supposed to be. Repeat until satisfied.

However, from other information I could limit the list of possible
projections, and that obviously helped.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 7/24/12 9:36 AM, Barry Rowlingson b.rowling...@lancaster.ac.uk
wrote:

On Tue, Jul 24, 2012 at 5:17 PM, Juan Tomas Sayago
juantomas.say...@gmail.com wrote:
 Hello,
 I have data in an unknown projection  but I know where more or less
where
 it is located, but I don't know from what to transform to what I want to
 use.

First of all it could be one of the global systems like UTM zones.
Find out which UTM zone you might expect it to be and see what the
coordinates are of the region you think it might be.

Also, make sure you haven't got a .prj file with your shapefile, or
that the CRS is hidden in some metadata file.

Or try looking up the geographical area here:

http://www.spatialreference.org/

 but you might have to be creative - searching for 'Great Britain'
doesnt find the British national grid system, but searching for
'British' does. That site is also good for checking UTM zones. Note
there's a lot of unofficial CRS in there.

 Another useful thing is to figure out if the coordinates are metres
or km or miles, feet and inches - how far apart are your points in the
coordinate system?

 If all that fails, if you have the location of some of your points in
a known coordinate system, you can probably work out the transform
yourself, but it's likely to be imprecise.

 Finally, try beating the person who gave you the data until they tell
you! Or being nice to them.

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

Re: [R-sig-Geo] shapefile on the wrongly placed

[MacQueen, Don]

May we assume you are using R to do this?

 - Read it into R using an appropriate function from an R package. (I use
readOGR from rgdal)

 - Make sure that it is assigned the correct projection.

 - Transform it to the desired projection using spTransform()

This has always worked for me.
The help page for spTransform has examples.
  

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 5/23/12 11:02 AM, Juan Tomas Sayago juantomas.say...@gmail.com
wrote:

Dear list
I have this shape file created on one location using a specific
projection
and when I try to put it into WGS1984 it appears in a different place
where
it should be. How do I fix it? What can I do?
Juan

-- 
Juan Tomás Sayago Gómez
Graduate Research Assistant
West Virginia University - RRI
886 Chestnut Ridge Road, Room 520
P.O. Box 6825
Morgantown, WV 26506-6825

   [[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] rgdal package for mac

[MacQueen, Don]

I have rgdal installed on two 10.6.8 machines -- one regrettably still
using R 2.14.1 and the other R 2.15.0.

Over time I have used both the CRAN Extras rgdal and the kyngchaos rgdal.
Currently I have the CRAN Extras version (0.7-8, 2012-01-18) installed,
but if there is now a version mismatch as Roger suggested might be the
case, I would try the kyngchaos version of rgdal.

My basic process for GIS support on my Macs is to
  1) download and install the necessary frameworks from kyngchaos
  2) install rgdal

 GIS software for Mac (Snow Leopard and Lion)
http://hub.qgis.org/projects/quantum-gis/wiki/Download#3-MacOS-X
http://www.kyngchaos.com/software/qgis
http://www.kyngchaos.com/software/frameworks

The following describes the 2.15.0 machine.


 listing of my GIS support directory, from kyngchaos qgis and
frameworks pages
GIS-support-2012-04[8]% ls -1 *.dmg *.zip
 FreeType_Framework-2.4.9-1.dmg
 GDAL_Complete-1.9.dmg
 GRASS-6.4.2-3-Snow.dmg
 GSL_Framework-1.15-2.dmg
 Qgis-1.7.4-4.dmg
 cairo_Framework-1.10.2-3a-snow.dmg
 rgdal-0.7.8-1.zip   (this is the kyngchaos version of rgdal)

 Then in R:

## loading rgdal
 require(rgdal)
Loading required package: rgdal
Loading required package: sp
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.9.0, released 2011/12/29
Path to GDAL shared files:
/Library/Frameworks/GDAL.framework/Versions/1.9/Resources/gdal
Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 470]
Path to PROJ.4 shared files: (autodetected)

 sessionInfo()
R version 2.15.0 (2012-03-30)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] rgdal_0.7-8 sp_0.9-98   rmacq_1.1-6

loaded via a namespace (and not attached):
[1] grid_2.15.0lattice_0.20-6 tools_2.15.0


 
Note that my sp is not (quite) up to date, as Roger mentioned sp 0.9.99
somewhere in this long email chain.


-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 5/18/12 9:17 AM, de Larminat, Pierre p...@mpifg.de wrote:

Le 18 mai 2012 à 18:11, Roger Bivand a écrit :

 On Fri, 18 May 2012, de Larminat, Pierre wrote:
 
 Le 18 mai 2012 à 17:52, Roger Bivand a écrit :
 
 On Fri, 18 May 2012, de Larminat, Pierre wrote:
 
 
 
 Le 18 mai 2012 à 17:14, Roger Bivand a écrit :
 
 On Fri, 18 May 2012, de Larminat, Pierre wrote:
 
 Hello,
 
 I have the same problem than:
 [R-sig-Geo] rgdal package for mac
 Gabriele Cozzi gab.cozzi at gmail.comhttp://gmail.com
 Thu May 17 10:33:21 CEST 2012
 
 I have already done:
 setRepositories(ind=1:2)
 install.packages(rgdal)
 as said on the CRAN page for rgdal.
 
 I have restarted the computer and R.
 When I entered the lines again, I received the same message that I
had received before:
 
 setRepositories(ind=1:2)
 install.packages(rgdal)
 --- SVP s?lectionner un miroir CRAN pour cette session ---
 essai de l'URL
'http://www.stats.ox.ac.uk/pub/RWin/bin/macosx/leopard/contrib/2.15/
rgdal_0.7-8.tgz'
 Content type 'application/x-gzip' length 13810017 bytes (13.2 Mb)
 URL ouverte
 ==
 downloaded 13.2 Mb
 
 
 Les packages binaires t?l?charg?s sont dans
 
/var/folders/zT/zTnmGA9q2RaOuE+F76E3ETI/-Tmp-//RtmpKnHXly/downlo
aded_packages
 library(rgdal)
 Le chargement a n?cessit? le package : sp
 
 This does seem to imply that sp is necessary, doesn't it? Have you
 installed sp? Messages are given to be read, rather than ignored.
 
 Roger
 
 Yes I had installed sp. The message merely states that I had not
charged the sp library and that R had to do it by itself. When I do
run *library(sp)* in time, I still get:
 
 OK. The next error message refers to your system having
libcurl.4.dylib
 provides version 6.0.0, but that the rgdal OSX binary requires =
7.0.0:
 
 
 Error in dyn.load(file, DLLpath = DLLpath, ...) :
 impossible de charger l'objet partag?
'/Library/Frameworks/R.framework/Versions/2.15/Resources/library/rgd
al/libs/i386/rgdal.so':
 
dlopen(/Library/Frameworks/R.framework/Versions/2.15/Resources/libra
ry/rgdal/libs/i386/rgdal.so, 6): Library not loaded:
/usr/lib/libcurl.4.dylib
 Referenced from:
/Library/Frameworks/R.framework/Versions/2.15/Resources/library/rgda
l/libs/i386/rgdal.so
 Reason: Incompatible library version: rgdal.so requires version
7.0.0 or later, but libcurl.4.dylib provides version 6.0.0
 
  here. Are you running an older OSX version?
 
 Roger
 
 I am running R 2.15 on MacOS 10.6.8 and I do not know what
libcurl.4.dylib is.
 
 Please ask on R-sig-mac, or get local help; users may not know what
system 
 resources are, but do suffer the consequences. The cause is probably
that 
 the rgdal CRAN Extras OSX binary was built on 10.7, and you are running
 10.6 - the background for building on a newer system was given in an
 earlier thread.
 
 Roger
 

OK, thank 

Re: [R-sig-Geo] Storing Multiple spplots in List

[MacQueen, Don]

Please make a small reproducible example.
There is no way anyone can figure out the problem with what you have
supplied.

What is telling you they are blank?
Do your spplot() commands produce a plot when executed outside the
function?

This works for me:

tmp - list()
length(tmp) - 2

## I used the first example in ?spplot to do this:
tmp[[1]] -  spplot(  ## args )
tmp[[2]] -  spplot(  ## args )

## then I get plots from any of these:
lapply(tmp, print)

## or
print( tmp[[1]])

## or even,
print(tmp)


-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 4/13/12 9:01 PM, Mike Sukmanowsky mike.sukmanow...@gmail.com wrote:

Hi all,

I have the function below where I am attempting to create a list of
spplots
for later use.  The problem is that the list come back blank for all
elements.  Why doesn't spplot return proper plots that could be used via
print()?
Please see the assignment to the plots list below - plots[[letter]] -
spplot()


plotFSAs - function(spatialData, field) {
plots - list()
colours - brewer.pal(8, Greens)
 for(letter in LETTERS) {
inLetter - substr(spatialData@data$CENSUS_FSA, 1, 1) %in% letter
fsa.spatial - subset(spatialData, inLetter)
 classes - classIntervals(fsa.spatial@data[,c(field)], length(colours),
style=equal)
 plots[[letter]] - spplot(fsa.spatial,
field,
col.regions=colours,
at=classes$brks,
col=black,
lwd=0.8,
main=paste(FSAs Beginning with , letter),
cex.main=0.3,
colorkey=list(
labels = list(
at=classes$brks,
labels=round(classes$brks, 2)
)
)
)
}
return(plots)
}


-- 
Mike Sukmanowsky

   [[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] subset() for spatialPointsDataFrame

[MacQueen, Don]

Try

  points.df[points.df$value  0, c('column1','column')]

(or possibly using   points.df@data$value instead of points.df$value

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 3/6/12 7:09 AM, Johannes Radinger jradin...@gmx.at wrote:

Hi,

I just wanted to know how I have to use the
subset() command for a SpatialPointsDataFrame correctly.

I'd like to get all records of the SpatialPointsDataFrame where the
column value is  0. Furthermore I want to select only two columns of
that dataframe with these records. The subset command would be:

subset(points.df, value  0 , select = c(column1, column2))

When I do that I get following error:

Error in subset.SpatialPointsDataFrame(points.df, value  0, select =
c(column1,  : 
  Object 'value' not found

I also tried:
subset(points.df, points.df[value]  0 , select = c(column1, column2))

Error in points.df[value]  0 :
  Vergleich (6) ist nur fur atomare und Listentypen moglich (sorry is
in German)...

Best regards,
Johannes

-- 

Jetzt informieren: http://mobile.1und1.de/?ac=OM.PW.PW003K20328T7073a

___
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] spplot handling of overlapping points?

[MacQueen, Don]

I have a SpatialPointsDataFrame object in which many points are very close
together, such that the markers (plotting characters) tend to overlap. (Of
course, this depends on marker size and the scale at which I plot; if I
zoom in there is less overlap.)

It appears that when markers overlap, spplot() places a marker associated
with the larger value after, and therefore on top of, a marker associated
with a smaller value.

Am I correct? Or more generally, what is the algorithm that determines the
order in which markers are added?


I have searched ?spplot and related help pages and haven't found an
explanation (at least, not yet).

I can add that I don't believe markers are placed in the order in which
they appear in the SpatialPointsDataFrame. I say this because when I
attempt to reproduce an spplot using base graphics plot() I have to sort
from smallest to largest value to succeed.

I can probably provide a small reproducible example if necessary, but I'm
hoping it's not necessary.

Here are my actual commands.

  tmps is the SpatialPointsDataFrame.
  tmp is coordinates(tmps)
(they have the same number of rows in the same order)

Note that I'm using the cuts argument to break a continuous variable
('cpm2') into bins. I have carefully matched the colors in tmps$col2 with
the cbin.cols object passed to spplot(), and the break points for the bin
boundaries, so I believe that everything else that could affect the final
appearance, other than the order in which the markers are placed, is
controlled.

#1 using spplot
spplot(tmps,c('cpm2'),
   key.space='right',
   legendEntries=cbin.lbls,
   cuts=cbin.brks,
   col.regions=cbin.cols,
   cex=0.4)

#2 using base graphics
plot(tmp[,1],tmp[,2], asp=1, cex=0.6, pch=16, col=tmps$col2)



Thanks
-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062

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

Re: [R-sig-Geo] write/readOGR issue with date / time strings

[MacQueen, Don]

You don't truly have a date/time field, even though it looks like it. You
have a factor. 

So the first thing I would try is to convert the date/time field from
factor to character using format() or as.character() before writing the
shapefile. (That is, do the conversion myself rather than relying on
writeOGR's internal conversion; see the help page for writeOGR).

I'm a little bothered by the fact that Lat and Lon are factors; they
shouldn't be. This suggests that somewhere in your data is at least one
row where lat or long is not a valid number. But writeOGR doesn't use
these particular Lat and Lon columns to write coordinates to the
shapefile, so maybe it doesn't matter.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 11/11/11 12:54 PM, Corrie Curtice corrie.curt...@duke.edu wrote:

Hello, 

Apologies if this has already been addressed somewhere, I did a brief
search
of archives but didn't find quite this issue.

I'm writing out an ESRI shapefile.  My spdf has a date/time field.  Looks
like this:

 head(spdfUTM@data)
   krillGMTtime  Lat  Lon
1   2010-05-12 12:34:21 -64.67655969 -62.15040195
21  2010-05-12 12:35:12  -64.6771229  -62.1519511
38  2010-05-12 12:36:02 -64.67775863 -62.15340614
57  2010-05-12 12:36:53 -64.67838269 -62.15494829
78  2010-05-12 12:37:43 -64.67901497 -62.15647203
100 2010-05-12 12:38:37 -64.67973667 -62.15773224
 
All three fields are reported to be factors by str.  writeOGR appears
happy:

 writeOGR(spdfUTM,dd,layer=krillPoints-UTM,driver=ESRI
Shapefile,verbose=TRUE,overwrite=TRUE)
$object_type
[1] SpatialPointsDataFrame
$output_dsn
[1] /users/corriecurtice/documents/Data_2010/Shapefiles/
$output_layer
[1] krillPoints-UTM
$output_diver
[1] ESRI Shapefile
$output_n
[1] 11179
$output_nfields
[1] 3
$output_fields
[1] krillGMTtime Lat  Lon
$output_fclasses
[1] 4 4 4
$dataset_options
NULL
$layer_options
NULL
Warning message:
In writeOGR(spdfUTM, dd, layer = krillPoints-UTM, driver = ESRI
Shapefile,  :
  existing layer removed

When I read it right back in again, the date/time field is NAs.  This is
also true if I load the shapefile into ArcMap.

 foo - readOGR(dd,layer=krillPoints-UTM,verbose=TRUE)
OGR data source with driver: ESRI Shapefile
Source: /users/corriecurtice/documents/Data_2010/Shapefiles/, layer:
krillPoints-UTM
with 11179 features and 3 fields
Feature type: wkbPoint with 2 dimensions

 head(foo@data)
  krillGMTti  Lat  Lon
1   NA -64.67655969 -62.15040195
2   NA  -64.6771229  -62.1519511
3   NA -64.67775863 -62.15340614
4   NA -64.67838269 -62.15494829
5   NA -64.67901497 -62.15647203
6   NA -64.67973667 -62.15773224

Thoughts? Am I missing something obvious?

I can break it up into separate date and time fields, but since the field
is
simply a factor anyways I'm not sure how that would help.

Cheers,

Corrie

---
Corrie Curtice
Research Analyst
Marine Geospatial Ecology Lab
Nicholas School of the Environment, Duke University
http://mgel.env.duke.edu
em: corrie.curt...@duke.edu
ph: 252-504-7538




   [[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

[R-sig-Geo] rbind of SpatialPolygonDataFrame objects is different with or without rgdal package

[MacQueen, Don]

It appears that when I rbind() a pair of SpatialPolygonDataFrame objects
when the rgdal package has been loaded the resulting projection string
gets an embedded space character, whereas without rgdal it does not. The
space character breaks further use of rbind(). I don't think it matters
which one is first in the search path (I've tried it both ways).



Here's an example.

I have three simple spatial polygon data frame objects; each consists of a
single polygon. They have the same projection. I can send them as an
attachment if requested.

## in a new R session (objects created in an earlier session):

 require(sp)
Loading required package: sp

 proj4string(ex1.sp) == proj4string(ex2.sp)
[1] TRUE

 proj4string(ex1.sp) == proj4string(ex3.sp)
[1] TRUE

 tmp1 - rbind(ex1.sp, spChFIDs(ex2.sp,'ex2'))

 proj4string(ex1.sp) == proj4string(tmp1)
[1] TRUE

 sessionInfo()
R version 2.13.1 (2011-07-08)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] sp_0.9-88

loaded via a namespace (and not attached):
[1] grid_2.13.1 lattice_0.19-30


 require(rgdal)
Loading required package: rgdal
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.8.0, released 2011/01/12
Path to GDAL shared files:
/Library/Frameworks/R.framework/Versions/2.13/Resources/library/rgdal/gdal
Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009, [PJ_VERSION: 470]
Path to PROJ.4 shared files:
/Library/Frameworks/R.framework/Versions/2.13/Resources/library/rgdal/proj



 tmp2 - rbind(ex1.sp, spChFIDs(ex2.sp,'ex2'))

 proj4string(ex1.sp) == proj4string(tmp2)
[1] FALSE

 proj4string(ex1.sp)
[1] +proj=lcc +lat_1=38.43 +lat_2=37.07
+lat_0=36.5 +lon_0=-120.5 +x_0=200.0001016 +y_0=50.0001016001
+ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs

## note the leading space character
 proj4string(tmp2)
[1]  +proj=lcc +lat_1=38.43 +lat_2=37.07
+lat_0=36.5 +lon_0=-120.5 +x_0=200.0001016 +y_0=50.0001016001
+ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs
+towgs84=0,0,0

 tmp3 - rbind(tmp2, spChFIDs(ex3.sp,'ex3'))
Error in checkCRSequal(dots) : coordinate reference systems differ


 sessionInfo()
R version 2.13.1 (2011-07-08)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] rgdal_0.7-1 sp_0.9-88

loaded via a namespace (and not attached):
[1] grid_2.13.1 lattice_0.19-30
 



 str(ex1.sp)
Formal class 'SpatialPolygonsDataFrame' [package sp] with 5 slots
  ..@ data   :'data.frame': 1 obs. of  1 variable:
  .. ..$ name: Factor w/ 1 level example1: 1
  ..@ polygons   :List of 1
  .. ..$ :Formal class 'Polygons' [package sp] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package sp] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 6223669 2081123
  .. .. .. .. .. .. ..@ area   : num 20287432
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:8, 1:2] 6221449 6226373 6225988
6223664 6223594 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt: num [1:2] 6223669 2081123
  .. .. .. ..@ ID   : chr 1
  .. .. .. ..@ area : num 20287432
  ..@ plotOrder  : int 1
  ..@ bbox   : num [1:2, 1:2] 6220680 2078695 6226373 2083511
  .. ..- attr(*, dimnames)=List of 2
  .. .. ..$ : chr [1:2] x y
  .. .. ..$ : chr [1:2] min max
  ..@ proj4string:Formal class 'CRS' [package sp] with 1 slots
  .. .. ..@ projargs: chr +proj=lcc +lat_1=38.43
+lat_2=37.07 +lat_0=36.5 +lon_0=-120.5 +x_0=200.0001016
+y_0=50.00010160| __truncated__





-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062

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