Domenico Giusti | 30 Mar 15:39 2015
Picon

R 3.1.3 - GNU/Linux - readOGR -- Layer not found

Dear list,

I cannot read some of the layers (table and view) from a
Postgres/PostGIS DB.

Error in ogrInfo.

> ogrListLayers("PG:dbname=dbname")

doesn't list all the tables and views. e.g. I have two views with both a
3D point geom, but only one of them is in the list.

What could be the problem?

Thanks,

> sessionInfo()
R version 3.1.3 (2015-03-09)
Platform: i486-pc-linux-gnu (32-bit)
Running under: Debian GNU/Linux 7 (wheezy)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
LC_TIME=en_US.UTF-8
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8
LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C

[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8
LC_IDENTIFICATION=C
(Continue reading)

Doan Thanh Thuy | 30 Mar 16:14 2015
Picon

Problem giving points back to original location in raster in R

Dear all,
I am doing Latin Hypercube Sampling (LHS) in R to find the best location to
sample new soil data points. Firstly, I have 5 raster maps of NDVI, Slope,
TWI, Altitude, Geology which were then imported in R. All those data were
collected into a data frame in R to be the input of LHS. The results of LHS
will be the list of points at which I should take soil sample.
However, the problem is that I don't know how to put those resultant points
back to the raster map. Can anyone help me with this kind of stupid
question please? I am at the very early stage in R :)
Thanks in advance!
Best regards,
Thuy

Doan Thanh Thuy
Department of Land Information System
Faculty of Land Management
Vietnam national University of Agriculture
Mobile: +841689686205
Email: doanthanhthuy209 <at> gmail.com, ThanhThuy.Doan <at> UGent.be

	[[alternative HTML version deleted]]
Leask, Graham | 29 Mar 12:50 2015
Picon
Picon

Heat mapping within shape file

What is a straightforward way to produce a heat map with the heat map strictly constrained within but
filling shape file borders?

If someone can illustrate an example in R I would be most grateful.

Best wishes

Graham
Mike | 28 Mar 19:01 2015
Picon

Complete traversals: Postman problem with shape files

Hi all, 

I'm relatively new to the listserv, but learning a lot.  Please excuse me if this is a semi-related question.

As part of a public health project to survey an area for tobacco retailers, we're attempting to investigate a systematic, complete, and reasonably efficient (not necessarily perfect) traversal of all the roads in a certain catchment area.  This is an application of the classic "postman" problem - traversing all the edges of a graph - related to the canonical Traveling Postman Problems / TSPs.   Specifically, we'd be taking regions of VA and breaking that up into volunteer catchment areas - a county or two - to drive the area.

I've a proof of concept plot and code below that takes a shape file of regions (in this case, of VA counties) and traverses them in a reasonably efficient way... but this uses centroids of shapes and creates a distance matrix off of that, instead of working with the road shape file itself, which is currently stumping me.

(Below is me muddling through a test case using the counties as a whole, code and picture.)

