Subscribe to R-sig-geo feed
This is an archive for R-sig-geo, not a forum. It is not possible to post through Nabble - you may not start a new thread nor follow up an existing thread. If you wish to post, but are not subscribed to the list through the list homepage, subscribe first (through the list homepage, not via Nabble) and post from your subscribed email address. Until 2015-06-20, subscribers could post through Nabble, but policy has been changed as too many non-subscribers misunderstood the interface.
Updated: 2 hours 2 min ago

Re: Easiest way to get a a subset from MODIS

Wed, 05/24/2017 - 23:34
Answering my own question.

I was able to do what I wanted using the MODIS package. The runGdal function does what I want. Thanks to those who coded it!

-A

From: Andy Bunn <[hidden email]<mailto:[hidden email]>>
Date: Wednesday, May 24, 2017 at 12:22 PM
To: R-sig-Geo <[hidden email]<mailto:[hidden email]>>
Subject: Easiest way to get a a subset from MODIS

Hi all, It's been years since I've used any MODIS data. Looking around I see there are many R packages that purport to interface nicely with MODIS servers and deliver you a subset of data.

I have an EarthData login and can (and have) download the tile I want via http://search.earthdata.nasa.gov but getting the hdf, converting to tif, clipping etc. is tedious and I'm wondering what the best way to do this in R is. So,  let's imagine I would like to get a class(raster) object of the "250m_16_days_NDVI" band from "MYD13Q1" product for a 10x10km area around the point c(68.75,161.42) for the period of nearest July 1, 2015. What package and approach is the easiest in R?

Is it something like this?

library(MODISTools)
dat2get <- data.frame(lat=68.75,long=161.42,
                      start.date=as.POSIXlt("2015-07-01"),
                      end.date=as.POSIXlt("2015-07-01"),
                      ID=1)

MODISSubsets(LoadDat = dat2get,
             Products = "MYD13Q1",
             Bands = "250m_16_days_NDVI",
             Size = c(10,10),
             StartDate = TRUE)

Regards, Andy



        [[alternative HTML version deleted]]

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

Re: Deal with multiple factorlevel in one grid square

Wed, 05/24/2017 - 19:28
Hiya,

Do you want a raster with attribute information (ie several attributes from your  shape file)? You can achieve this with a workflow that uses
raster::extract() to assigne unique IDs of polygon to raster cells
create a RAT for you raster
assign to the unique IDs of raster RAT the attributes of your polygon unique IDs (many to one relationship) with something like merge() or functions in libraries like diplyr or data.table

Cheers
Herry

-----Original Message-----
From: R-sig-Geo [mailto:[hidden email]] On Behalf Of Patrick Schratz
Sent: Wednesday, 24 May 2017 10:48 PM
To: Michael Sumner <[hidden email]>; Miriam Püts <[hidden email]>
Cc: [hidden email]
Subject: Re: [R-sig-Geo] Deal with multiple factorlevel in one grid square

Do you want to do what is called “zonal statistics” in ArcGIS with “majority" option?

You may check out this SO question and try the “mode” function within raster::extract() - maybe it does what you need.

(I’m also unsure if I understand the question correctly)

Cheers, Pat

PhD Student at Department of Geography - GIScience group Friedrich-Schiller-University Jena, Germany
Tel.: +49-3641-9-48973
Web: https://pat-s.github.io

