[R-sig-Geo] raster - unrotate?
Michael Sumner
mdsumner at gmail.com
Thu Jul 20 02:47:06 CEST 2017
That's great, thanks for the update! I didn't actually realize all that
oceandata stuff was on Thredds. You might try ROpenSci's rerddap package
too.
I will try that source with the in-dev https://Github.com/hypertidy/tidync
. Please check it out if you can, and are feeling brave!
(tidync aims to provide more general outputs than raster can, more easily
than the API wrappers.)
Cheers, Mike
On Thu, 20 Jul 2017, 10:39 Ben Tupper, <btupper at bigelow.org> wrote:
> Ahoy!
>
> My (current) solution is to user raster::merge(). Works a treat and fits
> my workflow nicely. Since I am obtaining the raster data using OPeNDAP it
> is easy to get the two regions, rotate the extent of the eastern region,
> and then merge. Below is a pseudo-example (that works!) as it crops a
> local file twice, once for each region, rather than using an OPeNDaP calls
> via `ncdf4::nccvar_get(start = blah, count = something)` to get the two
> regions.
>
>
> ### start
>
> library(raster)
> uri <- '
> https://oceandata.sci.gsfc.nasa.gov:443/opendap/MODISA/L3SMI/2016/001/A20160012016032.L3m_R32_SST_sst_9km.nc
> '
> R <- raster::raster(uri, varname = 'sst')
>
> R
> #class : RasterLayer
> #dimensions : 180, 360, 64800 (nrow, ncol, ncell)
> #resolution : 1, 1 (x, y)
> #extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
> #coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> #data source : in memory
> #names : layer
> #values : 1.482572e-05, 0.9999928 (min, max)
>
> par(mfrow = c(2,1))
> plot(R)
> lines(c(180, 100, 100, 180), c(80,80,0,0))
> lines(c(-180,-90,-90,-180), c(80,80,0,0))
>
> # eastern section
> eR <- raster::crop(R, c(100, 180, 0, 80))
>
> # rotate the extent
> eExtent <- raster::extent(eR)
> raster::extent(eR) <- c(-260, -180, eExtent[3:4])
>
> # western section
> wR <- raster::crop(R, c(-180, -90, 0, 80))
>
> # merge
> newR <- raster::merge(eR, wR)
> newR
>
> plot(newR)
>
> #### end
>
> Cheers,
> Ben
>
> > On Jun 22, 2017, at 7:26 PM, Michael Sumner <mdsumner at gmail.com> wrote:
> >
> >
> >
> >
> > On Thu, 22 Jun 2017 at 21:28 Ben Tupper <btupper at bigelow.org> 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 <mdsumner at gmail.com> 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 <btupper at bigelow.org> wrote:
> >> Hello,
> >>
> >> We have rasters that span [-180, 180] from NASA's Ocean Color (
> https://oceancolor.gsfc.nasa.gov/) datasets. We are trying to extract a
> region that spans 100E to 90W, that is 100 to -90. The region 'wraps'
> across the edges as shown by the plot at the address below.
> >>
> >> https://dl.dropboxusercontent.com/u/8433654/sst.png
> >>
> >>
> >> library(raster)
> >> uri <- '
> https://oceandata.sci.gsfc.nasa.gov:443/opendap/MODISA/L3SMI/2016/001/A20160012016032.L3m_R32_SST_sst_9km.nc
> '
> >> R <- raster::raster(uri, varname = 'sst')
> >>
> >> R
> >> #class : RasterLayer
> >> #dimensions : 180, 360, 64800 (nrow, ncol, ncell)
> >> #resolution : 1, 1 (x, y)
> >> #extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
> >> #coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> >> #data source : in memory
> >> #names : layer
> >> #values : 1.482572e-05, 0.9999928 (min, max)
> >>
> >> plot(R)
> >> lines(c(180, 100, 100, 180), c(80,80,0,0))
> >> lines(c(-180,-90,-90,-180), c(80,80,0,0))
> >>
> >> I see that there is the nice raster::rotate() function to rotate raster
> coordinates from [0,360] to [-180,180]. That would make extracting the
> region super easy if the inverse were available. Is there an equivalent
> way to implement the inverse or raster::rotate()? I could dig into the
> source of raster::rotate() to see how to code one up, but I hate like heck
> to reinvent the wheel.
> >>
> >> Thanks!
> >> Ben
> >>
> >>
> >> Ben Tupper
> >> Bigelow Laboratory for Ocean Sciences
> >> 60 Bigelow Drive, P.O. Box 380
> >> East Boothbay, Maine 04544
> >> http://www.bigelow.org
> >>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> R-sig-Geo at r-project.org
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >> --
> >> Dr. Michael Sumner
> >> Software and Database Engineer
> >> Australian Antarctic Division
> >> 203 Channel Highway
> >> Kingston Tasmania 7050 Australia
> >>
> >
> > Ben Tupper
> > Bigelow Laboratory for Ocean Sciences
> > 60 Bigelow Drive, P.O. Box 380
> > East Boothbay, Maine 04544
> > http://www.bigelow.org
> >
> >
> >
> > --
> > Dr. Michael Sumner
> > Software and Database Engineer
> > Australian Antarctic Division
> > 203 Channel Highway
> > Kingston Tasmania 7050 Australia
>
> Ben Tupper
> Bigelow Laboratory for Ocean Sciences
> 60 Bigelow Drive, P.O. Box 380
> East Boothbay, Maine 04544
> http://www.bigelow.org
>
> Ecocast Reports: http://seascapemodeling.org/ecocast.html
> Tick Reports: https://report.bigelow.org/tick/
> Jellyfish Reports: https://jellyfish.bigelow.org/jellyfish/
>
>
>
> --
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list