7

I'm quite okay at R programming, but fairly new in GIS (especially in R).

I have a set of coordinates and I'd like to create a custom shape around them. What I want to do step-by-step:

  1. Create a circle around the points and slice it to couple of sectors. The radius and the number of wedges are given manually at this point.
  2. Group the points by the wedges
  3. Calculate the radius of the given wedge based on the points characteristic of the points in the sector (I have the model for this)
  4. Set the radius sector by sector.
  5. Export the wedges as different regions of a polygon shapefile

enter image description here

Final results should look like this

Unfortunately, I'm stuck at the very beginning. I created a standard pie chart by using add.pie {mapplots}, but I wasn't able to convert its sectors to spatial objects.

  • 1
    Quite a few steps here. Break it down. First step might be "compute the coordinates of a wedge of a circle centred at (x,y), radius r, from angle theta1 to theta2 and create a SpatialPolygon object". Write the function definition and then we might be able to fill it in. That would be a whole question here - we like one Q and one A on SO... – Spacedman Oct 4 '15 at 18:46
  • @Spacedman Your are right, I give it a try. ;) – user1731038 Oct 4 '15 at 18:48
  • An example of this kind of calculation, involving the creation of wedge shapes according to attributes of a table and their conversion to shapefile format, appears at gis.stackexchange.com/questions/31294. – whuber Oct 6 '15 at 20:14
  • Very nice indeed! – mdsumner Oct 22 '15 at 0:34
1

In the interests of bending tools in ways they weren't designed for, you can get there quickly with the maptools and graticule packages.

(You should really think about doing this with a proper map projection, and in some ways that would be simpler.)

library(graticule)

## say we want a circle at this longlat point and with a radius in metres of "distm"

llpt <- c(147, -42)
distm <- 10000

## decide what distm means for you
offsetLL <-  c(cos(llpt[2] * pi/180), 1) /  ((6370997 * pi)/180/ distm)

## circle segments very much like a graticule at the pole
g <- graticule(seq(-180, 180, by = 45), c(-90, -90 + offsetLL[2]), xlim = c(-180, 180), 
  proj = "+proj=laea +lat_0=-90 +ellps=sphere", tiles = TRUE)

## smash the object and shift/scale 
library(maptools)

g1 <- elide(g, scale = mean(offsetLL))

g2 <- elide(g1, shift = llpt - offsetLL/2)
proj4string(g2) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
  • 1
    That's a very smart solution! However, I still could get it work. I've created graticule and shifted its centroid to the desired coordinates, but as soon as I reproject the LAEA CRS to WGS84 the graticule become so destorted, Am I missing something here? Actually, I'm planning to use UTM anyway. I don't know if it's make it easier though. – user1731038 Oct 15 '15 at 12:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.