On 24. May 2017, 12:35 +0200, Miriam Püts <[hidden email]>, wrote:
> Hi Mike,
>
> I will try to explain it a bit more in detail. Maybe it is easier to understand if I start from the end. In the end I would like to have an ASCII file to read into Ecospace, which has the same extensions and coordinates as other files I already created. This ASCII should contain information about the sediment type within each predefined cell. To create this ASCII file I have a shape file with the polygons representing the sediment type and my grid which I applied to other variables to have the same extend. Now I would like to create a grid containing the information on sediment. Here, per grid cell the sediment type which covers the most of the cell should be defined and connected with the coordinates for the grid cell.
>
> I hope this makes it more clear...
>
>
> <°)))>< >°)))>< >°)))>< >°)))><
>
> Miriam Püts
> Marine Lebende Ressourcen/ Marine Living Resources Thünen-Institut für
> Seefischerei/ Thünen Institute of Sea fisheries Palmaille 9
> 22767 Hamburg (Germany)
>
> Tel: +49 40 38905-105
> Mail: [hidden email]
>
>
> Von: "Michael Sumner" <[hidden email]
> An: "Miriam Püts" <[hidden email]>, [hidden email]
> Gesendet: Mittwoch, 24. Mai 2017 11:55:14
> Betreff: Re: [R-sig-Geo] Deal with multiple factorlevel in one grid
> square
>
>
>
> raster::extract(grid, poly, weights = TRUE) is a start. It returns a list which is painful to deal with at first, but can be collected into one data frame for standard summarizing.
>
> I'm still a bit confused about whether you want an estimate of a cell overlap in a polygon or something else.
>
> Cheers, Mike
> On Wed, 24 May 2017, 17:12 Miriam Püts < [ mailto:[hidden email] | [hidden email] ] > wrote:
>
>
> Hi everyone,
>
> I have the following problem: I have a personalized grid and a shape file with polygons representing sediment types. Now I would like to apply this grid to the Polygons to identify the sediment type most common within each grid. I tried it with rasterize, but here I can only chose last or first. Die you have any suggestions how I might get the sediment type for each raster cell?
>
> Thank you for your help!
>
> <°)))>< >°)))>< >°)))>< >°)))><
>
> Miriam Püts
> Marine Lebende Ressourcen/ Marine Living Resources Thünen-Institut für
> Seefischerei/ Thünen Institute of Sea fisheries Palmaille 9
> 22767 Hamburg (Germany)
>
> Tel: +49 40 38905-105
> Mail: [ mailto:[hidden email] | [hidden email] ]
>
> _______________________________________________
> R-sig-Geo mailing list
> [ mailto:[hidden email] | [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
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
        [[alternative HTML version deleted]]

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

Re: Recs for R tools that subset/tile raster bricks?

Wed, 05/24/2017 - 18:18
Thanks for your response, Mike! The files are ENVI images (*.dat). The
output of print() on one is copied below. The other images have slightly
different row/col sizes but the same number of layers.

Thanks for the tip on the vignettes. It looks like "Writing functions with
the ”raster” package" is the one you're referring to? blockSize() looks
helpful as I'm not sure what size would be appropriate to crop to.

--

class       : RasterBrick
dimensions  : 24390, 2032, 49560480, 373  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : 520554.6, 522586.6, 4766639, 4791029  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=11 +datum=NAD83 +units=m +no_defs
+ellps=GRS80 +towgs84=0,0,0
data source : E:\RCEW_CoReg_2014\183421\183421_Georegistered
names       :
Warp..Band.8.ang20140919t183421_corr_v1c_img...381.360260.Nanometers.,
Warp..Band.9.ang20140919t183421_corr_v1c_img...386.368927.Nanometers.,
Warp..Band.10.ang20140919t183421_corr_v1c_img...391.377594.Nanometers.,
Warp..Band.11.ang20140919t183421_corr_v1c_img...396.386292.Nanometers.,
Warp..Band.12.ang20140919t183421_corr_v1c_img...401.394958.Nanometers.,
Warp..Band.13.ang20140919t183421_corr_v1c_img...406.403625.Nanometers.,
Warp..Band.14.ang20140919t183421_corr_v1c_img...411.412292.Nanometers.,
Warp..Band.15.ang20140919t183421_corr_v1c_img...416.420959.Nanometers.,
Warp..Band.16.ang20140919t183421_corr_v1c_img...421.429626.Nanometers.,
Warp..Band.17.ang20140919t183421_corr_v1c_img...426.438293.Nanometers.,
Warp..Band.18.ang20140919t183421_corr_v1c_img...431.446960.Nanometers.,
Warp..Band.19.ang20140919t183421_corr_v1c_img...436.455627.Nanometers.,
Warp..Band.20.ang20140919t183421_corr_v1c_img...441.464294.Nanometers.,
Warp..Band.21.ang20140919t183421_corr_v1c_img...446.... <truncated>

Thanks,
Megan

On Wed, May 24, 2017 at 6:58 PM, Michael Sumner <[hidden email]> wrote:

> Do also see the raster vignettes, there is one on performance and using
> line- or block- based chunks for large jobs.
>
> Cheers, Mike
>
>
> On Thu, 25 May 2017, 08:32 Michael Sumner <[hidden email]> wrote:
>
>> Try crop() (or more abstractly extract(raster, cells) once you've
>> determined which cells x-y are needed).
>>
>> Getting the most efficient way to do this depends on the source of the
>> brick. It's not possible to know without knowing details on your setup. The
>> file format, method of creation and dimensions of the brick and sizes of
>> your subset query can all affect the best ways to go.
>>
>> Happy to help but it might take a bit of conversation to figure out
>> what's best. If you can share the output of print on the brick anf the
>> extent of a subset that is a good start.
>>
>> Cheers, Mike
>>
>> On Thu, 25 May 2017, 05:50 Megan Maloney <[hidden email]> wrote:
>>
>>> Hello all,
>>>
>>> Does anyone know of a good R tool to subset/tile raster bricks?
>>>
>>> I'd like to break a raster brick into smaller pieces to process. I'm
>>> finding functions to subset or tile a single raster layer, but not for
>>> multi-layer raster bricks.
>>>
>>> The brick input is hyperspectral flightlines - very large, 1 m^2
>>> resampling
>>> covering a watershed, 340 layers. My code currently takes in a raster
>>> brick, extracts values per cell, applies beta coefficients, and outputs a
>>> nitrogen value to a raster.
>>>
>>> I don't need the subsets/tiles to be a particular size or number, just
>>> smaller. I was warned against dropping a whole flightline into R. I'm
>>> new to
>>> R, so I'm not sure what limitations I'm running into except processing
>>> time.
>>> All advice is welcome.
>>>
>>> Thank you very much for your help.
>>> Megan
>>>
>>>         [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> 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
>>
>> --
> 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

Re: Recs for R tools that subset/tile raster bricks?

Wed, 05/24/2017 - 17:58
Do also see the raster vignettes, there is one on performance and using
line- or block- based chunks for large jobs.

Cheers, Mike

On Thu, 25 May 2017, 08:32 Michael Sumner <[hidden email]> wrote:

> Try crop() (or more abstractly extract(raster, cells) once you've
> determined which cells x-y are needed).
>
> Getting the most efficient way to do this depends on the source of the
> brick. It's not possible to know without knowing details on your setup. The
> file format, method of creation and dimensions of the brick and sizes of
> your subset query can all affect the best ways to go.
>
> Happy to help but it might take a bit of conversation to figure out what's
> best. If you can share the output of print on the brick anf the extent of a
> subset that is a good start.
>
> Cheers, Mike
>
> On Thu, 25 May 2017, 05:50 Megan Maloney <[hidden email]> wrote:
>
>> Hello all,
>>
>> Does anyone know of a good R tool to subset/tile raster bricks?
>>
>> I'd like to break a raster brick into smaller pieces to process. I'm
>> finding functions to subset or tile a single raster layer, but not for
>> multi-layer raster bricks.
>>
>> The brick input is hyperspectral flightlines - very large, 1 m^2
>> resampling
>> covering a watershed, 340 layers. My code currently takes in a raster
>> brick, extracts values per cell, applies beta coefficients, and outputs a
>> nitrogen value to a raster.
>>
>> I don't need the subsets/tiles to be a particular size or number, just
>> smaller. I was warned against dropping a whole flightline into R. I'm new
>> to
>> R, so I'm not sure what limitations I'm running into except processing
>> time.
>> All advice is welcome.
>>
>> Thank you very much for your help.
>> Megan
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> 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
>
> -- 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

Re: Recs for R tools that subset/tile raster bricks?

Wed, 05/24/2017 - 17:33
Try crop() (or more abstractly extract(raster, cells) once you've
determined which cells x-y are needed).

Getting the most efficient way to do this depends on the source of the
brick. It's not possible to know without knowing details on your setup. The
file format, method of creation and dimensions of the brick and sizes of
your subset query can all affect the best ways to go.

Happy to help but it might take a bit of conversation to figure out what's
best. If you can share the output of print on the brick anf the extent of a
subset that is a good start.

Cheers, Mike

On Thu, 25 May 2017, 05:50 Megan Maloney <[hidden email]> wrote:

> Hello all,
>
> Does anyone know of a good R tool to subset/tile raster bricks?
>
> I'd like to break a raster brick into smaller pieces to process. I'm
> finding functions to subset or tile a single raster layer, but not for
> multi-layer raster bricks.
>
> The brick input is hyperspectral flightlines - very large, 1 m^2 resampling
> covering a watershed, 340 layers. My code currently takes in a raster
> brick, extracts values per cell, applies beta coefficients, and outputs a
> nitrogen value to a raster.
>
> I don't need the subsets/tiles to be a particular size or number, just
> smaller. I was warned against dropping a whole flightline into R. I'm new
> to
> R, so I'm not sure what limitations I'm running into except processing
> time.
> All advice is welcome.
>
> Thank you very much for your help.
> Megan
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> 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

Recs for R tools that subset/tile raster bricks?

Wed, 05/24/2017 - 14:49
Hello all,

Does anyone know of a good R tool to subset/tile raster bricks?

I'd like to break a raster brick into smaller pieces to process. I'm
finding functions to subset or tile a single raster layer, but not for
multi-layer raster bricks.

The brick input is hyperspectral flightlines - very large, 1 m^2 resampling
covering a watershed, 340 layers. My code currently takes in a raster
brick, extracts values per cell, applies beta coefficients, and outputs a
nitrogen value to a raster.

I don't need the subsets/tiles to be a particular size or number, just
smaller. I was warned against dropping a whole flightline into R. I'm new to
R, so I'm not sure what limitations I'm running into except processing time.
All advice is welcome.

Thank you very much for your help.
Megan

        [[alternative HTML version deleted]]

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

Easiest way to get a a subset from MODIS

Wed, 05/24/2017 - 14:22
Hi all, It's been years since I've used any MODIS data. Looking around I see there are many R packages that purport to interface nicely with MODIS servers and deliver you a subset of data.

I have an EarthData login and can (and have) download the tile I want via http://search.earthdata.nasa.gov but getting the hdf, converting to tif, clipping etc. is tedious and I'm wondering what the best way to do this in R is. So,  let's imagine I would like to get a class(raster) object of the "250m_16_days_NDVI" band from "MYD13Q1" product for a 10x10km area around the point c(68.75,161.42) for the period of nearest July 1, 2015. What package and approach is the easiest in R?

Is it something like this?

library(MODISTools)
dat2get <- data.frame(lat=68.75,long=161.42,
                      start.date=as.POSIXlt("2015-07-01"),
                      end.date=as.POSIXlt("2015-07-01"),
                      ID=1)

MODISSubsets(LoadDat = dat2get,
             Products = "MYD13Q1",
             Bands = "250m_16_days_NDVI",
             Size = c(10,10),
             StartDate = TRUE)

Regards, Andy



        [[alternative HTML version deleted]]

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

Re: gstat or geoR create variogram with sets of bounding boxes

Wed, 05/24/2017 - 08:35
THAT WORK FINE,
THANKS!

2017-05-24 7:17 GMT-03:00 Dr. Benedikt Gräler <[hidden email]>:

> Dear Javier,
>
> if I got you right, you seem to be looking for a "pooled estimation of
> within-strata variograms" as the help page on the function "variogram()" of
> gstat calls it. At first, define an indicator for each of your points that
> identifies the ID of the polygon it belongs to, then include it as
> regressor and set dX=0.5:
>
> Try the example below for the meuse data set where only point-pairs of the
> same soil type (corresponding to your polygon Ids) are considered and then
> a pooled variogram is calculated.
>
> HTH,
>
>  Ben
>
>
> library(sp)
> library(gstat)
>
> data(meuse)
> coordinates(meuse) <- ~x+y
>
> spplot(meuse, "zinc")
> spplot(meuse, "soil")
>
> vAll <- variogram(zinc~soil, meuse)
> plot(vAll)
>
> vWithinSoils <- variogram(zinc~soil, meuse, dX=0.5)
> plot(vWithinSoils)
>
>
>
>
> On 23/05/2017 17:23, Javier Moreira wrote:
>
>> hi,
>> im trying to make a variogram that consider a limit made by a bounding
>> box,
>> but not just one, but dividing the area of the hole set of points in
>> different subzones (treatments in this case).
>> I have tried to do so with bounding box in geoR, but i cant do it.
>> Also, another way its to use anisotropy, but when try to put an alpha
>> parameter in geostat, doesnt make any change that without it.
>> bottom line, i have to perform a variogram that includes only the points
>> in
>> 1 direction, or, tell the variogram to limit the surface to a polygon (
>> several ones)
>>
>> if any one can help,
>> thanks
>>
>>
> --
> Dr. Benedikt Gräler
> 52°North Initiative for Geospatial Open Source Software GmbH
> Martin-Luther-King-Weg 24
> 48155 Muenster, Germany
>
> E-Mail: [hidden email]
> Fon: +49-(0)-251/396371-39
> Fax: +49-(0)-251/396371-11
>
> http://52north.org/
> Twitter: @FiveTwoN
>
> General Managers: Dr. Albert Remke, Dr. Andreas Wytzisk
> Local Court Muenster HRB 10849
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>


--
Javier Moreira de Souza
Ingeniero Agrónomo
099 406 006

        [[alternative HTML version deleted]]

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

Re: Deal with multiple factorlevel in one grid square

Wed, 05/24/2017 - 07:48
Do you want to do what is called “zonal statistics” in ArcGIS with “majority" option?

You may check out this SO question and try the “mode” function within raster::extract() - maybe it does what you need.

(I’m also unsure if I understand the question correctly)

Cheers, Pat

PhD Student at Department of Geography - GIScience group
Friedrich-Schiller-University Jena, Germany
Tel.: +49-3641-9-48973
Web: https://pat-s.github.io

On 24. May 2017, 12:35 +0200, Miriam Püts <[hidden email]>, wrote:
> Hi Mike,
>
> I will try to explain it a bit more in detail. Maybe it is easier to understand if I start from the end. In the end I would like to have an ASCII file to read into Ecospace, which has the same extensions and coordinates as other files I already created. This ASCII should contain information about the sediment type within each predefined cell. To create this ASCII file I have a shape file with the polygons representing the sediment type and my grid which I applied to other variables to have the same extend. Now I would like to create a grid containing the information on sediment. Here, per grid cell the sediment type which covers the most of the cell should be defined and connected with the coordinates for the grid cell.
>
> I hope this makes it more clear...
>
>
> <°)))>< >°)))>< >°)))>< >°)))><
>
> Miriam Püts
> Marine Lebende Ressourcen/ Marine Living Resources
> Thünen-Institut für Seefischerei/ Thünen Institute of Sea fisheries
> Palmaille 9
> 22767 Hamburg (Germany)
>
> Tel: +49 40 38905-105
> Mail: [hidden email]
>
>
> Von: "Michael Sumner" <[hidden email]
> An: "Miriam Püts" <[hidden email]>, [hidden email]
> Gesendet: Mittwoch, 24. Mai 2017 11:55:14
> Betreff: Re: [R-sig-Geo] Deal with multiple factorlevel in one grid square
>
>
>
> raster::extract(grid, poly, weights = TRUE) is a start. It returns a list which is painful to deal with at first, but can be collected into one data frame for standard summarizing.
>
> I'm still a bit confused about whether you want an estimate of a cell overlap in a polygon or something else.
>
> Cheers, Mike
> On Wed, 24 May 2017, 17:12 Miriam Püts < [ mailto:[hidden email] | [hidden email] ] > wrote:
>
>
> Hi everyone,
>
> I have the following problem: I have a personalized grid and a shape file with polygons representing sediment types. Now I would like to apply this grid to the Polygons to identify the sediment type most common within each grid. I tried it with rasterize, but here I can only chose last or first. Die you have any suggestions how I might get the sediment type for each raster cell?
>
> Thank you for your help!
>
> <°)))>< >°)))>< >°)))>< >°)))><
>
> Miriam Püts
> Marine Lebende Ressourcen/ Marine Living Resources
> Thünen-Institut für Seefischerei/ Thünen Institute of Sea fisheries
> Palmaille 9
> 22767 Hamburg (Germany)
>
> Tel: +49 40 38905-105
> Mail: [ mailto:[hidden email] | [hidden email] ]
>
> _______________________________________________
> R-sig-Geo mailing list
> [ mailto:[hidden email] | [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
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
        [[alternative HTML version deleted]]

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

Re: Deal with multiple factorlevel in one grid square

Wed, 05/24/2017 - 05:35
Hi Mike,

I will try to explain it a bit more in detail. Maybe it is easier to understand if I start from the end. In the end I would like to have an ASCII file to read into Ecospace, which has the same extensions and coordinates as other files I already created. This ASCII should contain information about the sediment type within each predefined cell. To create this ASCII file I have a shape file with the polygons representing the sediment type and my grid which I applied to other variables to have the same extend. Now I would like to create a grid containing the information on sediment. Here, per grid cell the sediment type which covers the most of the cell should be defined and connected with the coordinates for the grid cell.

I hope this makes it more clear...


<°)))>< >°)))>< >°)))>< >°)))><

Miriam Püts
Marine Lebende Ressourcen/ Marine Living Resources
Thünen-Institut für Seefischerei/ Thünen Institute of Sea fisheries
Palmaille 9
22767 Hamburg (Germany)

Tel: +49 40 38905-105
Mail: [hidden email]


Von: "Michael Sumner" <[hidden email]>
An: "Miriam Püts" <[hidden email]>, [hidden email]
Gesendet: Mittwoch, 24. Mai 2017 11:55:14
Betreff: Re: [R-sig-Geo] Deal with multiple factorlevel in one grid square



raster::extract(grid, poly, weights = TRUE) is a start. It returns a list which is painful to deal with at first, but can be collected into one data frame for standard summarizing.

I'm still a bit confused about whether you want an estimate of a cell overlap in a polygon or something else.

Cheers, Mike
On Wed, 24 May 2017, 17:12 Miriam Püts < [ mailto:[hidden email] | [hidden email] ] > wrote:


Hi everyone,

I have the following problem: I have a personalized grid and a shape file with polygons representing sediment types. Now I would like to apply this grid to the Polygons to identify the sediment type most common within each grid. I tried it with rasterize, but here I can only chose last or first. Die you have any suggestions how I might get the sediment type for each raster cell?

Thank you for your help!

<°)))>< >°)))>< >°)))>< >°)))><

Miriam Püts
Marine Lebende Ressourcen/ Marine Living Resources
Thünen-Institut für Seefischerei/ Thünen Institute of Sea fisheries
Palmaille 9
22767 Hamburg (Germany)

Tel: +49 40 38905-105
Mail: [ mailto:[hidden email] | [hidden email] ]