I know this is not at all a trivial problem.  Any guidance or expertise here?  I could make accessible more useful test cases/datasets (like a single county's roads) if that were helpful.  Or perhaps there are packages that do this out of the box besides what I'm using?

Best wishes and much respect,
Mike Fliss
UNC Epidemiology PhD student.




library("maptools")
library("spdep")
library("sp")
library("TSP")
#data("USCA312")

VAcounties = readShapeSpatial("E:/Mike/GIS/Base shape files/VA Counties/VA Counties.shp")
plot(VAcounties, border="grey")
neighbors <-poly2nb(VAcounties)
plot(neighbors, coordinates(VAcounties), line="grey", add=TRUE)

VApoints = read.csv("E:/Dropbox/Community/Driving VA/vapoints.csv")
VA.SPDF = SpatialPointsDataFrame(VApoints[,2:3], VApoints[1], proj4string = CRS("+proj=longlat"))

vamatrix <- read.csv("E:/Dropbox/Community/Driving VA/vacountymatrix.csv", stringsAsFactors=FALSE)
va2 = as.matrix(vamatrix)
dimnames(va2) = list(names(vamatrix), names(vamatrix))

vatsp = TSP(va2)
vatsp=insert_dummy(vatsp, label="cut")

#tour = solve_TSP(vatsp, method="repetitive_nn", control=list(start=1))
tour = solve_TSP(vatsp, method="2-opt", control=list(rep=100))
path = cut_tour(tour, "cut")

head(labels(path))

plot(VAcounties, border="grey")
plot(VA.SPDF, pch=16, cex=.4, col="red", add=T)
path_line = SpatialLines(list(Lines(list(Line(VA.SPDF[path,])), ID="1")))
plot(path_line, add=T, col="black")
points(VA.SPDF[c(head(path,1), tail(path,1)),], pch = 19, col = "black")

--
---
Mike Dolan Fliss, MSW
UNC-CH Epidemiology PhD student
NC public health advocate!
-------------------
"We work on ourselves in order to help others, 
but also we help others in order to work on ourselves."
  - Pema Chodron

“Upon this gifted age, in its dark hour
Rains from the sky a meteoric shower
Of facts…they lie, unquestioned, uncombined.
Wisdom enough to leach us of our ill
Is daily spun; but there exists no loom
To weave it into a fabric.”
  - Edna St. Vincent Millay, 1939

There are two kinds of people in the world:
Those who can extrapolate from incomplete data.
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo <at> r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Walter Anderson | 27 Mar 23:45 2015
Picon

Calculating area of polygons created from a spatial intersect

Hello all,

I am attempting to automate an analysis that I developed with ArcInfo
using R and the gdal and geos packages (or any other) if possible.

Here is the basic process

I have a shape file (lines) that defines the limits of all of the
projects with each project having a unique identifier.

I have another shape file (polys) that contains total population and
low income population and represent Census block groups.  This shape
file has an area field which has the acreage of the total block group.

Process

Step 1.
I then buffer these project lines to create a second shape file that
represents the 'footprint' of the project. (Creates polys).

Step 2.
In ArcInfo, I perform an intersection of the two shape files
(footprint and census blocks) and this creates a third shape file
which has a unique polygon for every project/census block intersection.

Step 3.
I then perform an area calculation (acres) on this new poly shape file
and use this calculated area divided by the original area of the
associated census block group to apportion the two population datum to
this new polygon.

Step 4.
Finally, I sum the two population datums for each of the projects from
the attribute table of this final shape file.

When I try to replicate the above procedure I run into a problem with
Step 2 when I use what I think is the appropriate command:

gIntersects(buffered_projects, census_blocks, byID=TRUE)

This command is producing a matrix of each project/census block
combination and only providing me a true/false indication.  Is there
any way to replicate the process from ArcInfo that I outlined above
within R?

Walter Anderson
Andrew Duff | 27 Mar 21:30 2015
Picon

Re: Help with latlong to UTM conversion when UTM zones are different

A number of field folks prefer UTM because 

-it matches legacy paper USGS quad map series traditionally used for field navigation
-units are in meters and can be used to gauge field distances from a coordinate readout 

> On Mar 27, 2015, at 10:48 AM, Michael Sumner <mdsumner <at> gmail.com> wrote:
> 
> There is no good natural reason to use UTM, it mistifies me why our
> community tolerates this bizarre default. I always use a local equal-area
> projection unless some other compromise dictates a different choice.
> Cheers, Mike
> 
> On Fri, 27 Mar 2015 21:28 Barry Rowlingson <b.rowlingson <at> lancaster.ac.uk>
> wrote:
> 
>> If you have lat-long data that crosses two UTM zones then its
>> generally okay to just pick *one* and transform all the points to
>> that. Use the one that has the most points in. Basically use the UTM
>> zones as guidelines to pick one UTM zone coordinate system. Unless
>> your data spans several zones and you want quite high accuracy of
>> distance measurements. Some points bleeding over into an adjacent zone
>> are no problem.
>> 
>> All projections are approximations to the earth's spheroid, so points
>> that are within a single UTM zone have some distortion in their
>> distance or angle relationships. Transforming points that are within
>> an adjacent UTM zone is just an extension of that distortion. You can
>> compute the precise distance error if you want for the furthest points
>> by comparing with the geodesic distance.
>> 
>> Alternatively you might find there is a coordinate system that spans
>> your dataset nicely - often when a country or an island or a region
>> crosses UTM zones there is an official coordinate system defined that
>> is used by the authorities there.
>> 
>> Also alternatively, there's nothing to stop you defining a transverse
>> mercator system based on the centre of your data.
>> 
>> Barry
>> 
>> 
>> 
>> On Fri, Mar 27, 2015 at 7:44 AM, moses selebatso <selebatsom <at> yahoo.co.uk>
>> wrote:
>>> Hello
>>> I have animal movement data that I have converted from Lat/Long to UTM,
>> unfortunately the data falls in two UTM zones (34S and 35S). For some
>> reason R cannot display both of them in the same window (the 35S data is
>> way off the expected location).
>>> The question is how do I convert the data such that R can correctly read
>> it?
>>> Moses SELEBATSO
>>> 
>>> (+267) 318 5219 (H) (+267) 716 393 70 (C)
>>>  (+267) 738 393 70 (C
>>>        [[alternative HTML version deleted]]
>>> 
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo <at> r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> 
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo <at> r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 
>    [[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo <at> r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Ben Tupper | 27 Mar 21:20 2015

Pole centric map with Raster?

Hello,

I'm learning how to use Raster* objects and ultimately hope to draw rasters from a (north) polar
perspective.  I can do something like this using tiling in ggplot2:

https://dl.dropboxusercontent.com/u/8433654/polar-map-with-matrix.png

I have hopes of doing similar using raster, sp and lattice.  But I'm drowning in information as neophytes are
prone to do.  Below is a self-contained example...  I'm not getting too far.  It creates two plots:

raster without using projectRaster: https://dl.dropboxusercontent.com/u/8433654/raster_no_project.png

raster using projectRaster: https://dl.dropboxusercontent.com/u/8433654/raster_with_project.png

Neither is what I would like.  How does one get the map centered on the pole and draw a (warped) raster using spplot?

Thanks!
Ben

##### START
library(raster)
library(maps)
library(maptools)
library(rgdal)

# create a Raster from 'volcano'
my_proj <- "+proj=stere +lat_0=90 +lon_0=-30 +x_0=0 +y_0=0" 
r <- raster(volcano,
   xmn = -90, xmx = 0, ymn = 50, ymx = 85, 
   crs = my_proj)
# class       : RasterLayer 
# dimensions  : 87, 61, 5307  (nrow, ncol, ncell)
# resolution  : 1.47541, 0.4022989  (x, y)
# extent      : -90, 0, 50, 85  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=stere +lat_0=90 +lon_0=-30 +x_0=0 +y_0=0 +ellps=WGS84 
# data source : in memory
# names       : layer 
# values      : 94, 195  (min, max)   

pr <- projectRaster(r, crs=my_proj)
# class       : RasterLayer 
# dimensions  : 123, 71, 8733  (nrow, ncol, ncell)
# resolution  : 1.48, 0.402  (x, y)
# extent      : -97.4, 7.68, 42.79, 92.236  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=stere +lat_0=90 +lon_0=-30 +x_0=0 +y_0=0 +ellps=WGS84 
# data source : in memory
# names       : layer 
# values      : 93.88372, 194.8761  (min, max)

# create a coastline and transform to the projection
mp <- map2SpatialLines(map(plot = FALSE, interior = FALSE), proj4string = CRS(my_proj))
mp <- spTransform(mp, CRSobj = CRS(my_proj))
# class       : SpatialLines 
# features    : 2904 
# extent      : -179.9572, 190.2908, -85.44308, 83.57391  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=stere +lat_0=90 +lon_0=-30 +x_0=0 +y_0=0 +ellps=WGS84 

# show the map with the 'unprojected' raster
plot_r <- spplot(r,
   sp.layout = list( sp.lines, mp, first = FALSE)
)

# show the map with the projected raster
plot_pr <- spplot(pr,
   sp.layout = list( sp.lines, mp, first = FALSE)   
)

print(plot_r)
print(plot_pr)
##### END

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

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages:
[1] rgdal_0.9-2     maptools_0.8-34 maps_2.3-9      raster_2.3-33   sp_1.0-17      

loaded via a namespace (and not attached):
[1] foreign_0.8-63  grid_3.1.2      lattice_0.20-30 tools_3.1.2 

Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org
moses selebatso | 27 Mar 08:44 2015
Picon

Help with latlong to UTM conversion when UTM zones are different

 Hello
I have animal movement data that I have converted from Lat/Long to UTM, unfortunately the data falls in two
UTM zones (34S and 35S). For some reason R cannot display both of them in the same window (the 35S data is way
off the expected location).
The question is how do I convert the data such that R can correctly read it?
Moses SELEBATSO

(+267) 318 5219 (H) (+267) 716 393 70 (C)
  (+267) 738 393 70 (C
	[[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo <at> r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Sara Rivero Calle | 26 Mar 16:58 2015
Picon
Picon

geom_points adjusting size for multiplot

Hi, 
 I am using ggplot2 and geometric points. I am trying to create a plot with 3 different panels, one for each
time period, I want the size of the points to be a function of the value in column “CB" and since there is a
lot of overlapping I used the alpha=CB so that its slightly transparent. CB goes from 0-420000. I used a
subset of the data frame for each time period and the problem is that each time period has a different
maximum (subset 1 ranges from 200 000, subset 2 has a max of 420 000 and subset 3 about 60 000. What I have so far
works but I want to have all 3 of them sharing the same legend and the size of the points be consistent for all 3
plots. How do I modify the breaks into 20 000, 50 000, 100 000, 200 000, 400 000 ?

library(maps)
library(ggplot2)
library(data frame)
library(spatstat)
library(sp)
library(map tools)

wrld <- map_data("world")
CPR<-read.csv("matchupP605.csv",header=TRUE)

subset1<-CPR[which(CPR$YearCollected>=1960 & CPR$YearCollected<1980& CPR$CB>0),]
 subset2 <-CPR[which(CPR$YearCollected>=1980 & CPR$YearCollected<2000& CPR$CB>0),]
 subset3 <-CPR[which(CPR$YearCollected>=2000 & CPR$YearCollected<2011& CPR$CB>0),]

source("http://peterhaschke.com/Code/multiplot.R") ## for the multiplot function.

	p6<-ggplot(subset1, aes(x = Longitude, y = Latitude))+ geom_point(aes(size=CB, alpha=CB),
colour="red")+ coord_equal() + geom_path(data=wrld, aes(x=long, y=lat, group = group))+
xlim(-80,10) + ylim(30,65) 	+ theme_bw(base_size = 12, base_family = "Helvetica")	# alpha modifies transparency

	p7<-ggplot(subset2, aes(x = Longitude, y = Latitude))+ geom_point(aes(size=CB, alpha=CB),
colour="red")+ coord_equal() + geom_path(data=wrld, aes(x=long, y=lat, group = group))+
xlim(-80,10) + ylim(30,65) + theme_bw(base_size = 12, base_family = "Helvetica")	

	p8<-ggplot(subset3, aes(x = Longitude, y = Latitude))+ geom_point(aes(size=CB, alpha=CB),
colour="red")+ coord_equal() + geom_path(data=wrld, aes(x=long, y=lat, group = group))+
xlim(-80,10) + ylim(30,65)+ theme_bw(base_size = 12, base_family = "Helvetica")

multiplot(p6, p7, p8)

Thanks in advance!
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo <at> r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Wolfgang Biener | 26 Mar 12:24 2015
Picon

dividing a line by points lying on it

Hi there,

I want to divide a line into n+1 lines. The shall be divided by points 
lying on it.
Since I don't have any idea how to solve can also post the starting 
point of the problem.

Thanks for your Help

library(sp)
Sl = SpatialLines(list(Lines(list(Line(cbind(c(1,2,2),c(1,1,2)))),
                              ID="a")))

set.seed(123)
sample.points <- spsample(Sl, 10, type="regular")

x11()
plot(Sl)
plot(sample.points, pch=1, col="red", add=T)

--

-- 
Wolfgang Biener,
Division Electrical Energy Systems
Fraunhofer-Institut für Solare Energiesysteme ISE
Heidenhofstrasse 2, 79110 Freiburg, Germany

Phone: +49 (0) 7 61/ 45 88-5893
wolfgang.biener <at> ise.fraunhofer.de
http://www.ise.fraunhofer.de

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo <at> r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Wolfgang | 26 Mar 12:23 2015
Picon

dividing a line by points lying on it

Hi there, 

I want to divide a line into n+1 lines. The shall be divided by points lying
on it. 
Since I don't have any idea how to solve can also post the starting point of
the problem. 

Thanks for your Help 

library(sp)
Sl = SpatialLines(list(Lines(list(Line(cbind(c(1,2,2),c(1,1,2)))),
                             ID="a")))

set.seed(123)
sample.points <- spsample(Sl, 10, type="regular")

x11()
plot(Sl)
plot(sample.points, pch=1, col="red", add=T)

--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/dividing-a-line-by-points-lying-on-it-tp7587958.html
Sent from the R-sig-geo mailing list archive at Nabble.com.

Gmane