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 47 min ago

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

netcdf file mapping and dimensions issue

Tue, 05/23/2017 - 08:58
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

        [[alternative HTML version deleted]]

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

Re: what Manifold does that Rgdal fails? Projection transformation R

Tue, 05/23/2017 - 04:11
Appreciate the feedback, cheers

(proj4 can be willed into use for decomposed sp/sf/etc but beware of
internal transpose assumption in current version, going rogue is generally
not advisable)


On Tue, 23 May 2017, 18:48 Slawomir Kozielec <[hidden email]> wrote:

> all, Thank you for your assistance,
> I resolved the issue, though I am not sure how to 'call off' my request
> for advice.
> Basically the MapInfo was projected to BNG (27700), so rgdal tried to
> transform already projected file. On Manifold I most likely and purely by
> chance set it to initial BNG and this is how it was produced correctly (i
> have no practical experience with Manifold - I got it once upon a time with
> intention to learn but never had time to). Probably Proj4 could have worked
> as well (as you provide the initial projection), but it refuses to work
> with spatial dataframes, especially lines (as points you can just extract
> as coords and transform).
> What I did was basically to replace by substitution (not transformation)
> the proj4file attribute to BNG and then transformed it - everything went
> smoothly from then on
>
> Mike, the file format is not my choice - I am preparing an application for
> my team and shp is our internal standard. Alas our GIS team leader took few
> months off, some project to clean oceans off plastic, and his deputees did
> not know how to arrange the file for me - so i took a 'raw' file from our
> vendor. The file was in TAB format. I am not disputing whatever is better
> as I am a proper layman in the order of datums, elipsums and all other
> Latin words I would love to know:)
>
>
> On Tue, May 23, 2017 at 8:40 AM, Michael Sumner <[hidden email]>
> wrote:
>
>> Quickly, there is no need to convert to shapefile. MapInfo TAB or MIF are
>> both vastly superior to SHP and are well supported by rgdal and sf.
>>
>> That takes out one weak link. I'm happy to explore with Manifold if you
>> provide reproducible examples. Also I would ask on georeference.org too,
>> though you'll get a very disparaging (of open source) perspective from some
>> of their staff. It works both ways.
>>
>> Cheers, Mike
>>
>> On Tue, 23 May 2017, 17:13 Roger Bivand <[hidden email]> wrote:
>>
>>> On Mon, 22 May 2017, Slawomir Kozielec wrote:
>>>
>>> > Dear friends,
>>> >
>>> > pretty new to the topic so let me know if i do not follow and
>>> reproducible
>>> > examples convention. I am OK with R, but spatial data is still terra
>>> > incognita
>>>
>>> Please find out about ellipsoids and datums (sorry, not plural data).
>>> They
>>> describe the assumed shape of the world, and the relative distance of a
>>> point on the surface to an assumed 3D origin (roughly).
>>>
>>> In one of your CRS, the ellipsoid is Airy, in another where a=b it is a
>>> sphere. In addition, EPSG 4269 presupposes NAD 1983, equivalent to GRS
>>> 1980 and effectively WGS84. What you do not have are +towgs84=
>>> definitions, but if these are not given, they are assumed by proj in
>>> rgdal
>>> to be WGS84.
>>>
>>> None of these are "correct", all are based on chosen assumptions. Why
>>> Manifold is assuming a spheroid is unknown. Define your datums correctly,
>>> and things should clear up.
>>>
>>> Hope this clarifies,
>>>
>>> Roger
>>>
>>> >
>>> > I got a Mapinfo file and i had to convert it to ESRI shp and then
>>> within r
>>> > into spatial lines data frame and use with leaflet+OS maps epsg:3857
>>> >
>>> > whatever i tried within R with spTransform (rgdal) it did not work -
>>> either
>>> > was not visible or was visible with offset (approx 2 meters). Proj4
>>> did not
>>> > want to work with SLDF.
>>> > I took the file into Manifold, converted it to 3857, returned to R,
>>> > transformed it to epsg:4269 and it works perfectly OK with leaflet +
>>> OS map
>>> > 3857. But i have to automate this process to work within R, without any
>>> > Manifold intervention.
>>> >
>>> > I checked how the files differed:
>>> > i selected one item (the same of course) from the file and to make it
>>> easy
>>> > to present created a centroid for it. Then I checked projection, xy,
>>> and
>>> > latlong.
>>> >
>>> > the item created from the file without using Manifold:
>>> > proj4string
>>> >
>>> >    "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
>>> > +y_0=-100000 +ellps=airy +units=m +no_defs"
>>> >
>>> > centroid:
>>> >
>>> >    528685 181836.2
>>> >
>>> > centroid after spTransform to 4269
>>> >
>>> >    -0.1450149 51.52028
>>> >
>>> > the item from the file pre-processed by Manifold:
>>> > proj4string
>>> >
>>> >    "+proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137
>>> > +units=m +no_defs"
>>> >
>>> > centroid:
>>> >
>>> >    -16321.64 6713938
>>> >
>>> > centroid after spTransform to 4269:
>>> >
>>> >    -0.1466198 51.52079
>>> >
>>> > the latter is correct of course. Just to emphasize - it is the same
>>> file
>>> > and the same item
>>> >
>>> > Few questions:
>>> >
>>> > 1. any idea what happens in Manifold that fails in rgdal?
>>> >
>>> > 2. is there any mechanism in R that I can apply to fix the error (to
>>> match
>>> > the 'with Manifold' output)? Maybe a kind of offset or work on prj
>>> files?
>>> >
>>> > I tried to use the latter CRS for transformation
>>> >
>>> >    P5 <-spTransform(ptp.points1P, CRS( "+proj=merc +lon_0=0 +lat_ts=0
>>> > +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs"))
>>> >
>>> > and it returned still incorrect values
>>> > the xy is t least inverted but shifted
>>> >
>>> >     -16142.98 6713847
>>> >
>>> > the latlong is of course incorrect, but the same as without any
>>> attempt to
>>> > transform
>>> >
>>> >    -0.1450149 51.52028
>>> >
>>> >
>>> > content of prj file not transformed - failing to correctly show
>>> >
>>> >    PROJCS["Transverse_Mercator",GEOGCS["GCS_Airy
>>> >
>>> 1830",DATUM["D_unknown",SPHEROID["airy",6377563.396,299.3249646]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],UNIT["Meter",1]]
>>> >
>>> > content of prj manifold transformed file - good to go
>>> >
>>> >    PROJCS["Mercator_2SP",GEOGCS["GCS_unnamed
>>> >
>>> ellipse",DATUM["D_unknown",SPHEROID["Unknown",6378137,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator"],PARAMETER["standard_parallel_1",0],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]
>>> >
>>> > thanks in advance
>>> >
>>> >       [[alternative HTML version deleted]]
>>> >
>>> > _______________________________________________
>>> > R-sig-Geo mailing list
>>> > [hidden email]
>>> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>> >
>>>
>>> --
>>> Roger Bivand
>>> Department of Economics, Norwegian School of Economics,
>>> Helleveien 30, N-5045 Bergen, Norway.
>>> voice: +47 55 95 93 55 <+47%2055%2095%2093%2055>; e-mail:
>>> [hidden email]
>>> Editor-in-Chief of The R Journal,
>>> https://journal.r-project.org/index.html
>>> http://orcid.org/0000-0003-2392-6140
>>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>>
>>> _______________________________________________
>>> 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: what Manifold does that Rgdal fails? Projection transformation R

Tue, 05/23/2017 - 03:48
all, Thank you for your assistance,
I resolved the issue, though I am not sure how to 'call off' my request for
advice.
Basically the MapInfo was projected to BNG (27700), so rgdal tried to
transform already projected file. On Manifold I most likely and purely by
chance set it to initial BNG and this is how it was produced correctly (i
have no practical experience with Manifold - I got it once upon a time with
intention to learn but never had time to). Probably Proj4 could have worked
as well (as you provide the initial projection), but it refuses to work
with spatial dataframes, especially lines (as points you can just extract
as coords and transform).
What I did was basically to replace by substitution (not transformation)
the proj4file attribute to BNG and then transformed it - everything went
smoothly from then on

Mike, the file format is not my choice - I am preparing an application for
my team and shp is our internal standard. Alas our GIS team leader took few
months off, some project to clean oceans off plastic, and his deputees did
not know how to arrange the file for me - so i took a 'raw' file from our
vendor. The file was in TAB format. I am not disputing whatever is better
as I am a proper layman in the order of datums, elipsums and all other
Latin words I would love to know:)


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

> Quickly, there is no need to convert to shapefile. MapInfo TAB or MIF are
> both vastly superior to SHP and are well supported by rgdal and sf.
>
> That takes out one weak link. I'm happy to explore with Manifold if you
> provide reproducible examples. Also I would ask on georeference.org too,
> though you'll get a very disparaging (of open source) perspective from some
> of their staff. It works both ways.
>
> Cheers, Mike
>
> On Tue, 23 May 2017, 17:13 Roger Bivand <[hidden email]> wrote:
>
>> On Mon, 22 May 2017, Slawomir Kozielec wrote:
>>
>> > Dear friends,
>> >
>> > pretty new to the topic so let me know if i do not follow and
>> reproducible
>> > examples convention. I am OK with R, but spatial data is still terra
>> > incognita
>>
>> Please find out about ellipsoids and datums (sorry, not plural data). They
>> describe the assumed shape of the world, and the relative distance of a
>> point on the surface to an assumed 3D origin (roughly).
>>
>> In one of your CRS, the ellipsoid is Airy, in another where a=b it is a
>> sphere. In addition, EPSG 4269 presupposes NAD 1983, equivalent to GRS
>> 1980 and effectively WGS84. What you do not have are +towgs84=
>> definitions, but if these are not given, they are assumed by proj in rgdal
>> to be WGS84.
>>
>> None of these are "correct", all are based on chosen assumptions. Why
>> Manifold is assuming a spheroid is unknown. Define your datums correctly,
>> and things should clear up.
>>
>> Hope this clarifies,
>>
>> Roger
>>
>> >
>> > I got a Mapinfo file and i had to convert it to ESRI shp and then
>> within r
>> > into spatial lines data frame and use with leaflet+OS maps epsg:3857
>> >
>> > whatever i tried within R with spTransform (rgdal) it did not work -
>> either
>> > was not visible or was visible with offset (approx 2 meters). Proj4 did
>> not
>> > want to work with SLDF.
>> > I took the file into Manifold, converted it to 3857, returned to R,
>> > transformed it to epsg:4269 and it works perfectly OK with leaflet + OS
>> map
>> > 3857. But i have to automate this process to work within R, without any
>> > Manifold intervention.
>> >
>> > I checked how the files differed:
>> > i selected one item (the same of course) from the file and to make it
>> easy
>> > to present created a centroid for it. Then I checked projection, xy, and
>> > latlong.
>> >
>> > the item created from the file without using Manifold:
>> > proj4string
>> >
>> >    "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
>> > +y_0=-100000 +ellps=airy +units=m +no_defs"
>> >
>> > centroid:
>> >
>> >    528685 181836.2
>> >
>> > centroid after spTransform to 4269
>> >
>> >    -0.1450149 51.52028
>> >
>> > the item from the file pre-processed by Manifold:
>> > proj4string
>> >
>> >    "+proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137
>> > +units=m +no_defs"
>> >
>> > centroid:
>> >
>> >    -16321.64 6713938
>> >
>> > centroid after spTransform to 4269:
>> >
>> >    -0.1466198 51.52079
>> >
>> > the latter is correct of course. Just to emphasize - it is the same file
>> > and the same item
>> >
>> > Few questions:
>> >
>> > 1. any idea what happens in Manifold that fails in rgdal?
>> >
>> > 2. is there any mechanism in R that I can apply to fix the error (to
>> match
>> > the 'with Manifold' output)? Maybe a kind of offset or work on prj
>> files?
>> >
>> > I tried to use the latter CRS for transformation
>> >
>> >    P5 <-spTransform(ptp.points1P, CRS( "+proj=merc +lon_0=0 +lat_ts=0
>> > +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs"))
>> >
>> > and it returned still incorrect values
>> > the xy is t least inverted but shifted
>> >
>> >     -16142.98 6713847
>> >
>> > the latlong is of course incorrect, but the same as without any attempt
>> to
>> > transform
>> >
>> >    -0.1450149 51.52028
>> >
>> >
>> > content of prj file not transformed - failing to correctly show
>> >
>> >    PROJCS["Transverse_Mercator",GEOGCS["GCS_Airy
>> > 1830",DATUM["D_unknown",SPHEROID["airy",6377563.396,
>> 299.3249646]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],
>> PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],
>> PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.
>> 9996012717],PARAMETER["false_easting",400000],PARAMETER["
>> false_northing",-100000],UNIT["Meter",1]]
>> >
>> > content of prj manifold transformed file - good to go
>> >
>> >    PROJCS["Mercator_2SP",GEOGCS["GCS_unnamed
>> > ellipse",DATUM["D_unknown",SPHEROID["Unknown",6378137,0]]
>> ,PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]
>> ,PROJECTION["Mercator"],PARAMETER["standard_parallel_
>> 1",0],PARAMETER["central_meridian",0],PARAMETER["false_
>> easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]
>> >
>> > thanks in advance
>> >
>> >       [[alternative HTML version deleted]]
>> >
>> > _______________________________________________
>> > R-sig-Geo mailing list
>> > [hidden email]
>> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >
>>
>> --
>> Roger Bivand
>> Department of Economics, Norwegian School of Economics,
>> Helleveien 30, N-5045 Bergen, Norway.
>> voice: +47 55 95 93 55 <+47%2055%2095%2093%2055>; e-mail:
>> [hidden email]
>> Editor-in-Chief of The R Journal, https://journal.r-project.org/
>> index.html
>> http://orcid.org/0000-0003-2392-6140
>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>
>> _______________________________________________
>> 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: what Manifold does that Rgdal fails? Projection transformation R

Tue, 05/23/2017 - 02:40
Quickly, there is no need to convert to shapefile. MapInfo TAB or MIF are
both vastly superior to SHP and are well supported by rgdal and sf.

That takes out one weak link. I'm happy to explore with Manifold if you
provide reproducible examples. Also I would ask on georeference.org too,
though you'll get a very disparaging (of open source) perspective from some
of their staff. It works both ways.

Cheers, Mike

On Tue, 23 May 2017, 17:13 Roger Bivand <[hidden email]> wrote:

> On Mon, 22 May 2017, Slawomir Kozielec wrote:
>
> > Dear friends,
> >
> > pretty new to the topic so let me know if i do not follow and
> reproducible
> > examples convention. I am OK with R, but spatial data is still terra
> > incognita
>
> Please find out about ellipsoids and datums (sorry, not plural data). They
> describe the assumed shape of the world, and the relative distance of a
> point on the surface to an assumed 3D origin (roughly).
>
> In one of your CRS, the ellipsoid is Airy, in another where a=b it is a
> sphere. In addition, EPSG 4269 presupposes NAD 1983, equivalent to GRS
> 1980 and effectively WGS84. What you do not have are +towgs84=
> definitions, but if these are not given, they are assumed by proj in rgdal
> to be WGS84.
>
> None of these are "correct", all are based on chosen assumptions. Why
> Manifold is assuming a spheroid is unknown. Define your datums correctly,
> and things should clear up.
>
> Hope this clarifies,
>
> Roger
>
> >
> > I got a Mapinfo file and i had to convert it to ESRI shp and then within
> r
> > into spatial lines data frame and use with leaflet+OS maps epsg:3857
> >
> > whatever i tried within R with spTransform (rgdal) it did not work -
> either
> > was not visible or was visible with offset (approx 2 meters). Proj4 did
> not
> > want to work with SLDF.
> > I took the file into Manifold, converted it to 3857, returned to R,
> > transformed it to epsg:4269 and it works perfectly OK with leaflet + OS
> map
> > 3857. But i have to automate this process to work within R, without any
> > Manifold intervention.
> >
> > I checked how the files differed:
> > i selected one item (the same of course) from the file and to make it
> easy
> > to present created a centroid for it. Then I checked projection, xy, and
> > latlong.
> >
> > the item created from the file without using Manifold:
> > proj4string
> >
> >    "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
> > +y_0=-100000 +ellps=airy +units=m +no_defs"
> >
> > centroid:
> >
> >    528685 181836.2
> >
> > centroid after spTransform to 4269
> >
> >    -0.1450149 51.52028
> >
> > the item from the file pre-processed by Manifold:
> > proj4string
> >
> >    "+proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137
> > +units=m +no_defs"
> >
> > centroid:
> >
> >    -16321.64 6713938
> >
> > centroid after spTransform to 4269:
> >
> >    -0.1466198 51.52079
> >
> > the latter is correct of course. Just to emphasize - it is the same file
> > and the same item
> >
> > Few questions:
> >
> > 1. any idea what happens in Manifold that fails in rgdal?
> >
> > 2. is there any mechanism in R that I can apply to fix the error (to
> match
> > the 'with Manifold' output)? Maybe a kind of offset or work on prj files?
> >
> > I tried to use the latter CRS for transformation
> >
> >    P5 <-spTransform(ptp.points1P, CRS( "+proj=merc +lon_0=0 +lat_ts=0
> > +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs"))
> >
> > and it returned still incorrect values
> > the xy is t least inverted but shifted
> >
> >     -16142.98 6713847
> >
> > the latlong is of course incorrect, but the same as without any attempt
> to
> > transform
> >
> >    -0.1450149 51.52028
> >
> >
> > content of prj file not transformed - failing to correctly show
> >
> >    PROJCS["Transverse_Mercator",GEOGCS["GCS_Airy
> >
> 1830",DATUM["D_unknown",SPHEROID["airy",6377563.396,299.3249646]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],UNIT["Meter",1]]
> >
> > content of prj manifold transformed file - good to go
> >
> >    PROJCS["Mercator_2SP",GEOGCS["GCS_unnamed
> >
> ellipse",DATUM["D_unknown",SPHEROID["Unknown",6378137,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator"],PARAMETER["standard_parallel_1",0],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]
> >
> > thanks in advance
> >
> >       [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: [hidden email]
> Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
> http://orcid.org/0000-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>
> _______________________________________________
> 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: what Manifold does that Rgdal fails? Projection transformation R

Tue, 05/23/2017 - 02:13
On Mon, 22 May 2017, Slawomir Kozielec wrote:

> Dear friends,
>
> pretty new to the topic so let me know if i do not follow and reproducible
> examples convention. I am OK with R, but spatial data is still terra
> incognita

Please find out about ellipsoids and datums (sorry, not plural data). They
describe the assumed shape of the world, and the relative distance of a
point on the surface to an assumed 3D origin (roughly).

In one of your CRS, the ellipsoid is Airy, in another where a=b it is a
sphere. In addition, EPSG 4269 presupposes NAD 1983, equivalent to GRS
1980 and effectively WGS84. What you do not have are +towgs84=
definitions, but if these are not given, they are assumed by proj in rgdal
to be WGS84.

None of these are "correct", all are based on chosen assumptions. Why
Manifold is assuming a spheroid is unknown. Define your datums correctly,
and things should clear up.

Hope this clarifies,

Roger

>
> I got a Mapinfo file and i had to convert it to ESRI shp and then within r
> into spatial lines data frame and use with leaflet+OS maps epsg:3857
>
> whatever i tried within R with spTransform (rgdal) it did not work - either
> was not visible or was visible with offset (approx 2 meters). Proj4 did not
> want to work with SLDF.
> I took the file into Manifold, converted it to 3857, returned to R,
> transformed it to epsg:4269 and it works perfectly OK with leaflet + OS map
> 3857. But i have to automate this process to work within R, without any
> Manifold intervention.
>
> I checked how the files differed:
> i selected one item (the same of course) from the file and to make it easy
> to present created a centroid for it. Then I checked projection, xy, and
> latlong.
>
> the item created from the file without using Manifold:
> proj4string
>
>    "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
> +y_0=-100000 +ellps=airy +units=m +no_defs"
>
> centroid:
>
>    528685 181836.2
>
> centroid after spTransform to 4269
>
>    -0.1450149 51.52028
>
> the item from the file pre-processed by Manifold:
> proj4string
>
>    "+proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137
> +units=m +no_defs"
>
> centroid:
>
>    -16321.64 6713938
>
> centroid after spTransform to 4269:
>
>    -0.1466198 51.52079
>
> the latter is correct of course. Just to emphasize - it is the same file
> and the same item
>
> Few questions:
>
> 1. any idea what happens in Manifold that fails in rgdal?
>
> 2. is there any mechanism in R that I can apply to fix the error (to match
> the 'with Manifold' output)? Maybe a kind of offset or work on prj files?
>
> I tried to use the latter CRS for transformation
>
>    P5 <-spTransform(ptp.points1P, CRS( "+proj=merc +lon_0=0 +lat_ts=0
> +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs"))
>
> and it returned still incorrect values
> the xy is t least inverted but shifted
>
>     -16142.98 6713847
>
> the latlong is of course incorrect, but the same as without any attempt to
> transform
>
>    -0.1450149 51.52028
>
>
> content of prj file not transformed - failing to correctly show
>
>    PROJCS["Transverse_Mercator",GEOGCS["GCS_Airy
> 1830",DATUM["D_unknown",SPHEROID["airy",6377563.396,299.3249646]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],UNIT["Meter",1]]
>
> content of prj manifold transformed file - good to go
>
>    PROJCS["Mercator_2SP",GEOGCS["GCS_unnamed
> ellipse",DATUM["D_unknown",SPHEROID["Unknown",6378137,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator"],PARAMETER["standard_parallel_1",0],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]
>
> thanks in advance
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Pages