_______________________________________________
R-sig-Geo mailing list
[ mailto:[hidden email] | [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


        [[alternative HTML version deleted]]

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

Re: gstat or geoR create variogram with sets of bounding boxes

Wed, 05/24/2017 - 05:17
Dear Javier,

if I got you right, you seem to be looking for a "pooled estimation of
within-strata variograms" as the help page on the function "variogram()"
of gstat calls it. At first, define an indicator for each of your points
that identifies the ID of the polygon it belongs to, then include it as
regressor and set dX=0.5:

Try the example below for the meuse data set where only point-pairs of
the same soil type (corresponding to your polygon Ids) are considered
and then a pooled variogram is calculated.

HTH,

  Ben


library(sp)
library(gstat)

data(meuse)
coordinates(meuse) <- ~x+y

spplot(meuse, "zinc")
spplot(meuse, "soil")

vAll <- variogram(zinc~soil, meuse)
plot(vAll)

vWithinSoils <- variogram(zinc~soil, meuse, dX=0.5)
plot(vWithinSoils)




On 23/05/2017 17:23, Javier Moreira wrote:
> hi,
> im trying to make a variogram that consider a limit made by a bounding box,
> but not just one, but dividing the area of the hole set of points in
> different subzones (treatments in this case).
> I have tried to do so with bounding box in geoR, but i cant do it.
> Also, another way its to use anisotropy, but when try to put an alpha
> parameter in geostat, doesnt make any change that without it.
> bottom line, i have to perform a variogram that includes only the points in
> 1 direction, or, tell the variogram to limit the surface to a polygon (
> several ones)
>
> if any one can help,
> thanks
> --
Dr. Benedikt Gräler
52°North Initiative for Geospatial Open Source Software GmbH
Martin-Luther-King-Weg 24
48155 Muenster, Germany

E-Mail: [hidden email]
Fon: +49-(0)-251/396371-39
Fax: +49-(0)-251/396371-11

http://52north.org/
Twitter: @FiveTwoN

General Managers: Dr. Albert Remke, Dr. Andreas Wytzisk
Local Court Muenster HRB 10849

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
b_graeler.vcf (435 bytes) Download Attachment

Re: Deal with multiple factorlevel in one grid square

Wed, 05/24/2017 - 04:55
raster::extract(grid, poly, weights = TRUE) is a start. It returns a list
which is painful to deal with at first, but can be collected into one data
frame for standard summarizing.

I'm still a bit confused about whether you want an estimate of a cell
overlap in a polygon or something else.

Cheers, Mike

On Wed, 24 May 2017, 17:12 Miriam Püts <[hidden email]> wrote:

> Hi everyone,
>
> I have the following problem: I have a personalized grid and a shape file
> with polygons representing sediment types. Now I would like to apply this
> grid to the Polygons to identify the sediment type most common within each
> grid. I tried it with rasterize, but here I can only chose last or first.
> Die you have any suggestions how I might get the sediment type for each
> raster cell?
>
> Thank you for your help!
>
> <°)))>< >°)))>< >°)))>< >°)))><
>
> Miriam Püts
> Marine Lebende Ressourcen/ Marine Living Resources
> Thünen-Institut für Seefischerei/ Thünen Institute of Sea fisheries
> Palmaille 9
> 22767 Hamburg (Germany)
>
> Tel:  +49 40 38905-105
> Mail: [hidden email]
>
> _______________________________________________
> 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

Re: netcdf file mapping and dimensions issue

Wed, 05/24/2017 - 04:42
Thanks a lot!  I need to fix the basics first though.  For example, even
though I can run the command and warp files at will, there is still an
issue I have yet to figure out.

The code below produces a nicely mapped file that looks like there is a
problem with the srcnodata or possibly dstnodata values as all of the pixel
values are less than zero for a scene that clearly has nonzero data.  I am
pretty stumped as I have looked at the documentation for the ocean color
data to get the srcnodata value of -32767s.  Also tried -32767, which is
shown as the fill value in SeaDAS band attributes table.  The link below
shows the result.  The data look quite like the original file as viewed in
SeaDAS but the values are clearly wrong.

