Re: [R-sig-Geo] raster - unrotate?

2017-06-26 Thread Ben Tupper
Hi,

Thanks for the offer to help.  I have been redirected for the time being to 
other tasks, but I expect I'll be back to this in a few weeks.  I think I will 
try your archived rotate() as I'm starting with whole-earth coverage - nothing 
tricky.

You have been going really deep with tabularaster.  I have been folding 
tidyverse into my 'stuff' over the last few months, but I am on a relatively 
short leash resource-wise so the effort feels feeble.  Most of my work uses 
rasters and dismo https://cran.r-project.org/web/packages/dismo/index.html with 
only minor use of non-raster spatial data.  That might be my lame excuse for 
not diving into sf yet; that and I haven't seen examples I can ape for getting 
sf and raster to play together (well, until tabularaster). 

Thanks!
Ben

> On Jun 22, 2017, at 7:26 PM, Michael Sumner  wrote:
> 
> 
> 
> 
> On Thu, 22 Jun 2017 at 21:28 Ben Tupper  > wrote:
> Hi,
> 
> Wow!  A treasure from the past!
> 
> Hmmm.   I see what you mean; it might be best to snip-clip-and-smush from the 
> native presentation.  We are testing out the performance of the dineof 
> function in the sinkr package 
> (https://github.com/marchtaylor/sinkr/blob/master/R/dineof.R 
> ) to interpolate 
> patchy chlorophyl data.  
> 
> Thanks for the tip!
> 
> 
> 
> For what it's worth, I'm happy to help, there's no one size fits all here. 
> The dateline is one of those "easy problems" that does not have a general fix 
> and will bite in many different contexts. :)
> 
> The best distillation of my tricks is in this project, but it does depend on 
> your workflow and the actual data in use. 
> 
> https://github.com/hypertidy/tabularaster 
> 
> 
> If you have performance issues raster's cell abstraction  will help, but they 
> are a bit old-school when it comes to multi-part geometries. (It's identical 
> to standard database indexing concepts as are all "spatial" optimization 
> tricks)
> 
> (I see a desperate need for a general API for this kind of problem so that we 
> can build tools from a general framework, which I'm working on - that's some 
> kind of excuse for why this and related projects are quite raw and 
> unfinished.)
> 
> Cheers, Mike.  
> 
> Cheers,
> Ben
> 
>> On Jun 22, 2017, at 4:46 AM, Michael Sumner > > wrote:
>> 
>> 
>> It used to do the inverse, and I preserved a copy of the old behaviour here: 
>> 
>> https://github.com/AustralianAntarcticDivision/raadtools/blob/master/R/utils.R#L9
>>  
>> 
>> 
>> You're probably best to learn from how it approaches it rather than just use 
>> it, since it's not a general capability, just works for those very specific 
>> cases. 
>> 
>> I just tested it and it still seems to work fine, I used 
>> 
>> oceandata.sci.gsfc.nasa.gov/MODISA/Mapped/Monthly/9km/chlor/A20021822002212.L3m_MO_CHL_chlor_a_9km.nc
>>  
>> 
>>  
>> 
>> obtainable with 
>> 
>> https://oceandata.sci.gsfc.nasa.gov/cgi/getfile/A20021822002212.L3m_MO_CHL_chlor_a_9km.nc
>>  
>> 
>>  
>> 
>> If you're extracting points you might as well just convert them into the 
>> "Atlantic" range, but it's painful to get that right for polygons or lines 
>> (and very hard to generalize once you get it working anyway). 
>> 
>> For simple regions like the one you show I'd probably split it into two 
>> extents() and use those - but depends on your workflow, all these options 
>> are useful here and there. 
>> 
>> Cheers, Mike. 
>> 
>> 
>> On Thu, 22 Jun 2017 at 18:30 Ben Tupper > > wrote:
>> Hello,
>> 
>> We have rasters that span [-180, 180] from NASA's Ocean Color 
>> (https://oceancolor.gsfc.nasa.gov/ )  
>> datasets.  We are trying to extract a region that spans 100E to 90W, that is 
>> 100 to -90.  The region 'wraps' across the edges as shown by the plot at the 
>> address below.
>> 
>> https://dl.dropboxusercontent.com/u/8433654/sst.png 
>> 
>> 
>> 
>> library(raster)
>> uri <- 
>> 'https://oceandata.sci.gsfc.nasa.gov:443/opendap/MODISA/L3SMI/2016/001/A20160012016032.L3m_R32_SST_sst_9km.nc
>>  
>> '
>> R <- raster::raster(uri, varname = 'sst')
>> 
>> R
>> #class   : RasterLayer
>> #dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
>> #resolution  : 1, 1  (x, y)
>> #extent  : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
>> #coord. ref. : 

Re: [R-sig-Geo] Using gIntersection to split polygons with lines

2017-06-26 Thread Edzer Pebesma


On 26/06/17 12:56, Frederico Mestre wrote:
> Hello,
> 
> I'm trying to use gIntersection() to cut a polygon ("SpatialPolygons") with
> a line ("SpatialLinesDataFrame").
> 
> I wanted to obtain a "SpatialPolygons" object with the polygons defined by
> the lines. This is what I did:
> 
> polygons1 <- gIntersection(spgeom1=mypolygon, spgeom2=mylines)
> 
> However I get an "SpatialLines" object as an output instead of a
> "SpatialPolygons".
> 
> Am I doing something wrong? Is there any alternative?

Package sf has a function st_split() that cuts a polygon with a line,
similar to postgis' function st_split.

To use it, you'd need to build sf from source after having liblwgeom
installed; the CRAN binary builds do not have support for liblwgeom.

See https://github.com/edzer/sfr/ for instructions, and
https://github.com/edzer/sfr/issues/360 for an example.

> 
> I provide this sample code:
> 
> #Creating the polygons
> coords <- rbind(c(1,2),c(2,2),c(2,1),c(1,1))
> pol1 <- Polygon(coords)
> pol2 <- Polygons(srl=list(pol1), ID=1)
> pol3 <- SpatialPolygons(Srl=list(pol2))
> 
> #Creating the line
> l1 <- Line(rbind(c(0.5,1.5),c(2.5,1.5)))
> l2 <- Lines(list(l1),ID=1)
> l3 <- SpatialLines(list(l2))
> line1 <- SpatialLinesDataFrame(l3,data=as.data.frame(matrix(ncol=2,nrow=1)))
> 
> #Ploting both
> plot(pol3)
> plot(line1, add=TRUE)
> 
> #Intersection
> inters1 <- gIntersection(spgeom1=pol3, spgeom2=line1)
> inters1 <- gIntersection(spgeom1=line1, spgeom2=pol3)
> plot(inters1)
> class(inters1)
> 
> It only shows the lines (the object is of class "SpatialLines").
> 
> Thanks,
> Frederico
> ᐧ
> 
>   [[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  (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/



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

[R-sig-Geo] Using gIntersection to split polygons with lines

2017-06-26 Thread Frederico Mestre
Hello,

I'm trying to use gIntersection() to cut a polygon ("SpatialPolygons") with
a line ("SpatialLinesDataFrame").

I wanted to obtain a "SpatialPolygons" object with the polygons defined by
the lines. This is what I did:

polygons1 <- gIntersection(spgeom1=mypolygon, spgeom2=mylines)

However I get an "SpatialLines" object as an output instead of a
"SpatialPolygons".

Am I doing something wrong? Is there any alternative?

I provide this sample code:

#Creating the polygons
coords <- rbind(c(1,2),c(2,2),c(2,1),c(1,1))
pol1 <- Polygon(coords)
pol2 <- Polygons(srl=list(pol1), ID=1)
pol3 <- SpatialPolygons(Srl=list(pol2))

#Creating the line
l1 <- Line(rbind(c(0.5,1.5),c(2.5,1.5)))
l2 <- Lines(list(l1),ID=1)
l3 <- SpatialLines(list(l2))
line1 <- SpatialLinesDataFrame(l3,data=as.data.frame(matrix(ncol=2,nrow=1)))

#Ploting both
plot(pol3)
plot(line1, add=TRUE)

#Intersection
inters1 <- gIntersection(spgeom1=pol3, spgeom2=line1)
inters1 <- gIntersection(spgeom1=line1, spgeom2=pol3)
plot(inters1)
class(inters1)

It only shows the lines (the object is of class "SpatialLines").

Thanks,
Frederico
ᐧ

[[alternative HTML version deleted]]

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