raster - unrotate?

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
7 messages Options
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

raster - unrotate?

Ben Tupper
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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: raster - unrotate?

Michael Sumner-2
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 <[hidden email]> 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
> [hidden email]
> 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

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: raster - unrotate?

Ben Tupper
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!

Cheers,
Ben

> On Jun 22, 2017, at 4:46 AM, Michael Sumner <[hidden email]> 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 <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 <http://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 <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 <[hidden email] <mailto:[hidden email]>> wrote:
> Hello,
>
> We have rasters that span [-180, 180] from NASA's Ocean Color (https://oceancolor.gsfc.nasa.gov/ <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 <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 <https://oceandata.sci.gsfc.nasa.gov/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 <http://www.bigelow.org/>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email] <mailto:[hidden email]>
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo <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




        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: raster - unrotate?

Michael Sumner-2
On Thu, 22 Jun 2017 at 21:28 Ben Tupper <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> 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
>> <https://oceandata.sci.gsfc.nasa.gov/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
>> [hidden email]
>> 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

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: raster - unrotate?

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 <[hidden email]> wrote:
>
>
>
>
> On Thu, 22 Jun 2017 at 21:28 Ben Tupper <[hidden email] <mailto:[hidden email]>> 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 <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 <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 <[hidden email] <mailto:[hidden email]>> 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 <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 <http://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 <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 <[hidden email] <mailto:[hidden email]>> wrote:
>> Hello,
>>
>> We have rasters that span [-180, 180] from NASA's Ocean Color (https://oceancolor.gsfc.nasa.gov/ <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 <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 <https://oceandata.sci.gsfc.nasa.gov/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 <http://www.bigelow.org/>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email] <mailto:[hidden email]>
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo <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 <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




        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: raster - unrotate?

Ben Tupper
In reply to this post by Michael Sumner-2
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 <[hidden email]> wrote:
>
>
>
>
> On Thu, 22 Jun 2017 at 21:28 Ben Tupper <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> 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
>> [hidden email]
>> 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/

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: raster - unrotate?

Michael Sumner-2
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, <[hidden email]> 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 <[hidden email]> wrote:
> >
> >
> >
> >
> > On Thu, 22 Jun 2017 at 21:28 Ben Tupper <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> 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
> >> [hidden email]
> >> 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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Loading...