r412<-raster('/Users/dmwarner/A2011202183000.L2_LAC_OC.x.nc',
var='geophysical_data/chl_ocx')
f <- "
https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/A2011202183000.L2_LAC_OC.x.nc"
download.file(f, basename(f), mode = "wb")
system('gdalwarp -srcnodata "None" -dstnodata None -srcnodata -32767s HDF5:"
A2011202183000.L2_LAC_OC.x.nc"://geophysical_data/Rrs_412 output.tif')
test<-raster('/Users/dmwarner/output.tif')
plot(test)

https://cdn.rawgit.com/dmwarn/Tethys/34794e44/warprrs412.png

On Tue, May 23, 2017 at 6:04 PM, Michael Sumner <[hidden email]> wrote:

> See gdalUtils for a wrapper, or ncdump to easily get the names of
> variables and other metadata. I have workers L2 here that might help you
> batch this and I'm happy to help:
>
> https://github.com/mdsumner/roc/blob/master/R/readL2.R
>
> We want to get deeper into this soon so I'll be looking at it.
>
> There's possibly a way to do all the sds in one GDAL call too.
>
> Cheers, Mike
>
> On Wed, 24 May 2017, 02:34 Warner, David <[hidden email]> wrote:
>
>> Now if I can just get calls to gdal functions to work from within R....
>>
>>
>> On Tue, May 23, 2017 at 11:39 AM, Michael Sumner <[hidden email]>
>> wrote:
>>
>>>
>>> I tried a few pretty unpromising things and then pointed gdalwarp at it.
>>> Seems fine, gdalwarp knows how to find geolocation arrays and use them, and
>>> it will auto-choose an output grid, and you can provide a custom one (i.e.
>>> with a local projection.
>>>
>>> The SDS (subdataset) naming thing is a bit awkward until you get used to
>>> it.  This was just GDAL 2.1.2 on Ubuntu.
>>>
>>>
>>> f <- "https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/
>>> A2011202183000.L2_LAC_OC.x.nc"
>>> download.file(f, basename(f), mode = "wb")
>>>
>>> system('gdalwarp HDF5:"A2011202183000.L2_LAC_OC.x.nc
>>> "://geophysical_data/chl_ocx chl_ocx .tif')
>>> Creating output file that is 783P x 530L.
>>> Processing input file HDF5:A2011202183000.L2_LAC_OC.x.nc
>>> ://geophysical_data/chl_ocx.
>>> 0...10...20...30...40...50...60...70...80...90...100 - done.
>>> r <- raster("chl_ocx ")
>>> r
>>> class       : RasterLayer
>>> dimensions  : 530, 783, 414990  (nrow, ncol, ncell)
>>> resolution  : 0.01224719, 0.01224719  (x, y)
>>> extent      : -86.99178, -77.40223, 40.73677, 47.22778  (xmin, xmax,
>>> ymin, ymax)
>>> coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
>>> +towgs84=0,0,0
>>>
>>> #plot(r)
>>> #library(mapdata)
>>> #map("worldHires", add = TRUE)
>>>
>>> Cheers, Mike.
>>>
>>>
>>> On Tue, 23 May 2017 at 23:58 Warner, David <[hidden email]> wrote:
>>>
>>>> Greetings all
>>>>
>>>> I am working nc files from the Ocean Biology Procesing Group.  The files
>>>> are MODIS aqua L2 ocean color files.
>>>>
>>>> These files are not mapped or perhaps some would say orthorectified such
>>>> that when you create a rasterLayer from one of the contained variables
>>>> and
>>>> plot it, the location and orientation of features is incorrect.  For
>>>> example, Georgian Bay, a large eastern subset of Lake Huron in the US
>>>> and
>>>> Canada, may well show up west of the lake.
>>>>
>>>> The dimensions of the data are read as column and row numbers rather
>>>> than
>>>> latitude and longitude.
>>>>
>>>> I am looking for a way to map these data in R other than my current
>>>> method.  The current method for a single nc file is shown below.  I run
>>>> this modified to fit into a parallelized foreach() loop that is used to
>>>> process folders of nc files.
>>>>
>>>> The  basic steps that map the nc file data are 1) create a data.frame
>>>> from
>>>> the nc files measured variable(s) and the latitude and longitude then
>>>> rasterize using vect2rast.SpatialPoints() and 2) resample the
>>>> rasterLayer
>>>> created from the data.frame using resample() to an nc file I reprojected
>>>> using SeaDAS.  This seems overly complicated considering the fact that I
>>>> start with a raster and it is pretty slow (30-40 seconds per nc file in
>>>> parallel, nearly two minutes per file in the script example here).
>>>>
>>>> I have two questions are 1) is there a way to eliminate the used of
>>>> vect2rast.SpatialPoints(), which is the slowest part of the process and
>>>> 2)
>>>> are there any other ways to perhaps speed the process up if I can't
>>>> eliminate vect2rast.SpatialPoints().  I have tried rasterize() and that
>>>> does not provide the correct product.  It recreates what I started with,
>>>> unmapped data.  What am I missing, or is this just a function of the
>>>> data I
>>>> am using?
>>>>
>>>>
>>>> Script is below.  I have included links to the two data files required
>>>> in
>>>> the process.  One is the nc file,
>>>> https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/
>>>> A2011202183000.L2_LAC_OC.x.nc
>>>> the other is the reprojected version of this file from SeaDAS.
>>>> https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/
>>>> A2011202183000.reprojected.tif
>>>> To use them you will have to change the path to the file locations.
>>>>
>>>> Thanks for your patience with my writing and poorly written code.  Any
>>>> ideas would be helpful.
>>>>
>>>> Dave
>>>>
>>>> library(raster)
>>>> library(plotKML)
>>>>
>>>> rasterOptions(maxmemory = 1e+09)
>>>> start <- Sys.time()
>>>> #Create rasterLayers from the bands of nc file that are required for
>>>> mapping, with r412 being the measured variable of interest
>>>> r412<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>>> var='geophysical_data/Rrs_412')
>>>> r443<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>>> var='geophysical_data/Rrs_443')
>>>> r469<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>>> var='geophysical_data/Rrs_469')
>>>> r488<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>>> var='geophysical_data/Rrs_488')
>>>> r531<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>>> var='geophysical_data/Rrs_531')
>>>> r547<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>>> var='geophysical_data/Rrs_547')
>>>> r555<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>>> var='geophysical_data/Rrs_555')
>>>> r645<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>>> var='geophysical_data/Rrs_645')
>>>> r667<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>>> var='geophysical_data/Rrs_667')
>>>> r678<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>>> var='geophysical_data/Rrs_678')
>>>>
>>>> longitude <- raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/
>>>> Huron/2011/
>>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>>> var="navigation_data/longitude")
>>>> latitude <- raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/
>>>> Huron/2011/
>>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>>> var="navigation_data/latitude")
>>>> r412  # shows that dimensions are row and column numbers, not lat and
>>>> long
>>>>
>>>>
>>>> #Import single reprojected band of the nc file as geotif, with
>>>> reprojection
>>>> done manually in SeaDAS software
>>>> correct<-raster("/Users/dmwarner/Documents/MODIS/OC/
>>>> correct/Huron/A2011202183000.reprojected.tif")
>>>>
>>>>
>>>> #now rasterize and map the nc file after
>>>>
>>>> d <- data.frame(x = values(longitude), y = values(latitude), r412 =
>>>> values(r412))
>>>> d$r412[is.na(d$r412)]<- 99999
>>>> coordinates(d)<-~x+y
>>>> dd<-vect2rast.SpatialPoints(d, fname ='r412', cell.size=0.009337697,
>>>> method='raster')
>>>> ras.412 <- raster(dd)
>>>> ras.412.crop<-crop(ras.412, myext)
>>>> #ras.412.crop[ras.412.crop<0]<-NA
>>>> ras.412.crop[ras.412.crop>90000]<- NA
>>>>
>>>>
>>>> d <- data.frame(x = values(longitude), y = values(latitude), r443 =
>>>> values(r443))
>>>> d$r443[is.na(d$r443)]<- 99999
>>>> coordinates(d)<-~x+y
>>>> dd<-vect2rast.SpatialPoints(d, fname ='r443', cell.size=0.009337697,
>>>> method='raster')
>>>> ras.443 <- raster(dd)
>>>> ras.443.crop<-crop(ras.443, myext)
>>>> ras.443.crop[ras.443.crop<0]<-NA
>>>> ras.443.crop[ras.443.crop>90000]<- NA
>>>>
>>>> d <- data.frame(x = values(longitude), y = values(latitude), r469 =
>>>> values(r469))
>>>> d$r469[is.na(d$r469)]<- 99999
>>>> coordinates(d)<-~x+y
>>>> dd<-vect2rast.SpatialPoints(d, fname ='r469', cell.size=0.009337697,
>>>> method='raster')
>>>> ras.469 <- raster(dd)
>>>> ras.469.crop<-crop(ras.469, myext)
>>>> ras.469.crop[ras.469.crop<0]<-NA
>>>> ras.469.crop[ras.469.crop>90000]<- NA
>>>>
>>>> d <- data.frame(x = values(longitude), y = values(latitude), r488 =
>>>> values(r488))
>>>> d$r488[is.na(d$r488)]<- 99999
>>>> coordinates(d)<-~x+y
>>>> dd<-vect2rast.SpatialPoints(d, fname ='r488', cell.size=0.009337697,
>>>> method='raster')
>>>> ras.488 <- raster(dd)
>>>> ras.488.crop<-crop(ras.488, myext)
>>>> ras.488.crop[ras.488.crop<0]<-NA
>>>> ras.488.crop[ras.488.crop>90000]<- NA
>>>>
>>>> d <- data.frame(x = values(longitude), y = values(latitude), r531 =
>>>> values(r531))
>>>> d$r531[is.na(d$r531)]<- 99999
>>>> coordinates(d)<-~x+y
>>>> dd<-vect2rast.SpatialPoints(d, fname ='r531', cell.size=0.009337697,
>>>> method='raster')
>>>> ras.531 <- raster(dd)
>>>> ras.531.crop<-crop(ras.531, myext)
>>>> ras.531.crop[ras.531.crop<0]<-NA
>>>> ras.531.crop[ras.531.crop>90000]<- NA
>>>>
>>>> d <- data.frame(x = values(longitude), y = values(latitude), r547 =
>>>> values(r547))
>>>> d$r547[is.na(d$r547)]<- 99999
>>>> coordinates(d)<-~x+y
>>>> dd<-vect2rast.SpatialPoints(d, fname ='r547', cell.size=0.009337697,
>>>> method='raster')
>>>> ras.547 <- raster(dd)
>>>> ras.547.crop<-crop(ras.547, myext)
>>>> ras.547.crop[ras.547.crop<0]<-NA
>>>> ras.547.crop[ras.547.crop>90000]<- NA
>>>>
>>>> d <- data.frame(x = values(longitude), y = values(latitude), r555 =
>>>> values(r555))
>>>> d$r555[is.na(d$r555)]<- 99999
>>>> coordinates(d)<-~x+y
>>>> dd<-vect2rast.SpatialPoints(d, fname ='r555', cell.size=0.009337697,
>>>> method='raster')
>>>> ras.555 <- raster(dd)
>>>> ras.555.crop<-crop(ras.555, myext)
>>>> ras.555.crop[ras.555.crop<0]<-NA
>>>> ras.555.crop[ras.555.crop>90000]<- NA
>>>>
>>>> d <- data.frame(x = values(longitude), y = values(latitude), r645 =
>>>> values(r645))
>>>> d$r645[is.na(d$r645)]<- 99999
>>>> coordinates(d)<-~x+y
>>>> dd<-vect2rast.SpatialPoints(d, fname ='r645', cell.size=0.009337697,
>>>> method='raster')
>>>> ras.645 <- raster(dd)
>>>> ras.645.crop<-crop(ras.645, myext)
>>>> ras.645.crop[ras.645.crop<0]<-NA
>>>> ras.645.crop[ras.645.crop>90000]<- NA
>>>>
>>>> d <- data.frame(x = values(longitude), y = values(latitude), r667 =
>>>> values(r667))
>>>> d$r667[is.na(d$r667)]<- 99999
>>>> coordinates(d)<-~x+y
>>>> dd<-vect2rast.SpatialPoints(d, fname ='r667', cell.size=0.009337697,
>>>> method='raster')
>>>> ras.667 <- raster(dd)
>>>> ras.667.crop<-crop(ras.667, myext)
>>>> ras.667.crop[ras.667.crop<0]<-NA
>>>> ras.667.crop[ras.667.crop>90000]<- NA
>>>>
>>>> d <- data.frame(x = values(longitude), y = values(latitude), r678 =
>>>> values(r678))
>>>> d$r678[is.na(d$r678)]<- 99999
>>>> coordinates(d)<-~x+y
>>>> dd<-vect2rast.SpatialPoints(d, fname ='r678', cell.size=0.009337697,
>>>> method='raster')
>>>> ras.678 <- raster(dd)
>>>> ras.678.crop<-crop(ras.678, myext)
>>>> ras.678.crop[ras.678.crop<0]<-NA
>>>> ras.678.crop[ras.678.crop>90000]<- NA
>>>>
>>>> myext <- extent(-84.9, -79.6, 42.9, 46.6) # Approx. extent that will
>>>> cover
>>>> all of Lake Huron
>>>> correct.crop<-crop(correct, myext)
>>>> ras.412.r<-resample(ras.412.crop, correct.crop, method='bilinear')
>>>> ras.443.r<-resample(ras.443.crop, correct.crop, method='bilinear')
>>>> ras.469.r<-resample(ras.469.crop, correct.crop, method='bilinear')
>>>> ras.488.r<-resample(ras.488.crop, correct.crop, method='bilinear')
>>>> ras.531.r<-resample(ras.531.crop, correct.crop, method='bilinear')
>>>> ras.547.r<-resample(ras.547.crop, correct.crop, method='bilinear')
>>>> ras.555.r<-resample(ras.555.crop, correct.crop, method='bilinear')
>>>> ras.645.r<-resample(ras.645.crop, correct.crop, method='bilinear')
>>>> ras.667.r<-resample(ras.667.crop, correct.crop, method='bilinear')
>>>> ras.678.r<-resample(ras.678.crop, correct.crop, method='bilinear')
>>>> End<-Sys.time()
>>>> difftime(End, start)
>>>>
>>>> #plot original unmapped nc band Rrs_412 versus mapped version for gross
>>>> comparison
>>>> par(mfrow=c(2,1), mar=c(2,2,2,2))
>>>> my.orig.extent<-extent(100,475,200,600)
>>>> r412.crop<-crop(r412, my.orig.extent)
>>>> plot(r412.crop, main='Original unmapped')
>>>> plot(ras.412.r, main='Mapped in R')
>>>>
>>>> --
>>>> David Warner
>>>> Research Fisheries Biologist
>>>> U.S.G.S. Great Lakes Science Center
>>>> 1451 Green Road
>>>> Ann Arbor, MI 48105
>>>> 734-214-9392 <(734)%20214-9392>
>>>>
>>>>         [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> 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
>>>
>>>
>>
>>
>> --
>> David Warner
>> Research Fisheries Biologist
>> U.S.G.S. Great Lakes Science Center
>> 1451 Green Road
>> Ann Arbor, MI 48105
>> 734-214-9392
>>
> --
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> Kingston Tasmania 7050 Australia
>
>

--
David Warner
Research Fisheries Biologist
U.S.G.S. Great Lakes Science Center
1451 Green Road
Ann Arbor, MI 48105
734-214-9392

        [[alternative HTML version deleted]]

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

Re: Problem in ENFA: infinite or missing values

Wed, 05/24/2017 - 04:29
Can you you put together a reproducible example ( reprex)? I'm happy to
help but the set up is too onerous.

I'd recommend using raster to count points and for gridded data generally
but the advice is perhaps too abstract without an example. cellFromXY gives
the right output for a point in cell count as a grouped aggregate summary,
 but it's a module requiring to be embedded in a workflow (which means it's
super powerful, efficient and flexible) rather than a high level solution
requiring no further composition.

Cheers, Mike


On Fri, 19 May 2017, 16:51 Perry Beasley-hall <
[hidden email]> wrote:

> Hello everyone,
> I'm trying to do an ENFA using the package adehabitatHS. I'm following the
> example code and am creating the 'pr' file using:
>
> pr <- slot(count.points(species, map), "data")[,1]
>
> Where 'species' is a .csv file with two columns of latitude and longitude,
> and 'map' is a SpatialGridDataFrame file of 23 .asc inputs. For reference,
> I'm importing my .asc files using readGDAL and then using the cbind command
> to put them into the 'map' file.
>
> I encounter the following warning message when trying to create 'pr':
>
> 1: In points2grid(points, tolerance, round) :
>   grid has empty column/rows in dimension 1
> 2: In points2grid(points, tolerance, round) :
>   grid has empty column/rows in dimension 2
>
> I assumed this would be alright, but I'm getting the subsequent error when
> I attempt the enfa command:
>
> enfa1 <- enfa(pc, pr, scannf = FALSE)
>
> Error in eigen(Se) : infinite or missing values in 'x'
>
> Some searching online has told me that this might mean not all of my rows
> in the 'species' .csv correspond to variables in the 'map' file, but I
> don't think this is the case for my dataset.
>
> Can anyone help me with this issue? Many thanks in advance.
> ____
>
> Perry Beasley-Hall
> PhD Candidate
> School of Life and Environmental Sciences
>
> THE UNIVERSITY OF SYDNEY
> Room 329, Edgeworth David Building A11
> Sydney NSW 2006, Australia
> E-mail: [hidden email]
> Phone: +61 420 494 584
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> 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

Deal with multiple factorlevel in one grid square

Wed, 05/24/2017 - 02:12
Hi everyone,

I have the following problem: I have a personalized grid and a shape file with polygons representing sediment types. Now I would like to apply this grid to the Polygons to identify the sediment type most common within each grid. I tried it with rasterize, but here I can only chose last or first. Die you have any suggestions how I might get the sediment type for each raster cell?

Thank you for your help!

<°)))>< >°)))>< >°)))>< >°)))><

Miriam Püts
Marine Lebende Ressourcen/ Marine Living Resources
Thünen-Institut für Seefischerei/ Thünen Institute of Sea fisheries
Palmaille 9
22767 Hamburg (Germany)

Tel:  +49 40 38905-105
Mail: [hidden email]

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

Re: netcdf file mapping and dimensions issue

Tue, 05/23/2017 - 17:04
See gdalUtils for a wrapper, or ncdump to easily get the names of variables
and other metadata. I have workers L2 here that might help you batch this
and I'm happy to help:

https://github.com/mdsumner/roc/blob/master/R/readL2.R

We want to get deeper into this soon so I'll be looking at it.

There's possibly a way to do all the sds in one GDAL call too.

Cheers, Mike

On Wed, 24 May 2017, 02:34 Warner, David <[hidden email]> wrote:

> Now if I can just get calls to gdal functions to work from within R....
>
>
> On Tue, May 23, 2017 at 11:39 AM, Michael Sumner <[hidden email]>
> wrote:
>
>>
>> I tried a few pretty unpromising things and then pointed gdalwarp at it.
>> Seems fine, gdalwarp knows how to find geolocation arrays and use them, and
>> it will auto-choose an output grid, and you can provide a custom one (i.e.
>> with a local projection.
>>
>> The SDS (subdataset) naming thing is a bit awkward until you get used to
>> it.  This was just GDAL 2.1.2 on Ubuntu.
>>
>>
>> f <- "
>> https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/A2011202183000.L2_LAC_OC.x.nc
>> "
>> download.file(f, basename(f), mode = "wb")
>>
>> system('gdalwarp HDF5:"A2011202183000.L2_LAC_OC.x.nc"://geophysical_data/chl_ocx
>> chl_ocx .tif')
>> Creating output file that is 783P x 530L.
>> Processing input file HDF5:A2011202183000.L2_LAC_OC.x.nc
>> ://geophysical_data/chl_ocx.
>> 0...10...20...30...40...50...60...70...80...90...100 - done.
>> r <- raster("chl_ocx ")
>> r
>> class       : RasterLayer
>> dimensions  : 530, 783, 414990  (nrow, ncol, ncell)
>> resolution  : 0.01224719, 0.01224719  (x, y)
>> extent      : -86.99178, -77.40223, 40.73677, 47.22778  (xmin, xmax,
>> ymin, ymax)
>> coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
>> +towgs84=0,0,0
>>
>> #plot(r)
>> #library(mapdata)
>> #map("worldHires", add = TRUE)
>>
>> Cheers, Mike.
>>
>>
>> On Tue, 23 May 2017 at 23:58 Warner, David <[hidden email]> wrote:
>>
>>> Greetings all
>>>
>>> I am working nc files from the Ocean Biology Procesing Group.  The files
>>> are MODIS aqua L2 ocean color files.
>>>
>>> These files are not mapped or perhaps some would say orthorectified such
>>> that when you create a rasterLayer from one of the contained variables
>>> and
>>> plot it, the location and orientation of features is incorrect.  For
>>> example, Georgian Bay, a large eastern subset of Lake Huron in the US and
>>> Canada, may well show up west of the lake.
>>>
>>> The dimensions of the data are read as column and row numbers rather than
>>> latitude and longitude.
>>>
>>> I am looking for a way to map these data in R other than my current
>>> method.  The current method for a single nc file is shown below.  I run
>>> this modified to fit into a parallelized foreach() loop that is used to
>>> process folders of nc files.
>>>
>>> The  basic steps that map the nc file data are 1) create a data.frame
>>> from
>>> the nc files measured variable(s) and the latitude and longitude then
>>> rasterize using vect2rast.SpatialPoints() and 2) resample the rasterLayer
>>> created from the data.frame using resample() to an nc file I reprojected
>>> using SeaDAS.  This seems overly complicated considering the fact that I
>>> start with a raster and it is pretty slow (30-40 seconds per nc file in
>>> parallel, nearly two minutes per file in the script example here).
>>>
>>> I have two questions are 1) is there a way to eliminate the used of
>>> vect2rast.SpatialPoints(), which is the slowest part of the process and
>>> 2)
>>> are there any other ways to perhaps speed the process up if I can't
>>> eliminate vect2rast.SpatialPoints().  I have tried rasterize() and that
>>> does not provide the correct product.  It recreates what I started with,
>>> unmapped data.  What am I missing, or is this just a function of the
>>> data I
>>> am using?
>>>
>>>
>>> Script is below.  I have included links to the two data files required in
>>> the process.  One is the nc file,
>>>
>>> https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/A2011202183000.L2_LAC_OC.x.nc
>>> the other is the reprojected version of this file from SeaDAS.
>>>
>>> https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/A2011202183000.reprojected.tif
>>> To use them you will have to change the path to the file locations.
>>>
>>> Thanks for your patience with my writing and poorly written code.  Any
>>> ideas would be helpful.
>>>
>>> Dave
>>>
>>> library(raster)
>>> library(plotKML)
>>>
>>> rasterOptions(maxmemory = 1e+09)
>>> start <- Sys.time()
>>> #Create rasterLayers from the bands of nc file that are required for
>>> mapping, with r412 being the measured variable of interest
>>> r412<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>> var='geophysical_data/Rrs_412')
>>> r443<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>> var='geophysical_data/Rrs_443')
>>> r469<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>> var='geophysical_data/Rrs_469')
>>> r488<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>> var='geophysical_data/Rrs_488')
>>> r531<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>> var='geophysical_data/Rrs_531')
>>> r547<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>> var='geophysical_data/Rrs_547')
>>> r555<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>> var='geophysical_data/Rrs_555')
>>> r645<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>> var='geophysical_data/Rrs_645')
>>> r667<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>> var='geophysical_data/Rrs_667')
>>> r678<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>> var='geophysical_data/Rrs_678')
>>>
>>> longitude <-
>>> raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>> var="navigation_data/longitude")
>>> latitude <-
>>> raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>>> var="navigation_data/latitude")
>>> r412  # shows that dimensions are row and column numbers, not lat and
>>> long
>>>
>>>
>>> #Import single reprojected band of the nc file as geotif, with
>>> reprojection
>>> done manually in SeaDAS software
>>> correct<-raster("/Users/dmwarner/Documents/MODIS/OC/
>>> correct/Huron/A2011202183000.reprojected.tif")
>>>
>>>
>>> #now rasterize and map the nc file after
>>>
>>> d <- data.frame(x = values(longitude), y = values(latitude), r412 =
>>> values(r412))
>>> d$r412[is.na(d$r412)]<- 99999
>>> coordinates(d)<-~x+y
>>> dd<-vect2rast.SpatialPoints(d, fname ='r412', cell.size=0.009337697,
>>> method='raster')
>>> ras.412 <- raster(dd)
>>> ras.412.crop<-crop(ras.412, myext)
>>> #ras.412.crop[ras.412.crop<0]<-NA
>>> ras.412.crop[ras.412.crop>90000]<- NA
>>>
>>>
>>> d <- data.frame(x = values(longitude), y = values(latitude), r443 =
>>> values(r443))
>>> d$r443[is.na(d$r443)]<- 99999
>>> coordinates(d)<-~x+y
>>> dd<-vect2rast.SpatialPoints(d, fname ='r443', cell.size=0.009337697,
>>> method='raster')
>>> ras.443 <- raster(dd)
>>> ras.443.crop<-crop(ras.443, myext)
>>> ras.443.crop[ras.443.crop<0]<-NA
>>> ras.443.crop[ras.443.crop>90000]<- NA
>>>
>>> d <- data.frame(x = values(longitude), y = values(latitude), r469 =
>>> values(r469))
>>> d$r469[is.na(d$r469)]<- 99999
>>> coordinates(d)<-~x+y
>>> dd<-vect2rast.SpatialPoints(d, fname ='r469', cell.size=0.009337697,
>>> method='raster')
>>> ras.469 <- raster(dd)
>>> ras.469.crop<-crop(ras.469, myext)
>>> ras.469.crop[ras.469.crop<0]<-NA
>>> ras.469.crop[ras.469.crop>90000]<- NA
>>>
>>> d <- data.frame(x = values(longitude), y = values(latitude), r488 =
>>> values(r488))
>>> d$r488[is.na(d$r488)]<- 99999
>>> coordinates(d)<-~x+y
>>> dd<-vect2rast.SpatialPoints(d, fname ='r488', cell.size=0.009337697,
>>> method='raster')
>>> ras.488 <- raster(dd)
>>> ras.488.crop<-crop(ras.488, myext)
>>> ras.488.crop[ras.488.crop<0]<-NA
>>> ras.488.crop[ras.488.crop>90000]<- NA
>>>
>>> d <- data.frame(x = values(longitude), y = values(latitude), r531 =
>>> values(r531))
>>> d$r531[is.na(d$r531)]<- 99999
>>> coordinates(d)<-~x+y
>>> dd<-vect2rast.SpatialPoints(d, fname ='r531', cell.size=0.009337697,
>>> method='raster')
>>> ras.531 <- raster(dd)
>>> ras.531.crop<-crop(ras.531, myext)
>>> ras.531.crop[ras.531.crop<0]<-NA
>>> ras.531.crop[ras.531.crop>90000]<- NA
>>>
>>> d <- data.frame(x = values(longitude), y = values(latitude), r547 =
>>> values(r547))
>>> d$r547[is.na(d$r547)]<- 99999
>>> coordinates(d)<-~x+y
>>> dd<-vect2rast.SpatialPoints(d, fname ='r547', cell.size=0.009337697,
>>> method='raster')
>>> ras.547 <- raster(dd)
>>> ras.547.crop<-crop(ras.547, myext)
>>> ras.547.crop[ras.547.crop<0]<-NA
>>> ras.547.crop[ras.547.crop>90000]<- NA
>>>
>>> d <- data.frame(x = values(longitude), y = values(latitude), r555 =
>>> values(r555))
>>> d$r555[is.na(d$r555)]<- 99999
>>> coordinates(d)<-~x+y
>>> dd<-vect2rast.SpatialPoints(d, fname ='r555', cell.size=0.009337697,
>>> method='raster')
>>> ras.555 <- raster(dd)
>>> ras.555.crop<-crop(ras.555, myext)
>>> ras.555.crop[ras.555.crop<0]<-NA
>>> ras.555.crop[ras.555.crop>90000]<- NA
>>>
>>> d <- data.frame(x = values(longitude), y = values(latitude), r645 =
>>> values(r645))
>>> d$r645[is.na(d$r645)]<- 99999
>>> coordinates(d)<-~x+y
>>> dd<-vect2rast.SpatialPoints(d, fname ='r645', cell.size=0.009337697,
>>> method='raster')
>>> ras.645 <- raster(dd)
>>> ras.645.crop<-crop(ras.645, myext)
>>> ras.645.crop[ras.645.crop<0]<-NA
>>> ras.645.crop[ras.645.crop>90000]<- NA
>>>
>>> d <- data.frame(x = values(longitude), y = values(latitude), r667 =
>>> values(r667))
>>> d$r667[is.na(d$r667)]<- 99999
>>> coordinates(d)<-~x+y
>>> dd<-vect2rast.SpatialPoints(d, fname ='r667', cell.size=0.009337697,
>>> method='raster')
>>> ras.667 <- raster(dd)
>>> ras.667.crop<-crop(ras.667, myext)
>>> ras.667.crop[ras.667.crop<0]<-NA
>>> ras.667.crop[ras.667.crop>90000]<- NA
>>>
>>> d <- data.frame(x = values(longitude), y = values(latitude), r678 =
>>> values(r678))
>>> d$r678[is.na(d$r678)]<- 99999
>>> coordinates(d)<-~x+y
>>> dd<-vect2rast.SpatialPoints(d, fname ='r678', cell.size=0.009337697,
>>> method='raster')
>>> ras.678 <- raster(dd)
>>> ras.678.crop<-crop(ras.678, myext)
>>> ras.678.crop[ras.678.crop<0]<-NA
>>> ras.678.crop[ras.678.crop>90000]<- NA
>>>
>>> myext <- extent(-84.9, -79.6, 42.9, 46.6) # Approx. extent that will
>>> cover
>>> all of Lake Huron
>>> correct.crop<-crop(correct, myext)
>>> ras.412.r<-resample(ras.412.crop, correct.crop, method='bilinear')
>>> ras.443.r<-resample(ras.443.crop, correct.crop, method='bilinear')
>>> ras.469.r<-resample(ras.469.crop, correct.crop, method='bilinear')
>>> ras.488.r<-resample(ras.488.crop, correct.crop, method='bilinear')
>>> ras.531.r<-resample(ras.531.crop, correct.crop, method='bilinear')
>>> ras.547.r<-resample(ras.547.crop, correct.crop, method='bilinear')
>>> ras.555.r<-resample(ras.555.crop, correct.crop, method='bilinear')
>>> ras.645.r<-resample(ras.645.crop, correct.crop, method='bilinear')
>>> ras.667.r<-resample(ras.667.crop, correct.crop, method='bilinear')
>>> ras.678.r<-resample(ras.678.crop, correct.crop, method='bilinear')
>>> End<-Sys.time()
>>> difftime(End, start)
>>>
>>> #plot original unmapped nc band Rrs_412 versus mapped version for gross
>>> comparison
>>> par(mfrow=c(2,1), mar=c(2,2,2,2))
>>> my.orig.extent<-extent(100,475,200,600)
>>> r412.crop<-crop(r412, my.orig.extent)
>>> plot(r412.crop, main='Original unmapped')
>>> plot(ras.412.r, main='Mapped in R')
>>>
>>> --
>>> David Warner
>>> Research Fisheries Biologist
>>> U.S.G.S. Great Lakes Science Center
>>> 1451 Green Road
>>> Ann Arbor, MI 48105
>>> 734-214-9392 <(734)%20214-9392>
>>>
>>>         [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> 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
>>
>>
>
>
> --
> David Warner
> Research Fisheries Biologist
> U.S.G.S. Great Lakes Science Center
> 1451 Green Road
> Ann Arbor, MI 48105
> 734-214-9392
> --
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

Re: netcdf file mapping and dimensions issue

Tue, 05/23/2017 - 11:33
Now if I can just get calls to gdal functions to work from within R....


On Tue, May 23, 2017 at 11:39 AM, Michael Sumner <[hidden email]> wrote:

>
> I tried a few pretty unpromising things and then pointed gdalwarp at it.
> Seems fine, gdalwarp knows how to find geolocation arrays and use them, and
> it will auto-choose an output grid, and you can provide a custom one (i.e.
> with a local projection.
>
> The SDS (subdataset) naming thing is a bit awkward until you get used to
> it.  This was just GDAL 2.1.2 on Ubuntu.
>
>
> f <- "https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/
> A2011202183000.L2_LAC_OC.x.nc"
> download.file(f, basename(f), mode = "wb")
>
> system('gdalwarp HDF5:"A2011202183000.L2_LAC_OC.x.nc"://geophysical_data/chl_ocx
> chl_ocx .tif')
> Creating output file that is 783P x 530L.
> Processing input file HDF5:A2011202183000.L2_LAC_OC.x.nc
> ://geophysical_data/chl_ocx.
> 0...10...20...30...40...50...60...70...80...90...100 - done.
> r <- raster("chl_ocx ")
> r
> class       : RasterLayer
> dimensions  : 530, 783, 414990  (nrow, ncol, ncell)
> resolution  : 0.01224719, 0.01224719  (x, y)
> extent      : -86.99178, -77.40223, 40.73677, 47.22778  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
> +towgs84=0,0,0
>
> #plot(r)
> #library(mapdata)
> #map("worldHires", add = TRUE)
>
> Cheers, Mike.
>
>
> On Tue, 23 May 2017 at 23:58 Warner, David <[hidden email]> wrote:
>
>> Greetings all
>>
>> I am working nc files from the Ocean Biology Procesing Group.  The files
>> are MODIS aqua L2 ocean color files.
>>
>> These files are not mapped or perhaps some would say orthorectified such
>> that when you create a rasterLayer from one of the contained variables and
>> plot it, the location and orientation of features is incorrect.  For
>> example, Georgian Bay, a large eastern subset of Lake Huron in the US and
>> Canada, may well show up west of the lake.
>>
>> The dimensions of the data are read as column and row numbers rather than
>> latitude and longitude.
>>
>> I am looking for a way to map these data in R other than my current
>> method.  The current method for a single nc file is shown below.  I run
>> this modified to fit into a parallelized foreach() loop that is used to
>> process folders of nc files.
>>
>> The  basic steps that map the nc file data are 1) create a data.frame from
>> the nc files measured variable(s) and the latitude and longitude then
>> rasterize using vect2rast.SpatialPoints() and 2) resample the rasterLayer
>> created from the data.frame using resample() to an nc file I reprojected
>> using SeaDAS.  This seems overly complicated considering the fact that I
>> start with a raster and it is pretty slow (30-40 seconds per nc file in
>> parallel, nearly two minutes per file in the script example here).
>>
>> I have two questions are 1) is there a way to eliminate the used of
>> vect2rast.SpatialPoints(), which is the slowest part of the process and 2)
>> are there any other ways to perhaps speed the process up if I can't
>> eliminate vect2rast.SpatialPoints().  I have tried rasterize() and that
>> does not provide the correct product.  It recreates what I started with,
>> unmapped data.  What am I missing, or is this just a function of the data
>> I
>> am using?
>>
>>
>> Script is below.  I have included links to the two data files required in
>> the process.  One is the nc file,
>> https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/
>> A2011202183000.L2_LAC_OC.x.nc
>> the other is the reprojected version of this file from SeaDAS.
>> https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/
>> A2011202183000.reprojected.tif
>> To use them you will have to change the path to the file locations.
>>
>> Thanks for your patience with my writing and poorly written code.  Any
>> ideas would be helpful.
>>
>> Dave
>>
>> library(raster)
>> library(plotKML)
>>
>> rasterOptions(maxmemory = 1e+09)
>> start <- Sys.time()
>> #Create rasterLayers from the bands of nc file that are required for
>> mapping, with r412 being the measured variable of interest
>> r412<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_412')
>> r443<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_443')
>> r469<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_469')
>> r488<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_488')
>> r531<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_531')
>> r547<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_547')
>> r555<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_555')
>> r645<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_645')
>> r667<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_667')
>> r678<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_678')
>>
>> longitude <- raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/
>> Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var="navigation_data/longitude")
>> latitude <- raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/
>> Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var="navigation_data/latitude")
>> r412  # shows that dimensions are row and column numbers, not lat and long
>>
>>
>> #Import single reprojected band of the nc file as geotif, with
>> reprojection
>> done manually in SeaDAS software
>> correct<-raster("/Users/dmwarner/Documents/MODIS/OC/
>> correct/Huron/A2011202183000.reprojected.tif")
>>
>>
>> #now rasterize and map the nc file after
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r412 =
>> values(r412))
>> d$r412[is.na(d$r412)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r412', cell.size=0.009337697,
>> method='raster')
>> ras.412 <- raster(dd)
>> ras.412.crop<-crop(ras.412, myext)
>> #ras.412.crop[ras.412.crop<0]<-NA
>> ras.412.crop[ras.412.crop>90000]<- NA
>>
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r443 =
>> values(r443))
>> d$r443[is.na(d$r443)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r443', cell.size=0.009337697,
>> method='raster')
>> ras.443 <- raster(dd)
>> ras.443.crop<-crop(ras.443, myext)
>> ras.443.crop[ras.443.crop<0]<-NA
>> ras.443.crop[ras.443.crop>90000]<- NA
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r469 =
>> values(r469))
>> d$r469[is.na(d$r469)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r469', cell.size=0.009337697,
>> method='raster')
>> ras.469 <- raster(dd)
>> ras.469.crop<-crop(ras.469, myext)
>> ras.469.crop[ras.469.crop<0]<-NA
>> ras.469.crop[ras.469.crop>90000]<- NA
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r488 =
>> values(r488))
>> d$r488[is.na(d$r488)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r488', cell.size=0.009337697,
>> method='raster')
>> ras.488 <- raster(dd)
>> ras.488.crop<-crop(ras.488, myext)
>> ras.488.crop[ras.488.crop<0]<-NA
>> ras.488.crop[ras.488.crop>90000]<- NA
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r531 =
>> values(r531))
>> d$r531[is.na(d$r531)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r531', cell.size=0.009337697,
>> method='raster')
>> ras.531 <- raster(dd)
>> ras.531.crop<-crop(ras.531, myext)
>> ras.531.crop[ras.531.crop<0]<-NA
>> ras.531.crop[ras.531.crop>90000]<- NA
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r547 =
>> values(r547))
>> d$r547[is.na(d$r547)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r547', cell.size=0.009337697,
>> method='raster')
>> ras.547 <- raster(dd)
>> ras.547.crop<-crop(ras.547, myext)
>> ras.547.crop[ras.547.crop<0]<-NA
>> ras.547.crop[ras.547.crop>90000]<- NA
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r555 =
>> values(r555))
>> d$r555[is.na(d$r555)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r555', cell.size=0.009337697,
>> method='raster')
>> ras.555 <- raster(dd)
>> ras.555.crop<-crop(ras.555, myext)
>> ras.555.crop[ras.555.crop<0]<-NA
>> ras.555.crop[ras.555.crop>90000]<- NA
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r645 =
>> values(r645))
>> d$r645[is.na(d$r645)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r645', cell.size=0.009337697,
>> method='raster')
>> ras.645 <- raster(dd)
>> ras.645.crop<-crop(ras.645, myext)
>> ras.645.crop[ras.645.crop<0]<-NA
>> ras.645.crop[ras.645.crop>90000]<- NA
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r667 =
>> values(r667))
>> d$r667[is.na(d$r667)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r667', cell.size=0.009337697,
>> method='raster')
>> ras.667 <- raster(dd)
>> ras.667.crop<-crop(ras.667, myext)
>> ras.667.crop[ras.667.crop<0]<-NA
>> ras.667.crop[ras.667.crop>90000]<- NA
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r678 =
>> values(r678))
>> d$r678[is.na(d$r678)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r678', cell.size=0.009337697,
>> method='raster')
>> ras.678 <- raster(dd)
>> ras.678.crop<-crop(ras.678, myext)
>> ras.678.crop[ras.678.crop<0]<-NA
>> ras.678.crop[ras.678.crop>90000]<- NA
>>
>> myext <- extent(-84.9, -79.6, 42.9, 46.6) # Approx. extent that will cover
>> all of Lake Huron
>> correct.crop<-crop(correct, myext)
>> ras.412.r<-resample(ras.412.crop, correct.crop, method='bilinear')
>> ras.443.r<-resample(ras.443.crop, correct.crop, method='bilinear')
>> ras.469.r<-resample(ras.469.crop, correct.crop, method='bilinear')
>> ras.488.r<-resample(ras.488.crop, correct.crop, method='bilinear')
>> ras.531.r<-resample(ras.531.crop, correct.crop, method='bilinear')
>> ras.547.r<-resample(ras.547.crop, correct.crop, method='bilinear')
>> ras.555.r<-resample(ras.555.crop, correct.crop, method='bilinear')
>> ras.645.r<-resample(ras.645.crop, correct.crop, method='bilinear')
>> ras.667.r<-resample(ras.667.crop, correct.crop, method='bilinear')
>> ras.678.r<-resample(ras.678.crop, correct.crop, method='bilinear')
>> End<-Sys.time()
>> difftime(End, start)
>>
>> #plot original unmapped nc band Rrs_412 versus mapped version for gross
>> comparison
>> par(mfrow=c(2,1), mar=c(2,2,2,2))
>> my.orig.extent<-extent(100,475,200,600)
>> r412.crop<-crop(r412, my.orig.extent)
>> plot(r412.crop, main='Original unmapped')
>> plot(ras.412.r, main='Mapped in R')
>>
>> --
>> David Warner
>> Research Fisheries Biologist
>> U.S.G.S. Great Lakes Science Center
>> 1451 Green Road
>> Ann Arbor, MI 48105
>> 734-214-9392 <(734)%20214-9392>
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> 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
>
>

--
David Warner
Research Fisheries Biologist
U.S.G.S. Great Lakes Science Center
1451 Green Road
Ann Arbor, MI 48105
734-214-9392

        [[alternative HTML version deleted]]

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

Re: netcdf file mapping and dimensions issue

Tue, 05/23/2017 - 10:43
This looks promising!  Thank you!
Dave

On Tue, May 23, 2017 at 11:39 AM, Michael Sumner <[hidden email]> wrote:

>
> I tried a few pretty unpromising things and then pointed gdalwarp at it.
> Seems fine, gdalwarp knows how to find geolocation arrays and use them, and
> it will auto-choose an output grid, and you can provide a custom one (i.e.
> with a local projection.
>
> The SDS (subdataset) naming thing is a bit awkward until you get used to
> it.  This was just GDAL 2.1.2 on Ubuntu.
>
>
> f <- "https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/
> A2011202183000.L2_LAC_OC.x.nc"
> download.file(f, basename(f), mode = "wb")
>
> system('gdalwarp HDF5:"A2011202183000.L2_LAC_OC.x.nc"://geophysical_data/chl_ocx
> chl_ocx .tif')
> Creating output file that is 783P x 530L.
> Processing input file HDF5:A2011202183000.L2_LAC_OC.x.nc
> ://geophysical_data/chl_ocx.
> 0...10...20...30...40...50...60...70...80...90...100 - done.
> r <- raster("chl_ocx ")
> r
> class       : RasterLayer
> dimensions  : 530, 783, 414990  (nrow, ncol, ncell)
> resolution  : 0.01224719, 0.01224719  (x, y)
> extent      : -86.99178, -77.40223, 40.73677, 47.22778  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
> +towgs84=0,0,0
>
> #plot(r)
> #library(mapdata)
> #map("worldHires", add = TRUE)
>
> Cheers, Mike.
>
>
> On Tue, 23 May 2017 at 23:58 Warner, David <[hidden email]> wrote:
>
>> Greetings all
>>
>> I am working nc files from the Ocean Biology Procesing Group.  The files
>> are MODIS aqua L2 ocean color files.
>>
>> These files are not mapped or perhaps some would say orthorectified such
>> that when you create a rasterLayer from one of the contained variables and
>> plot it, the location and orientation of features is incorrect.  For
>> example, Georgian Bay, a large eastern subset of Lake Huron in the US and
>> Canada, may well show up west of the lake.
>>
>> The dimensions of the data are read as column and row numbers rather than
>> latitude and longitude.
>>
>> I am looking for a way to map these data in R other than my current
>> method.  The current method for a single nc file is shown below.  I run
>> this modified to fit into a parallelized foreach() loop that is used to
>> process folders of nc files.
>>
>> The  basic steps that map the nc file data are 1) create a data.frame from
>> the nc files measured variable(s) and the latitude and longitude then
>> rasterize using vect2rast.SpatialPoints() and 2) resample the rasterLayer
>> created from the data.frame using resample() to an nc file I reprojected
>> using SeaDAS.  This seems overly complicated considering the fact that I
>> start with a raster and it is pretty slow (30-40 seconds per nc file in
>> parallel, nearly two minutes per file in the script example here).
>>
>> I have two questions are 1) is there a way to eliminate the used of
>> vect2rast.SpatialPoints(), which is the slowest part of the process and 2)
>> are there any other ways to perhaps speed the process up if I can't
>> eliminate vect2rast.SpatialPoints().  I have tried rasterize() and that
>> does not provide the correct product.  It recreates what I started with,
>> unmapped data.  What am I missing, or is this just a function of the data
>> I
>> am using?
>>
>>
>> Script is below.  I have included links to the two data files required in
>> the process.  One is the nc file,
>> https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/
>> A2011202183000.L2_LAC_OC.x.nc
>> the other is the reprojected version of this file from SeaDAS.
>> https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/
>> A2011202183000.reprojected.tif
>> To use them you will have to change the path to the file locations.
>>
>> Thanks for your patience with my writing and poorly written code.  Any
>> ideas would be helpful.
>>
>> Dave
>>
>> library(raster)
>> library(plotKML)
>>
>> rasterOptions(maxmemory = 1e+09)
>> start <- Sys.time()
>> #Create rasterLayers from the bands of nc file that are required for
>> mapping, with r412 being the measured variable of interest
>> r412<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_412')
>> r443<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_443')
>> r469<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_469')
>> r488<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_488')
>> r531<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_531')
>> r547<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_547')
>> r555<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_555')
>> r645<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_645')
>> r667<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_667')
>> r678<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var='geophysical_data/Rrs_678')
>>
>> longitude <- raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/
>> Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var="navigation_data/longitude")
>> latitude <- raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/
>> Huron/2011/
>> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
>> var="navigation_data/latitude")
>> r412  # shows that dimensions are row and column numbers, not lat and long
>>
>>
>> #Import single reprojected band of the nc file as geotif, with
>> reprojection
>> done manually in SeaDAS software
>> correct<-raster("/Users/dmwarner/Documents/MODIS/OC/
>> correct/Huron/A2011202183000.reprojected.tif")
>>
>>
>> #now rasterize and map the nc file after
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r412 =
>> values(r412))
>> d$r412[is.na(d$r412)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r412', cell.size=0.009337697,
>> method='raster')
>> ras.412 <- raster(dd)
>> ras.412.crop<-crop(ras.412, myext)
>> #ras.412.crop[ras.412.crop<0]<-NA
>> ras.412.crop[ras.412.crop>90000]<- NA
>>
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r443 =
>> values(r443))
>> d$r443[is.na(d$r443)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r443', cell.size=0.009337697,
>> method='raster')
>> ras.443 <- raster(dd)
>> ras.443.crop<-crop(ras.443, myext)
>> ras.443.crop[ras.443.crop<0]<-NA
>> ras.443.crop[ras.443.crop>90000]<- NA
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r469 =
>> values(r469))
>> d$r469[is.na(d$r469)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r469', cell.size=0.009337697,
>> method='raster')
>> ras.469 <- raster(dd)
>> ras.469.crop<-crop(ras.469, myext)
>> ras.469.crop[ras.469.crop<0]<-NA
>> ras.469.crop[ras.469.crop>90000]<- NA
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r488 =
>> values(r488))
>> d$r488[is.na(d$r488)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r488', cell.size=0.009337697,
>> method='raster')
>> ras.488 <- raster(dd)
>> ras.488.crop<-crop(ras.488, myext)
>> ras.488.crop[ras.488.crop<0]<-NA
>> ras.488.crop[ras.488.crop>90000]<- NA
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r531 =
>> values(r531))
>> d$r531[is.na(d$r531)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r531', cell.size=0.009337697,
>> method='raster')
>> ras.531 <- raster(dd)
>> ras.531.crop<-crop(ras.531, myext)
>> ras.531.crop[ras.531.crop<0]<-NA
>> ras.531.crop[ras.531.crop>90000]<- NA
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r547 =
>> values(r547))
>> d$r547[is.na(d$r547)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r547', cell.size=0.009337697,
>> method='raster')
>> ras.547 <- raster(dd)
>> ras.547.crop<-crop(ras.547, myext)
>> ras.547.crop[ras.547.crop<0]<-NA
>> ras.547.crop[ras.547.crop>90000]<- NA
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r555 =
>> values(r555))
>> d$r555[is.na(d$r555)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r555', cell.size=0.009337697,
>> method='raster')
>> ras.555 <- raster(dd)
>> ras.555.crop<-crop(ras.555, myext)
>> ras.555.crop[ras.555.crop<0]<-NA
>> ras.555.crop[ras.555.crop>90000]<- NA
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r645 =
>> values(r645))
>> d$r645[is.na(d$r645)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r645', cell.size=0.009337697,
>> method='raster')
>> ras.645 <- raster(dd)
>> ras.645.crop<-crop(ras.645, myext)
>> ras.645.crop[ras.645.crop<0]<-NA
>> ras.645.crop[ras.645.crop>90000]<- NA
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r667 =
>> values(r667))
>> d$r667[is.na(d$r667)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r667', cell.size=0.009337697,
>> method='raster')
>> ras.667 <- raster(dd)
>> ras.667.crop<-crop(ras.667, myext)
>> ras.667.crop[ras.667.crop<0]<-NA
>> ras.667.crop[ras.667.crop>90000]<- NA
>>
>> d <- data.frame(x = values(longitude), y = values(latitude), r678 =
>> values(r678))
>> d$r678[is.na(d$r678)]<- 99999
>> coordinates(d)<-~x+y
>> dd<-vect2rast.SpatialPoints(d, fname ='r678', cell.size=0.009337697,
>> method='raster')
>> ras.678 <- raster(dd)
>> ras.678.crop<-crop(ras.678, myext)
>> ras.678.crop[ras.678.crop<0]<-NA
>> ras.678.crop[ras.678.crop>90000]<- NA
>>
>> myext <- extent(-84.9, -79.6, 42.9, 46.6) # Approx. extent that will cover
>> all of Lake Huron
>> correct.crop<-crop(correct, myext)
>> ras.412.r<-resample(ras.412.crop, correct.crop, method='bilinear')
>> ras.443.r<-resample(ras.443.crop, correct.crop, method='bilinear')
>> ras.469.r<-resample(ras.469.crop, correct.crop, method='bilinear')
>> ras.488.r<-resample(ras.488.crop, correct.crop, method='bilinear')
>> ras.531.r<-resample(ras.531.crop, correct.crop, method='bilinear')
>> ras.547.r<-resample(ras.547.crop, correct.crop, method='bilinear')
>> ras.555.r<-resample(ras.555.crop, correct.crop, method='bilinear')
>> ras.645.r<-resample(ras.645.crop, correct.crop, method='bilinear')
>> ras.667.r<-resample(ras.667.crop, correct.crop, method='bilinear')
>> ras.678.r<-resample(ras.678.crop, correct.crop, method='bilinear')
>> End<-Sys.time()
>> difftime(End, start)
>>
>> #plot original unmapped nc band Rrs_412 versus mapped version for gross
>> comparison
>> par(mfrow=c(2,1), mar=c(2,2,2,2))
>> my.orig.extent<-extent(100,475,200,600)
>> r412.crop<-crop(r412, my.orig.extent)
>> plot(r412.crop, main='Original unmapped')
>> plot(ras.412.r, main='Mapped in R')
>>
>> --
>> David Warner
>> Research Fisheries Biologist
>> U.S.G.S. Great Lakes Science Center
>> 1451 Green Road
>> Ann Arbor, MI 48105
>> 734-214-9392 <(734)%20214-9392>
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> 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
>
>

--
David Warner
Research Fisheries Biologist
U.S.G.S. Great Lakes Science Center
1451 Green Road
Ann Arbor, MI 48105
734-214-9392

        [[alternative HTML version deleted]]

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

Re: netcdf file mapping and dimensions issue

Tue, 05/23/2017 - 10:39
I tried a few pretty unpromising things and then pointed gdalwarp at it.
Seems fine, gdalwarp knows how to find geolocation arrays and use them, and
it will auto-choose an output grid, and you can provide a custom one (i.e.
with a local projection.

The SDS (subdataset) naming thing is a bit awkward until you get used to
it.  This was just GDAL 2.1.2 on Ubuntu.


f <- "
https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/A2011202183000.L2_LAC_OC.x.nc"
download.file(f, basename(f), mode = "wb")

system('gdalwarp
HDF5:"A2011202183000.L2_LAC_OC.x.nc"://geophysical_data/chl_ocx
chl_ocx .tif')
Creating output file that is 783P x 530L.
Processing input file HDF5:A2011202183000.L2_LAC_OC.x.nc
://geophysical_data/chl_ocx.
0...10...20...30...40...50...60...70...80...90...100 - done.
r <- raster("chl_ocx ")
r
class       : RasterLayer
dimensions  : 530, 783, 414990  (nrow, ncol, ncell)
resolution  : 0.01224719, 0.01224719  (x, y)
extent      : -86.99178, -77.40223, 40.73677, 47.22778  (xmin, xmax, ymin,
ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0

#plot(r)
#library(mapdata)
#map("worldHires", add = TRUE)

Cheers, Mike.


On Tue, 23 May 2017 at 23:58 Warner, David <[hidden email]> wrote:

> Greetings all
>
> I am working nc files from the Ocean Biology Procesing Group.  The files
> are MODIS aqua L2 ocean color files.
>
> These files are not mapped or perhaps some would say orthorectified such
> that when you create a rasterLayer from one of the contained variables and
> plot it, the location and orientation of features is incorrect.  For
> example, Georgian Bay, a large eastern subset of Lake Huron in the US and
> Canada, may well show up west of the lake.
>
> The dimensions of the data are read as column and row numbers rather than
> latitude and longitude.
>
> I am looking for a way to map these data in R other than my current
> method.  The current method for a single nc file is shown below.  I run
> this modified to fit into a parallelized foreach() loop that is used to
> process folders of nc files.
>
> The  basic steps that map the nc file data are 1) create a data.frame from
> the nc files measured variable(s) and the latitude and longitude then
> rasterize using vect2rast.SpatialPoints() and 2) resample the rasterLayer
> created from the data.frame using resample() to an nc file I reprojected
> using SeaDAS.  This seems overly complicated considering the fact that I
> start with a raster and it is pretty slow (30-40 seconds per nc file in
> parallel, nearly two minutes per file in the script example here).
>
> I have two questions are 1) is there a way to eliminate the used of
> vect2rast.SpatialPoints(), which is the slowest part of the process and 2)
> are there any other ways to perhaps speed the process up if I can't
> eliminate vect2rast.SpatialPoints().  I have tried rasterize() and that
> does not provide the correct product.  It recreates what I started with,
> unmapped data.  What am I missing, or is this just a function of the data I
> am using?
>
>
> Script is below.  I have included links to the two data files required in
> the process.  One is the nc file,
> https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/A2011202183000.L2_LAC_OC.x.nc
> the other is the reprojected version of this file from SeaDAS.
>
> https://cdn.rawgit.com/dmwarn/Tethys/5ce3aad2/A2011202183000.reprojected.tif
> To use them you will have to change the path to the file locations.
>
> Thanks for your patience with my writing and poorly written code.  Any
> ideas would be helpful.
>
> Dave
>
> library(raster)
> library(plotKML)
>
> rasterOptions(maxmemory = 1e+09)
> start <- Sys.time()
> #Create rasterLayers from the bands of nc file that are required for
> mapping, with r412 being the measured variable of interest
> r412<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
> var='geophysical_data/Rrs_412')
> r443<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
> var='geophysical_data/Rrs_443')
> r469<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
> var='geophysical_data/Rrs_469')
> r488<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
> var='geophysical_data/Rrs_488')
> r531<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
> var='geophysical_data/Rrs_531')
> r547<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
> var='geophysical_data/Rrs_547')
> r555<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
> var='geophysical_data/Rrs_555')
> r645<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
> var='geophysical_data/Rrs_645')
> r667<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
> var='geophysical_data/Rrs_667')
> r678<-raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
> var='geophysical_data/Rrs_678')
>
> longitude <-
> raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
> var="navigation_data/longitude")
> latitude <- raster('/Users/dmwarner/Documents/MODIS/MODIS/L2v73/Huron/2011/
> A2011202183000.L2_LAC_OC.x.nc <http://a2011202183000.l2_lac_oc.x.nc/>',
> var="navigation_data/latitude")
> r412  # shows that dimensions are row and column numbers, not lat and long
>
>
> #Import single reprojected band of the nc file as geotif, with reprojection
> done manually in SeaDAS software
> correct<-raster("/Users/dmwarner/Documents/MODIS/OC/
> correct/Huron/A2011202183000.reprojected.tif")
>
>
> #now rasterize and map the nc file after
>
> d <- data.frame(x = values(longitude), y = values(latitude), r412 =
> values(r412))
> d$r412[is.na(d$r412)]<- 99999
> coordinates(d)<-~x+y
> dd<-vect2rast.SpatialPoints(d, fname ='r412', cell.size=0.009337697,
> method='raster')
> ras.412 <- raster(dd)
> ras.412.crop<-crop(ras.412, myext)
> #ras.412.crop[ras.412.crop<0]<-NA
> ras.412.crop[ras.412.crop>90000]<- NA
>
>
> d <- data.frame(x = values(longitude), y = values(latitude), r443 =
> values(r443))
> d$r443[is.na(d$r443)]<- 99999
> coordinates(d)<-~x+y
> dd<-vect2rast.SpatialPoints(d, fname ='r443', cell.size=0.009337697,
> method='raster')
> ras.443 <- raster(dd)
> ras.443.crop<-crop(ras.443, myext)
> ras.443.crop[ras.443.crop<0]<-NA
> ras.443.crop[ras.443.crop>90000]<- NA
>
> d <- data.frame(x = values(longitude), y = values(latitude), r469 =
> values(r469))
> d$r469[is.na(d$r469)]<- 99999
> coordinates(d)<-~x+y
> dd<-vect2rast.SpatialPoints(d, fname ='r469', cell.size=0.009337697,
> method='raster')
> ras.469 <- raster(dd)
> ras.469.crop<-crop(ras.469, myext)
> ras.469.crop[ras.469.crop<0]<-NA
> ras.469.crop[ras.469.crop>90000]<- NA
>
> d <- data.frame(x = values(longitude), y = values(latitude), r488 =
> values(r488))
> d$r488[is.na(d$r488)]<- 99999
> coordinates(d)<-~x+y
> dd<-vect2rast.SpatialPoints(d, fname ='r488', cell.size=0.009337697,
> method='raster')
> ras.488 <- raster(dd)
> ras.488.crop<-crop(ras.488, myext)
> ras.488.crop[ras.488.crop<0]<-NA
> ras.488.crop[ras.488.crop>90000]<- NA
>
> d <- data.frame(x = values(longitude), y = values(latitude), r531 =
> values(r531))
> d$r531[is.na(d$r531)]<- 99999
> coordinates(d)<-~x+y
> dd<-vect2rast.SpatialPoints(d, fname ='r531', cell.size=0.009337697,
> method='raster')
> ras.531 <- raster(dd)
> ras.531.crop<-crop(ras.531, myext)
> ras.531.crop[ras.531.crop<0]<-NA
> ras.531.crop[ras.531.crop>90000]<- NA
>
> d <- data.frame(x = values(longitude), y = values(latitude), r547 =
> values(r547))
> d$r547[is.na(d$r547)]<- 99999
> coordinates(d)<-~x+y
> dd<-vect2rast.SpatialPoints(d, fname ='r547', cell.size=0.009337697,
> method='raster')
> ras.547 <- raster(dd)
> ras.547.crop<-crop(ras.547, myext)
> ras.547.crop[ras.547.crop<0]<-NA
> ras.547.crop[ras.547.crop>90000]<- NA
>
> d <- data.frame(x = values(longitude), y = values(latitude), r555 =
> values(r555))
> d$r555[is.na(d$r555)]<- 99999
> coordinates(d)<-~x+y
> dd<-vect2rast.SpatialPoints(d, fname ='r555', cell.size=0.009337697,
> method='raster')
> ras.555 <- raster(dd)
> ras.555.crop<-crop(ras.555, myext)
> ras.555.crop[ras.555.crop<0]<-NA
> ras.555.crop[ras.555.crop>90000]<- NA
>
> d <- data.frame(x = values(longitude), y = values(latitude), r645 =
> values(r645))
> d$r645[is.na(d$r645)]<- 99999
> coordinates(d)<-~x+y
> dd<-vect2rast.SpatialPoints(d, fname ='r645', cell.size=0.009337697,
> method='raster')
> ras.645 <- raster(dd)
> ras.645.crop<-crop(ras.645, myext)
> ras.645.crop[ras.645.crop<0]<-NA
> ras.645.crop[ras.645.crop>90000]<- NA
>
> d <- data.frame(x = values(longitude), y = values(latitude), r667 =
> values(r667))
> d$r667[is.na(d$r667)]<- 99999
> coordinates(d)<-~x+y
> dd<-vect2rast.SpatialPoints(d, fname ='r667', cell.size=0.009337697,
> method='raster')
> ras.667 <- raster(dd)
> ras.667.crop<-crop(ras.667, myext)
> ras.667.crop[ras.667.crop<0]<-NA
> ras.667.crop[ras.667.crop>90000]<- NA
>
> d <- data.frame(x = values(longitude), y = values(latitude), r678 =
> values(r678))
> d$r678[is.na(d$r678)]<- 99999
> coordinates(d)<-~x+y
> dd<-vect2rast.SpatialPoints(d, fname ='r678', cell.size=0.009337697,
> method='raster')
> ras.678 <- raster(dd)
> ras.678.crop<-crop(ras.678, myext)
> ras.678.crop[ras.678.crop<0]<-NA
> ras.678.crop[ras.678.crop>90000]<- NA
>
> myext <- extent(-84.9, -79.6, 42.9, 46.6) # Approx. extent that will cover
> all of Lake Huron
> correct.crop<-crop(correct, myext)
> ras.412.r<-resample(ras.412.crop, correct.crop, method='bilinear')
> ras.443.r<-resample(ras.443.crop, correct.crop, method='bilinear')
> ras.469.r<-resample(ras.469.crop, correct.crop, method='bilinear')
> ras.488.r<-resample(ras.488.crop, correct.crop, method='bilinear')
> ras.531.r<-resample(ras.531.crop, correct.crop, method='bilinear')
> ras.547.r<-resample(ras.547.crop, correct.crop, method='bilinear')
> ras.555.r<-resample(ras.555.crop, correct.crop, method='bilinear')
> ras.645.r<-resample(ras.645.crop, correct.crop, method='bilinear')
> ras.667.r<-resample(ras.667.crop, correct.crop, method='bilinear')
> ras.678.r<-resample(ras.678.crop, correct.crop, method='bilinear')
> End<-Sys.time()
> difftime(End, start)
>
> #plot original unmapped nc band Rrs_412 versus mapped version for gross
> comparison
> par(mfrow=c(2,1), mar=c(2,2,2,2))
> my.orig.extent<-extent(100,475,200,600)
> r412.crop<-crop(r412, my.orig.extent)
> plot(r412.crop, main='Original unmapped')
> plot(ras.412.r, main='Mapped in R')
>
> --
> David Warner
> Research Fisheries Biologist
> U.S.G.S. Great Lakes Science Center
> 1451 Green Road
> Ann Arbor, MI 48105
> 734-214-9392 <(734)%20214-9392>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> 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

gstat or geoR create variogram with sets of bounding boxes

Tue, 05/23/2017 - 10:23
hi,
im trying to make a variogram that consider a limit made by a bounding box,
but not just one, but dividing the area of the hole set of points in
different subzones (treatments in this case).
I have tried to do so with bounding box in geoR, but i cant do it.
Also, another way its to use anisotropy, but when try to put an alpha
parameter in geostat, doesnt make any change that without it.
bottom line, i have to perform a variogram that includes only the points in
1 direction, or, tell the variogram to limit the surface to a polygon (
several ones)

if any one can help,
thanks

--
Javier Moreira de Souza
Ingeniero Agrónomo
099 406 006

        [[alternative HTML version deleted]]

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

Pages