from:"Antonio Silva"

Re: [R-sig-Geo] extract data from multilayer nc file

2019-05-30 Thread Antonio Silva

Thanks Thiago

Good idea to change stack to brick with level command.

Best regards,

A.O.


Em qui, 30 de mai de 2019 às 17:03, Thiago V. dos Santos <
thi_vel...@yahoo.com.br> escreveu:

> Alternatively, there is a "level" argument for the brick command that
> allows accessing the 4th dimension variable of the netcdf file - usually
> height for atmosphere-related variables or depth for surface- or
> ocean-related variables.
>
> Depending on your purpose, you might want to access a specific depth (say
> the first one):
>
> rst.data <- brick('~/Downloads/
> dataset-armor-3d-rep-monthly_1559213337443.nc', varname='so', level=1)
>
> Or you could loop over the depths:
>
> # Loop over depths
> for (i in 1:2){
>   rst.data <- brick('~/Downloads/
> dataset-armor-3d-rep-monthly_1559213337443.nc', varname='so', level=i)
>
>   # Do whatever you need to do with each layer
> }
>
> See ?brick
>
> Greetings,
>  -- Thiago V. dos Santos
>
> Postdoctoral Research Fellow
> Department of Climate and Space Science and Engineering
> University of Michigan
>
>
> On Thursday, May 30, 2019, 3:42:01 PM EDT, Antonio Silva <
> aolinto@gmail.com> wrote:
>
>
> Thanks a lot Sarah. I ran your examples. Since I use Linux Ubuntu command
> lines just worked fine.
>
> I will learn more about NCO.
>
> Once more thanks for your help.
>
> All the best
>
> A.O.
>
> Em qui, 30 de mai de 2019 às 13:12, Sarah Goslee 
> escreveu:
>
> > Hi Antonio,
> >
> > I work with NetCDF files a lot, and I find it enormously easier to
> > split up multidimensional files before importing them into R.
> >
> > I use NCO for this (http://nco.sourceforge.net/) but I'm sure there
> > are other options. I'm afraid I have no idea how well that works on
> > windows.
> >
> > For your sample file, here's how you could split it by depth and by
> > variable.
> >
> > ncks -d depth,0,0 -v to dataset-armor-3d-rep-monthly_1559213337443.nc
> > vo50.nc
> > ncks -d depth,0,0 -v so dataset-armor-3d-rep-monthly_1559213337443.nc
> > so50.nc
> >
> > ncks -d depth,1,1 -v to dataset-armor-3d-rep-monthly_1559213337443.nc
> > vo75.nc
> > ncks -d depth,1,1 -v so dataset-armor-3d-rep-monthly_1559213337443.nc
> > so75.nc
> >
> > > library(raster)
> > Loading required package: sp
> > > vo75 <- stack("vo75.nc")
> > Loading required namespace: ncdf4
> > > dim(vo75)
> > [1] 29 41 12
> >
> > Sarah
> >
> > On Thu, May 30, 2019 at 10:24 AM Antonio Silva 
> > wrote:
> > >
> > > Hello all,
> > >
> > > I want to extract monthly sea water temperature and salinity (50 and
> 75 m
> > > depth) for a given area.
> > >
> > > To do so I downloaded dataset-armor-3d-rep-monthly (product
> > > MULTIOBS_GLO_PHY_REP_015_002) from Copernicus Marine environment
> > monitoring
> > > service (http://marine.copernicus.eu/). A sample file can be found at
> > > https://app.box.com/s/f8yyf3jn9kye1htqc7aktjk38jp39yrm.
> > >
> > > I also prepared a shapefile delimiting the area of interest (
> > > https://app.box.com/s/ciq76w7918zir3hu89l8d3y9dn3csuxp).
> > >
> > > Some time ago I dealt with monthly salinity data files, that is, with
> > only
> > > one variable. I prepared a script and made it available at
> > > https://gist.github.com/aolinto/59e43cdbac3ff3e0cd5eb2eab5efcbc4.
> > >
> > > Now I have a nc file with several variables (so and to) and each one
> has
> > 4
> > > dimensions (time (for each month), lat, lon and depth (50 and 75 m)).
> > >
> > > For a file with 3 dimensions per variable I write
> > > rst.data <- stack("dataset-sss-ssd-rep-monthly_sos_2008_2017.nc
> > ",varname="sos")
> > >
> > >
> > > But this turn I receive the message
> > > > rst.data <- stack("dataset-armor-3d-rep-monthly_1559213337443.nc
> > ",varname="to")
> > > # indicate the nc file
> > > Warning messages:
> > > 1: In .stackCDF(x, varname = varname, bands = bands) :
> > >  to has 4 dimensions, I am using the last one
> > > 2: In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
> > >  "level" set to 1 (there are 2 levels)
> > >
> > > I tried several commands (bands, layers) but I couldn't make a stack to
> > 50
> > > and another to 75 data.
> > >
> > > I appreciate any help.
> > >
> > > Best regards,
> > >
> > > Antônio Olinto
> > > Instituto de Pesca (Fisheries Institute)
> > > São Paulo, Brasil
> > >
> > >[[alternative HTML version deleted]]
> > >
> > > ___
> > > R-sig-Geo mailing list
> > > R-sig-Geo@r-project.org
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> >
> >
> > --
> > Sarah Goslee (she/her)
> > http://www.numberwright.com
>
> >
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] extract data from multilayer nc file

2019-05-30 Thread Antonio Silva

Thanks a lot Sarah. I ran your examples. Since I use Linux Ubuntu command
lines just worked fine.

I will learn more about NCO.

Once more thanks for your help.

All the best

A.O.

Em qui, 30 de mai de 2019 às 13:12, Sarah Goslee 
escreveu:

> Hi Antonio,
>
> I work with NetCDF files a lot, and I find it enormously easier to
> split up multidimensional files before importing them into R.
>
> I use NCO for this (http://nco.sourceforge.net/) but I'm sure there
> are other options. I'm afraid I have no idea how well that works on
> windows.
>
> For your sample file, here's how you could split it by depth and by
> variable.
>
> ncks -d depth,0,0 -v to dataset-armor-3d-rep-monthly_1559213337443.nc
> vo50.nc
> ncks -d depth,0,0 -v so dataset-armor-3d-rep-monthly_1559213337443.nc
> so50.nc
>
> ncks -d depth,1,1 -v to dataset-armor-3d-rep-monthly_1559213337443.nc
> vo75.nc
> ncks -d depth,1,1 -v so dataset-armor-3d-rep-monthly_1559213337443.nc
> so75.nc
>
> > library(raster)
> Loading required package: sp
> > vo75 <- stack("vo75.nc")
> Loading required namespace: ncdf4
> > dim(vo75)
> [1] 29 41 12
>
> Sarah
>
> On Thu, May 30, 2019 at 10:24 AM Antonio Silva 
> wrote:
> >
> > Hello all,
> >
> > I want to extract monthly sea water temperature and salinity (50 and 75 m
> > depth) for a given area.
> >
> > To do so I downloaded dataset-armor-3d-rep-monthly (product
> > MULTIOBS_GLO_PHY_REP_015_002) from Copernicus Marine environment
> monitoring
> > service (http://marine.copernicus.eu/). A sample file can be found at
> > https://app.box.com/s/f8yyf3jn9kye1htqc7aktjk38jp39yrm.
> >
> > I also prepared a shapefile delimiting the area of interest (
> > https://app.box.com/s/ciq76w7918zir3hu89l8d3y9dn3csuxp).
> >
> > Some time ago I dealt with monthly salinity data files, that is, with
> only
> > one variable. I prepared a script and made it available at
> > https://gist.github.com/aolinto/59e43cdbac3ff3e0cd5eb2eab5efcbc4.
> >
> > Now I have a nc file with several variables (so and to) and each one has
> 4
> > dimensions (time (for each month), lat, lon and depth (50 and 75 m)).
> >
> > For a file with 3 dimensions per variable I write
> > rst.data <- stack("dataset-sss-ssd-rep-monthly_sos_2008_2017.nc
> ",varname="sos")
> >
> >
> > But this turn I receive the message
> > > rst.data <- stack("dataset-armor-3d-rep-monthly_1559213337443.nc
> ",varname="to")
> > # indicate the nc file
> > Warning messages:
> > 1: In .stackCDF(x, varname = varname, bands = bands) :
> >   to has 4 dimensions, I am using the last one
> > 2: In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
> >   "level" set to 1 (there are 2 levels)
> >
> > I tried several commands (bands, layers) but I couldn't make a stack to
> 50
> > and another to 75 data.
> >
> > I appreciate any help.
> >
> > Best regards,
> >
> > Antônio Olinto
> > Instituto de Pesca (Fisheries Institute)
> > São Paulo, Brasil
> >
> > [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
> --
> Sarah Goslee (she/her)
> http://www.numberwright.com
>

[[alternative HTML version deleted]]

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

[R-sig-Geo] extract data from multilayer nc file

2019-05-30 Thread Antonio Silva

Hello all,

I want to extract monthly sea water temperature and salinity (50 and 75 m
depth) for a given area.

To do so I downloaded dataset-armor-3d-rep-monthly (product
MULTIOBS_GLO_PHY_REP_015_002) from Copernicus Marine environment monitoring
service (http://marine.copernicus.eu/). A sample file can be found at
https://app.box.com/s/f8yyf3jn9kye1htqc7aktjk38jp39yrm.

I also prepared a shapefile delimiting the area of interest (
https://app.box.com/s/ciq76w7918zir3hu89l8d3y9dn3csuxp).

Some time ago I dealt with monthly salinity data files, that is, with only
one variable. I prepared a script and made it available at
https://gist.github.com/aolinto/59e43cdbac3ff3e0cd5eb2eab5efcbc4.

Now I have a nc file with several variables (so and to) and each one has 4
dimensions (time (for each month), lat, lon and depth (50 and 75 m)).

For a file with 3 dimensions per variable I write
rst.data <- stack("dataset-sss-ssd-rep-monthly_sos_2008_2017.nc",varname="sos")


But this turn I receive the message
> rst.data <- 
> stack("dataset-armor-3d-rep-monthly_1559213337443.nc",varname="to")
# indicate the nc file
Warning messages:
1: In .stackCDF(x, varname = varname, bands = bands) :
  to has 4 dimensions, I am using the last one
2: In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
  "level" set to 1 (there are 2 levels)

I tried several commands (bands, layers) but I couldn't make a stack to 50
and another to 75 data.

I appreciate any help.

Best regards,

Antônio Olinto
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil

[[alternative HTML version deleted]]

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

[R-sig-Geo] select points out of a polygon

2019-03-26 Thread Antonio Silva

Dear list members

I have a spatial point data frame (spt.points) and a spatial polygon data
frame (spt.poly).

With spt.points[spt.poly,] I can select the points that overlap the polygon.

How can I get the points that do not overlap the polygon?

Thanks a lot.

All the best,

Antônio Olinto Ávila da Silva
São Paulo, Brasil

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] how to get an attribute from a raster layer

2018-12-28 Thread Antonio Silva

Dear Ben

Thanks again. I didn't find this function in my search.

Happy new year to all.

Antonio Olinto

Em sex, 28 de dez de 2018 às 13:19, Ben Tupper 
escreveu:

> Hi,
>
> You are looking for raster::getZ() but be aware that the Z value can be of
> any type - in this case I would guess it is a Date class so you will need
> to cast it to character using format() to get a label you desire.
>
>
> https://www.rdocumentation.org/packages/raster/versions/2.8-4/topics/z-values
>
>
> Cheers,
> Ben
>
> On Dec 28, 2018, at 10:04 AM, Antonio Silva  wrote:
>
> Dear list members
>
> I need help to do something that should be easy but I could not execute.
>
> I just want to plot a raster layer and use "z-value" as the title.
>
> rst.data[[1]]
> class   : RasterLayer
> band: 1  (of  4  bands)
> dimensions  : 33, 45, 1485  (nrow, ncol, ncell)
> resolution  : 0.25, 0.25  (x, y)
> extent  : -50.25, -39, -29, -20.75  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> data source : /mnt/.../dataset-sss-ssd-rep-monthly_2016-11-12_sos.nc
> names   : X2016.11.15.1
> z-value : 2016-11-15
> zvar: sos
>
> Is it possible to extract the z-value value directly (2016-11-15)?
>
> The best I could do was
>
> plot(rst.data[[1]],main=substr(names(rst.data[[1]]),2,11))
>
> Thanks a lot
>
> Antonio Olinto
> Fisheries Institute
> São Paulo, Brazil
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
> Ben Tupper
> Bigelow Laboratory for Ocean Sciences
> 60 Bigelow Drive, P.O. Box 380
> East Boothbay, Maine 04544
> http://www.bigelow.org
>
> Ecological Forecasting: https://eco.bigelow.org/
>

[[alternative HTML version deleted]]

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

[R-sig-Geo] how to get an attribute from a raster layer

2018-12-28 Thread Antonio Silva

Dear list members

I need help to do something that should be easy but I could not execute.

I just want to plot a raster layer and use "z-value" as the title.

rst.data[[1]]
class   : RasterLayer
band: 1  (of  4  bands)
dimensions  : 33, 45, 1485  (nrow, ncol, ncell)
resolution  : 0.25, 0.25  (x, y)
extent  : -50.25, -39, -29, -20.75  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : /mnt/.../dataset-sss-ssd-rep-monthly_2016-11-12_sos.nc
names   : X2016.11.15.1
z-value : 2016-11-15
zvar: sos

Is it possible to extract the z-value value directly (2016-11-15)?

The best I could do was

plot(rst.data[[1]],main=substr(names(rst.data[[1]]),2,11))

Thanks a lot

Antonio Olinto
Fisheries Institute
São Paulo, Brazil

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] get data from nc file

2018-12-27 Thread Antonio Silva

Thanks for you nice idea and code Ben. I will use them for sure.

Have a wonderful 2019. All the best

Antonio Olinto
Fisheries Institute
Brazil


Em qua, 26 de dez de 2018 às 22:01, Ben Tupper 
escreveu:

> Hi,
>
> .  I don't think there is any other way as the attributes are sort of
> buried in the string; that's unfortunate.  I guess you could at least make
> a reusable function assuming you'll be doing this again or looking to pull
> other attributes.  Something like this...
>
>
> #' Extract one of the GLobal Attributes of a TRMM NetCDF as a named vector
> #'
> #' @param nc the ncdf4 object
> #' @param name the name of the global attribute
> #' @param sep the separator used to delimit fields in the attribute
> #' @return named character vector of attributes
> nc_att_split <- function(nc, name = "FileHeader", sep = ";\n"){
>
> a1 <- ncdf4::ncatt_get(nc, 0)[[name[1]]]
> if (is.null(a1)) return(a1)
>
> a2 <- strsplit(a1,";\n", fixed = TRUE)[[1]]
> aa <- strsplit(a2, "=", fixed = TRUE)
>
> x <- sapply(aa,
> function(s) x = if(length(s) <= 1) "" else s[2]
> )
> names(x) <- sapply(aa,
> function(s) x = if(length(s) <= 1) "unknown" else s[1]
> )
>
> x
> }
>
>
> nc <- ncdf4::nc_open("3B43.20080101.7A.HDF.nc")
> x <- nc_att_split(nc)
> as.Date(x[['StartGranuleDateTime']], format = "%Y-%m-%dT%H:%M:%OSZ")
> [1] "2008-01-01"
>
>
> Cheers,
> Ben
>
> > On Dec 26, 2018, at 3:42 PM, Antonio Silva 
> wrote:
> >
> > Dear list members
> >
> > I downloaded some nc files with precipitation data from
> > https://pmm.nasa.gov/data-access/downloads/trmm (Level 3 3B43:
> > Multisatellite Precipitation). For the image link see the global
> attribute
> > "history" (below).
> >
> > With ncdf4::nc_open I cloud open the file (nc.data <- nc_open("
> > 3B43.20080101.7A.HDF.nc")
> >
> > I want to extract the "StartGranuleDateTime" but it is inside the global
> > attribute FileHeader (see below).
> >
> > With ncatt_get(nc.data,0,"FileHeader")$value I got
> > [1]
> >
> "AlgorithmID=3B43;\nAlgorithmVersion=3B43_7.0;\nFileName=3B43.20080101.7A.HDF;\nGenerationDateTime=2012-11-29T19:12:01.000Z;\nStartGranuleDateTime=2008-01-01T00:00:00.000Z;\nStopGranuleDateTime=2008-01-31T23:59:59.999Z;\nGranuleNumber=;\nNumberOfSwaths=0;\nNumberOfGrids=1;\nGranuleStart=;\nTimeInterval=MONTH;\nProcessingSystem=PPS;\nProductVersion=7A;\nMissingData=;\n"
> >
> > Is there any way to extract only the string "2008-01-01T00:00:00.000Z"?
> >
> > The best I could do was
> >
> as.Date(substr(strsplit(ncatt_get(nc.data,0,"FileHeader")$value,";\n")[[1]][5],22,45),"%Y-%m-%dT%H:%M:%OSZ")
> >
> > but probably, I suppose, there must be a more direct way of getting the
> > data. I appreciate any suggestions.
> >
> > Best regards,
> >
> > Antonio Olinto
> > Fisheries Institute
> > Brazil
> >
> > nc.data
> > File 3B43.20080101.7A.HDF.nc (NC_FORMAT_CLASSIC):
> >
> > 1 variables (excluding dimension variables):
> >float precipitation[nlat,nlon]
> >units: mm/hr
> >coordinates: nlon nlat
> >_FillValue: -.900390625
> >
> > 2 dimensions:
> >nlon  Size:33
> >long_name: longitude
> >standard_name: longitude
> >units: degrees_east
> >nlat  Size:41
> >long_name: latitude
> >standard_name: latitude
> >units: degrees_north
> >
> >5 global attributes:
> >Grid.GridHeader: BinMethod=ARITHMETIC_MEAN;
> > Registration=CENTER;
> > LatitudeResolution=0.25;
> > LongitudeResolution=0.25;
> > NorthBoundingCoordinate=50;
> > SouthBoundingCoordinate=-50;
> > EastBoundingCoordinate=180;
> > WestBoundingCoordinate=-180;
> > Origin=SOUTHWEST;
> >
> >FileHeader: AlgorithmID=3B43;
> > AlgorithmVersion=3B43_7.0;
> > FileName=3B43.20080101.7A.HDF;
> > GenerationDateTime=2012-11-29T19:12:01.000Z;
> > StartGranuleDateTime=2008-01-01T00:00:00.000Z;
> > StopGranuleDateTime=2008-01-31T23:59:59.999Z;
> > GranuleNumber=;
> > NumberOfSwaths=0;
> > NumberOfGrids=1;
> > GranuleStart=;
> > TimeInterval=MONTH;
> > ProcessingSystem=PPS;
>

[R-sig-Geo] get data from nc file

2018-12-26 Thread Antonio Silva

Dear list members

I downloaded some nc files with precipitation data from
https://pmm.nasa.gov/data-access/downloads/trmm (Level 3 3B43:
Multisatellite Precipitation). For the image link see the global attribute
"history" (below).

With ncdf4::nc_open I cloud open the file (nc.data <- nc_open("
3B43.20080101.7A.HDF.nc")

I want to extract the "StartGranuleDateTime" but it is inside the global
attribute FileHeader (see below).

With ncatt_get(nc.data,0,"FileHeader")$value I got
[1]
"AlgorithmID=3B43;\nAlgorithmVersion=3B43_7.0;\nFileName=3B43.20080101.7A.HDF;\nGenerationDateTime=2012-11-29T19:12:01.000Z;\nStartGranuleDateTime=2008-01-01T00:00:00.000Z;\nStopGranuleDateTime=2008-01-31T23:59:59.999Z;\nGranuleNumber=;\nNumberOfSwaths=0;\nNumberOfGrids=1;\nGranuleStart=;\nTimeInterval=MONTH;\nProcessingSystem=PPS;\nProductVersion=7A;\nMissingData=;\n"

Is there any way to extract only the string "2008-01-01T00:00:00.000Z"?

The best I could do was
as.Date(substr(strsplit(ncatt_get(nc.data,0,"FileHeader")$value,";\n")[[1]][5],22,45),"%Y-%m-%dT%H:%M:%OSZ")

but probably, I suppose, there must be a more direct way of getting the
data. I appreciate any suggestions.

Best regards,

Antonio Olinto
Fisheries Institute
Brazil

nc.data
File 3B43.20080101.7A.HDF.nc (NC_FORMAT_CLASSIC):

 1 variables (excluding dimension variables):
float precipitation[nlat,nlon]
units: mm/hr
coordinates: nlon nlat
_FillValue: -.900390625

 2 dimensions:
nlon  Size:33
long_name: longitude
standard_name: longitude
units: degrees_east
nlat  Size:41
long_name: latitude
standard_name: latitude
units: degrees_north

5 global attributes:
Grid.GridHeader: BinMethod=ARITHMETIC_MEAN;
Registration=CENTER;
LatitudeResolution=0.25;
LongitudeResolution=0.25;
NorthBoundingCoordinate=50;
SouthBoundingCoordinate=-50;
EastBoundingCoordinate=180;
WestBoundingCoordinate=-180;
Origin=SOUTHWEST;

FileHeader: AlgorithmID=3B43;
AlgorithmVersion=3B43_7.0;
FileName=3B43.20080101.7A.HDF;
GenerationDateTime=2012-11-29T19:12:01.000Z;
StartGranuleDateTime=2008-01-01T00:00:00.000Z;
StopGranuleDateTime=2008-01-31T23:59:59.999Z;
GranuleNumber=;
NumberOfSwaths=0;
NumberOfGrids=1;
GranuleStart=;
TimeInterval=MONTH;
ProcessingSystem=PPS;
ProductVersion=7A;
MissingData=;

FileInfo: DataFormatVersion=m;
TKCodeBuildVersion=1;
MetadataVersion=m;
FormatPackage=HDF Version 4.2 Release 4, January 25, 2009;
BlueprintFilename=TRMM.V7.3B43.blueprint.xml;
BlueprintVersion=BV_13;
TKIOVersion=1.6;
MetadataStyle=PVL;
EndianType=LITTLE_ENDIAN;

GridHeader: BinMethod=ARITHMETIC_MEAN;
Registration=CENTER;
LatitudeResolution=0.25;
LongitudeResolution=0.25;
NorthBoundingCoordinate=50;
SouthBoundingCoordinate=-50;
EastBoundingCoordinate=180;
WestBoundingCoordinate=-180;
Origin=SOUTHWEST;

history: 2018-12-26 17:57:56 GMT Hyrax-1.13.4
https://disc2.gesdisc.eosdis.nasa.gov:443/opendap/TRMM_L3/TRMM_3B43.7/2008/3B43.20080101.7A.HDF.nc?precipitation[604:636][3:43],nlat[3:43],nlon[604:636]

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] subsetting a spatial polygons

2018-09-12 Thread Antonio Silva

Thanks Michael and Tim, thanks for the attention

The first thing I did was I try the subsets using subset.

subset(polys,
coordinates(polys)[,1]>=-46 & coordinates(polys)[,1]<=-45 &
coordinates(polys)[,2]>=-25 & coordinates(polys)[,2]<=-24)

But it seems to work only with Spatial*DataFrame objects. Then I (unluckily)
gave up the x-y approach.

With Mike's advice I tried again as

plot(
polys[
coordinates(polys)[,1]>=-46 & coordinates(polys)[,1]<=-45 &
coordinates(polys)[,2]>=-25 & coordinates(polys)[,2]<=-24,]
,add=T,col="green")

and I got what I wanted in a more "elegant" fashion.

Once more thanks a lot for the attention.

-- 
Antônio Olinto Ávila da Silva
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil

Em qua, 12 de set de 2018 às 05:05, Michael Sumner 
escreveu:

> crop doesn't select it actually cuts on the boundary, for simple shapes I
> use coordinates (for centroids in easy matrix) and select with [ using
> tests on x and y.
>
> Cheers, Mike
>
> On Wed, Sep 12, 2018, 15:29 Tim Appelhans  wrote:
>
> > Antonio,
> >
> > I am not sure why you think that your solution is not very elegant.
> > In case you want to have more visual control over the subsetting, you
> > could try mapedit:
> >
> > library(mapedit)
> > myselection = selectFeatures(polys, mode = "draw")
> >
> > which will let you draw a e.g. rectangle and only return those features
> > that intersect it.
> >
> > Best
> > Tim
> >
> >
> > On 09/12/2018 04:18 AM, Antonio Silva wrote:
> > > Dear list users
> > >
> > > I have a SpatialPolygons with several squares. How to subset it to have
> > > only the squares between given latitudes and longitudes?
> > >
> > > In the example
> > >
> > > library(sp)
> > > library(raster)
> > > grd <-
> > >
> >
> GridTopology(cellcentre.offset=c(-47.75,-25.416667),cellsize=c(10/60,10/60),cells.dim=c(23,12))
> > > polys <- as.SpatialPolygons.GridTopology(grd)
> > > proj4string(polys) <- CRS("+proj=longlat +datum=WGS84")
> > > plot(polys,axes=T)
> > >
> > > How to select only the squares, let's say, between 24-25°S and 45-46°W?
> > >
> > > The farthest I went was:
> > >
> > > e <- extent(-45.9,-45.1,-24.9,-24.1) # which is not very elegant
> > > mask <- crop(polys,e)
> > > polys2 <- polys[mask,]
> > > plot(polys2,add=T,col="green")
> > >
> > > Thanks a lot. Best regards
> > >
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> --
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> Kingston Tasmania 7050 Australia
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

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

[R-sig-Geo] subsetting a spatial polygons

2018-09-11 Thread Antonio Silva

Dear list users

I have a SpatialPolygons with several squares. How to subset it to have
only the squares between given latitudes and longitudes?

In the example

library(sp)
library(raster)
grd <-
GridTopology(cellcentre.offset=c(-47.75,-25.416667),cellsize=c(10/60,10/60),cells.dim=c(23,12))
polys <- as.SpatialPolygons.GridTopology(grd)
proj4string(polys) <- CRS("+proj=longlat +datum=WGS84")
plot(polys,axes=T)

How to select only the squares, let's say, between 24-25°S and 45-46°W?

The farthest I went was:

e <- extent(-45.9,-45.1,-24.9,-24.1) # which is not very elegant
mask <- crop(polys,e)
polys2 <- polys[mask,]
plot(polys2,add=T,col="green")

Thanks a lot. Best regards

-- 
Antônio Olinto Ávila da Silva
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] drawing a polygon from several lines

2018-08-24 Thread Antonio Silva

It's my pleasure to share Lulla. Thanks for the suggestion, I will adopt it.

Regards,

Antonio

Em sex, 24 de ago de 2018 às 10:58, Vijay Lulla 
escreveu:

> Very nice, Antonio!!  Thanks for sharing your script.  The only minor
> suggestion I have to offer is that you can make your "cut isobath" portion
> easier to read (at least IMO) by using `with` and `between` (either
> dplyr::between or data.table::between).  Below is how I would've done it.
>
> limmin <- with(list(coords = coordinates(isomin)[[1]][[1]]),
>subset(coords, between(coords[, 1], lonmin, lonmax) &
>   between(coords[, 2], latmin, latmax)))
>
> Again, thanks for sharing your problem and solution.  I found it
> instructive!
> Cordially,
> Vijay.
>
>
> On Fri, Aug 24, 2018 at 9:16 AM Antonio Silva 
> wrote:
>
>> Dear colleagues
>>
>> I finally could select the squares (cells) under the polygon. I extract
>> the
>> points from the selected isobaths that fell within lat / lon limits and
>> with them I created a spatial polygon.
>>
>> Bellow I share the script I wrote. Considering that I don't know much
>> about
>> mapping in R, I would appreciate to hear suggestions to improve the code.
>>
>> Best regards.
>>
>> --
>> Antônio Olinto Ávila da Silva
>> Instituto de Pesca (Fisheries Institute)
>> São Paulo, Brasil
>>
>> ===
>>
>> library(rgdal)
>> library(rgeos)
>> rm(list = ls())
>>
>> # import shapes
>> # isobath
>> https://www.dropbox.com/sh/us859uyijns36r9/AAC9N4OSwweLuW-fdsBWUy-ra?dl=0
>> isobaths <- readOGR(".","isobath")
>> isobaths <- spTransform(isobaths, CRS("+proj=longlat +datum=WGS84"))
>> proj4string(isobaths)
>> summary(isobaths)
>>
>> # create a spatial grid and polygons
>> grd <-
>>
>> GridTopology(cellcentre.offset=c(-47.75,-25.416667),cellsize=c(10/60,10/60),cells.dim=c(23,12))
>> squares <- as.SpatialPolygons.GridTopology(grd)
>> proj4string(squares) <- CRS("+proj=longlat +datum=WGS84")
>> summary(squares)
>> IDs <- sapply(slot(squares, "polygons"), function(x) slot(x, "ID"))
>>
>> # set limits
>> latmin <- -24.8
>> latmax <- -24.02
>> lonmin <- -47.2
>> lonmax <- -44.8
>> depmin <- -25
>> depmax <- -60
>>
>> # bounding box
>> mat.area <-
>>
>> matrix(c(lonmin,latmax,lonmax,latmax,lonmax,latmin,lonmin,latmin),ncol=2,byrow=T)
>> spp.area <- SpatialPolygons(list(Polygons(list(Polygon(mat.area)) ,
>> ID='1')))
>> proj4string(spp.area) <- CRS("+proj=longlat +datum=WGS84")
>>
>> # plot polygons isolines and bounding box
>> plot(squares, axes=T)
>> #~ text(coordinates(squares),labels=IDs, cex=0.7)
>> plot(isobaths,add=T)
>> plot(spp.area,border="red",add=T)
>>
>> # select isobaths
>> summary(isobaths)
>> isomin <- subset(isobaths, ID %in% depmin)
>> proj4string(isomin) <- CRS("+proj=longlat +datum=WGS84")
>> isomax <- subset(isobaths, ID %in% depmax)
>> proj4string(isomax) <- CRS("+proj=longlat +datum=WGS84")
>>
>> plot(isomin,add=T,col="deepskyblue3",lwd=2)
>> plot(isomax,add=T,col="dodgerblue4",lwd=2)
>>
>> # cut isobaths
>> limmin <- subset(coordinates(isomin)[[1]][[1]],
>>   coordinates(isomin)[[1]][[1]][,2]<=latmax &
>> coordinates(isomin)[[1]][[1]][,2]>=latmin &
>>   coordinates(isomin)[[1]][[1]][,1]<=lonmax &
>> coordinates(isomin)[[1]][[1]][,1]>=lonmin)
>> limmax <- subset(coordinates(isomax)[[1]][[1]],
>>   coordinates(isomax)[[1]][[1]][,2]<=latmax &
>> coordinates(isomax)[[1]][[1]][,2]>=latmin &
>>   coordinates(isomax)[[1]][[1]][,1]<=lonmax &
>> coordinates(isomax)[[1]][[1]][,1]>=lonmin)
>>
>> points(limmin,col="chartreuse3")
>> points(limmax,col="chartreuse4")
>>
>> # create the polygon
>> spp.farea <-
>>
>> SpatialPolygons(list(Polygons(list(Polygon(rbind(limmin,limmax))),ID='1')))
>> proj4string(spp.farea) <- CRS("+proj=longlat +datum=WGS84")
>>
>> plot(spp.farea,col="chocolate3",lwd=2,add=T)
>>
>> # select the squares under the polygon
>> fareasq <- squares[spp.farea,]
>> fareace <- gCentroid(fareasq,byid=TRUE)
>>
>> # final plot
>> plot(fareasq,axes=T)
>> points(fareace,pch=10,col="darkgreen",cex=4)
>> text(coordinates(squares), labels=IDs, cex=0.7)
>> plot(isomin,add=T,col="deepskyblue3",lwd=2)
>> plot(isomax,add=T,col="dodgerblue4",lwd=2)
>> plot(spp.area,border="red",add=T,lty=2)
>>
>> [[alternative HTML version deleted]]
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
> --
> Vijay Lulla
>
> Assistant Professor,
> Dept. of Geography, IUPUI
> 425 University Blvd, CA-207C.
> Indianapolis, IN-46202
> vlu...@iupui.edu
> ORCID: https://orcid.org/-0002-0823-2522
> Webpage: http://vijaylulla.com
>

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] drawing a polygon from several lines

2018-08-24 Thread Antonio Silva

Dear colleagues

I finally could select the squares (cells) under the polygon. I extract the
points from the selected isobaths that fell within lat / lon limits and
with them I created a spatial polygon.

Bellow I share the script I wrote. Considering that I don't know much about
mapping in R, I would appreciate to hear suggestions to improve the code.

Best regards.

-- 
Antônio Olinto Ávila da Silva
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil

===

library(rgdal)
library(rgeos)
rm(list = ls())

# import shapes
# isobath
https://www.dropbox.com/sh/us859uyijns36r9/AAC9N4OSwweLuW-fdsBWUy-ra?dl=0
isobaths <- readOGR(".","isobath")
isobaths <- spTransform(isobaths, CRS("+proj=longlat +datum=WGS84"))
proj4string(isobaths)
summary(isobaths)

# create a spatial grid and polygons
grd <-
GridTopology(cellcentre.offset=c(-47.75,-25.416667),cellsize=c(10/60,10/60),cells.dim=c(23,12))
squares <- as.SpatialPolygons.GridTopology(grd)
proj4string(squares) <- CRS("+proj=longlat +datum=WGS84")
summary(squares)
IDs <- sapply(slot(squares, "polygons"), function(x) slot(x, "ID"))

# set limits
latmin <- -24.8
latmax <- -24.02
lonmin <- -47.2
lonmax <- -44.8
depmin <- -25
depmax <- -60

# bounding box
mat.area <-
matrix(c(lonmin,latmax,lonmax,latmax,lonmax,latmin,lonmin,latmin),ncol=2,byrow=T)
spp.area <- SpatialPolygons(list(Polygons(list(Polygon(mat.area)) ,
ID='1')))
proj4string(spp.area) <- CRS("+proj=longlat +datum=WGS84")

# plot polygons isolines and bounding box
plot(squares, axes=T)
#~ text(coordinates(squares),labels=IDs, cex=0.7)
plot(isobaths,add=T)
plot(spp.area,border="red",add=T)

# select isobaths
summary(isobaths)
isomin <- subset(isobaths, ID %in% depmin)
proj4string(isomin) <- CRS("+proj=longlat +datum=WGS84")
isomax <- subset(isobaths, ID %in% depmax)
proj4string(isomax) <- CRS("+proj=longlat +datum=WGS84")

plot(isomin,add=T,col="deepskyblue3",lwd=2)
plot(isomax,add=T,col="dodgerblue4",lwd=2)

# cut isobaths
limmin <- subset(coordinates(isomin)[[1]][[1]],
  coordinates(isomin)[[1]][[1]][,2]<=latmax &
coordinates(isomin)[[1]][[1]][,2]>=latmin &
  coordinates(isomin)[[1]][[1]][,1]<=lonmax &
coordinates(isomin)[[1]][[1]][,1]>=lonmin)
limmax <- subset(coordinates(isomax)[[1]][[1]],
  coordinates(isomax)[[1]][[1]][,2]<=latmax &
coordinates(isomax)[[1]][[1]][,2]>=latmin &
  coordinates(isomax)[[1]][[1]][,1]<=lonmax &
coordinates(isomax)[[1]][[1]][,1]>=lonmin)

points(limmin,col="chartreuse3")
points(limmax,col="chartreuse4")

# create the polygon
spp.farea <-
SpatialPolygons(list(Polygons(list(Polygon(rbind(limmin,limmax))),ID='1')))
proj4string(spp.farea) <- CRS("+proj=longlat +datum=WGS84")

plot(spp.farea,col="chocolate3",lwd=2,add=T)

# select the squares under the polygon
fareasq <- squares[spp.farea,]
fareace <- gCentroid(fareasq,byid=TRUE)

# final plot
plot(fareasq,axes=T)
points(fareace,pch=10,col="darkgreen",cex=4)
text(coordinates(squares), labels=IDs, cex=0.7)
plot(isomin,add=T,col="deepskyblue3",lwd=2)
plot(isomax,add=T,col="dodgerblue4",lwd=2)
plot(spp.area,border="red",add=T,lty=2)

[[alternative HTML version deleted]]

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

[R-sig-Geo] drawing a polygon from several lines

2018-08-22 Thread Antonio Silva

Hello all,

I want to select the cells of a grid that in a range of latitudes and
depths.
Latitudes limits are indicated by two lines and the depths by two isobaths
from a shapefile.
The shapefile can be downloaded at
https://www.dropbox.com/sh/us859uyijns36r9/AAC9N4OSwweLuW-fdsBWUy-ra?dl=0

Below are the steps I've already taken. I could not find a way to create a
polygon from the lines of depth and latitude. With this polygon I could
easily select the overlapped cells.

I would really appreciate any help. Best regards


library(rgeos)
library(sp)
rm(list = ls())

# import shape
isobaths <- readOGR(".","isobath")
isobaths <- spTransform(bath, CRS("+proj=longlat +datum=WGS84"))
proj4string(isobaths)
summary(isobaths)

# create a spatial grid and polygons
grd <-
GridTopology(cellcentre.offset=c(-47.75,-25.416667),cellsize=c(10/60,10/60),cells.dim=c(23,12))
polys <- as.SpatialPolygons.GridTopology(grd)
proj4string(polys) <- CRS("+proj=longlat +datum=WGS84")
summary(polys)
IDs <- sapply(slot(polys, "polygons"), function(x) slot(x, "ID"))

# create spatial lines
# southern limit
x <- c(-47.5,-45.5)
y <- c(-24.8,-24.8)
ls <- SpatialLines(list(Lines(Line(cbind(x,y)),
ID="b")),proj4string=CRS("+proj=longlat +datum=WGS84"))

# northern limit
x <- c(-46.2,-44.5)
y <- c(-24.02,-24.02)
ln <- SpatialLines(list(Lines(Line(cbind(x,y)),
ID="c")),proj4string=CRS("+proj=longlat +datum=WGS84"))

# plot isolines polygons and lines
plot(polys)
axis(1);axis(4)
text(coordinates(polys), labels=IDs, cex=0.7)
plot(isobaths,add=T)
plot(ls,add=T,col="blue",lwd=2)
plot(ln,add=T,col="orange",lwd=2)

# select isobaths
summary(isobaths)
dmin <- subset(isobaths, ID %in% "-25")
proj4string(dmin) <- CRS("+proj=longlat +datum=WGS84")

dmax <- subset(isobaths, ID %in% "-60")
proj4string(dmax) <- CRS("+proj=longlat +datum=WGS84")

plot(dmin,add=T,col="red",lwd=2)
plot(dmax,add=T,col="red",lwd=2)

# from here I'd like to have a polygon to select the cells.
na.omit(over(polys,dmin))

-- 
Antonio Olinto Ávila da Silva
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] how to create several polygons from a list of vertices

2018-08-15 Thread Antonio Silva

Hi Vijay, Michael and Obrl

Thanks for the answers, finally I could do the job.

I learned a lot too!

Best regards

Antonio Olinto

Em Ter, 14 de ago de 2018 11:19 PM, Vijay Lulla 
escreveu:

> Maybe you can try this then?
>
> polys <- SpatialPolygons(lapply(unique(vertices$cod),
>  function(x) {
>Polygons(list(Polygon(vertices[vertices$cod ==
> x, 1:2])), ID=x)
>  } ))
> HTH,
> Vijay.
>
> On Tue, Aug 14, 2018 at 6:03 PM Antonio Silva 
> wrote:
>
>> Thanks Lulla,
>>
>> Nice solution. I could also export it as a shapefile after transforming
>> it to a spatial polygon dataframe.
>>
>> The problem is that I could not "individualize" the squares in a
>> multipart layer. They all have the same ID. I tried to change this without
>> success: "Single ID required".
>>
>> The attribute table of the shapefile should have 6 lines in my example
>> and not only one.
>>
>> Any other option?
>>
>> Thanks again,
>>
>> Antonio Olinto
>>
>>
>> Em ter, 14 de ago de 2018 às 18:10, Vijay Lulla 
>> escreveu:
>>
>>> Maybe something like this?
>>>
>>> poly <- SpatialPolygons(list(Polygons(tapply(seq_len(nrow(vertices)),
>>>  vertices$cod,
>>>          function(x)
>>> Polygon(vertices[x,1:2])), ID="1")),
>>> proj4string=CRS("+proj=longlat +ellps=WGS84
>>> +datum=WGS84 +no_defs"))
>>>
>>>
>>> On Tue, Aug 14, 2018 at 4:17 PM Antonio Silva 
>>> wrote:
>>>
>>>> Hi,
>>>>
>>>> I have a data.frame with the vertices (lon / lat) and codes from several
>>>> squares (more than 500 in the real dataset).
>>>> I want to create an object with these polygons (squares) and after this
>>>> export it as a shapefile.
>>>> With the script below I can draw one square.
>>>> library(sp)
>>>> P1 = Polygon(vertices[1:4,1:2])
>>>> Ps1 = SpatialPolygons(list(Polygons(list(P1), ID = "1")),
>>>> proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
>>>> plot(Ps1, axes = TRUE)
>>>>
>>>> Now I'm trying to create one object with all squares at once.
>>>> Is it possible?
>>>>
>>>> Thanks a lot,
>>>>
>>>> Antônio Olinto
>>>>
>>>> sample data:vertices
>>>>lon lat cod
>>>> 1  -33 -23   1
>>>> 2  -32 -23   1
>>>> 3  -32 -22   1
>>>> 4  -33 -22   1
>>>> 5  -32 -23   2
>>>> 6  -31 -23   2
>>>> 7  -31 -22   2
>>>> 8  -32 -22   2
>>>> 9  -31 -23   3
>>>> 10 -30 -23   3
>>>> 11 -30 -22   3
>>>> 12 -31 -22   3
>>>> 13 -33 -22   4
>>>> 14 -32 -22   4
>>>> 15 -32 -21   4
>>>> 16 -33 -21   4
>>>> 17 -32 -22   5
>>>> 18 -31 -22   5
>>>> 19 -31 -21   5
>>>> 20 -32 -21   5
>>>> 21 -31 -22   6
>>>> 22 -30 -22   6
>>>> 23 -30 -21   6
>>>> 24 -31 -21   6
>>>>
>>>> [[alternative HTML version deleted]]
>>>>
>>>> ___
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo@r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>
> --
> Vijay Lulla
>
> Assistant Professor,
> Dept. of Geography, IUPUI
> 425 University Blvd, CA-207C.
> Indianapolis, IN-46202
> vlu...@iupui.edu
> ORCID: https://orcid.org/-0002-0823-2522
> Webpage: http://vijaylulla.com
>

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] how to create several polygons from a list of vertices

2018-08-14 Thread Antonio Silva

Thanks Lulla,

Nice solution. I could also export it as a shapefile after transforming it
to a spatial polygon dataframe.

The problem is that I could not "individualize" the squares in a multipart
layer. They all have the same ID. I tried to change this without success:
"Single ID required".

The attribute table of the shapefile should have 6 lines in my example and
not only one.

Any other option?

Thanks again,

Antonio Olinto


Em ter, 14 de ago de 2018 às 18:10, Vijay Lulla 
escreveu:

> Maybe something like this?
>
> poly <- SpatialPolygons(list(Polygons(tapply(seq_len(nrow(vertices)),
>  vertices$cod,
>  function(x)
> Polygon(vertices[x,1:2])), ID="1")),
> proj4string=CRS("+proj=longlat +ellps=WGS84
> +datum=WGS84 +no_defs"))
>
>
> On Tue, Aug 14, 2018 at 4:17 PM Antonio Silva 
> wrote:
>
>> Hi,
>>
>> I have a data.frame with the vertices (lon / lat) and codes from several
>> squares (more than 500 in the real dataset).
>> I want to create an object with these polygons (squares) and after this
>> export it as a shapefile.
>> With the script below I can draw one square.
>> library(sp)
>> P1 = Polygon(vertices[1:4,1:2])
>> Ps1 = SpatialPolygons(list(Polygons(list(P1), ID = "1")),
>> proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
>> plot(Ps1, axes = TRUE)
>>
>> Now I'm trying to create one object with all squares at once.
>> Is it possible?
>>
>> Thanks a lot,
>>
>> Antônio Olinto
>>
>> sample data:vertices
>>lon lat cod
>> 1  -33 -23   1
>> 2  -32 -23   1
>> 3  -32 -22   1
>> 4  -33 -22   1
>> 5  -32 -23   2
>> 6  -31 -23   2
>> 7  -31 -22   2
>> 8  -32 -22   2
>> 9  -31 -23   3
>> 10 -30 -23   3
>> 11 -30 -22   3
>> 12 -31 -22   3
>> 13 -33 -22   4
>> 14 -32 -22   4
>> 15 -32 -21   4
>> 16 -33 -21   4
>> 17 -32 -22   5
>> 18 -31 -22   5
>> 19 -31 -21   5
>> 20 -32 -21   5
>> 21 -31 -22   6
>> 22 -30 -22   6
>> 23 -30 -21   6
>> 24 -31 -21   6
>>
>> [[alternative HTML version deleted]]
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>

[[alternative HTML version deleted]]

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

[R-sig-Geo] how to create several polygons from a list of vertices

2018-08-14 Thread Antonio Silva

Hi,

I have a data.frame with the vertices (lon / lat) and codes from several
squares (more than 500 in the real dataset).
I want to create an object with these polygons (squares) and after this
export it as a shapefile.
With the script below I can draw one square.
library(sp)
P1 = Polygon(vertices[1:4,1:2])
Ps1 = SpatialPolygons(list(Polygons(list(P1), ID = "1")),
proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
plot(Ps1, axes = TRUE)

Now I'm trying to create one object with all squares at once.
Is it possible?

Thanks a lot,

Antônio Olinto

sample data:vertices
   lon lat cod
1  -33 -23   1
2  -32 -23   1
3  -32 -22   1
4  -33 -22   1
5  -32 -23   2
6  -31 -23   2
7  -31 -22   2
8  -32 -22   2
9  -31 -23   3
10 -30 -23   3
11 -30 -22   3
12 -31 -22   3
13 -33 -22   4
14 -32 -22   4
15 -32 -21   4
16 -33 -21   4
17 -32 -22   5
18 -31 -22   5
19 -31 -21   5
20 -32 -21   5
21 -31 -22   6
22 -30 -22   6
23 -30 -21   6
24 -31 -21   6

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] make a raster from Aquarius files

2017-12-01 Thread Antonio Silva

Thanks Ben and Michael for the attention

The file I'm trying to rasterize can be downloaded at
https://oceancolor.gsfc.nasa.gov/cgi/l3

I selected the Aquarius sea surface salinity smoothed file from June 2015
(SMI HDF).

As Ben pointed it is also available from OBPG Nasa Ocean Color site at
https://oceandata.sci.gsfc.nasa.gov/Aquarius/Mapped/Monthly/1deg/V5.0_SSS/

or directly from
https://oceandata.sci.gsfc.nasa.gov/cgi/getfile/Q20151522015181.L3m_MO_SCISM_V5.0_SSS_1deg.bz2

Thanks again

Antonio Olinto Avila da Silva
Fisheries Institute
São Paulo, Brasil


2017-12-01 12:35 GMT-02:00 Ben Tupper <btup...@bigelow.org>:

> Hi,
>
> Those Aquarius files look quite different from others I have used from
> OBPG (mostly MODISA and SeaWiFS).
>
> As a short-cut alternative, you could read all of the values into a matrix
> and make a global raster from that.  An example is shown below.  You can
> also subset the extraction as you have shown, but it might be easier to
> subset after making the raster using raster::crop()
>
> CHeers,
> Ben
>
>
> ### START
>
> library(raster)
> library(ncdf4)
> library(rasterVis)
>
> filename = "Q20151522015181.L3m_MO_SCIA_V5.0_SSS_1deg"
>
> nc = ncdf4::nc_open(filename)
> m = ncdf4::ncvar_get(nc, 'l3m_data')
> ncdf4::nc_close(nc)
> r = raster::raster(t(m))
> rasterVis::levelplot(r)
>
> ### END
>
> > On Dec 1, 2017, at 6:24 AM, Michael Sumner <mdsum...@gmail.com> wrote:
> >
> > That file name does not correspond to the standard patterns used by the
> > oceancolor site. All the L3m products from there are now (NetCDF 4.0) .nc
> > and so will work fine with raster/ncdf4.  (Some years ago they were HDF4
> -
> > without an extension, as the shortcuts in the image thumbnails hints
> > (SMI/HDF and BIN/HDF - SMI/L3m standard mapped image in your case).
> >
> > I think you've got some other provider's version of a file, but there's
> not
> > enough information here to know where you got it or what form it's in.
> I'm
> > happy to look if you can point us to the source of
> > Q20151522015181.L3m_MO_SCISM_V5.0_SSS_1deg.
> >
> > But otherwise, can you share with us the output of
> >
> > nc.data<-nc_open("Q20151522015181.L3m_MO_SCISM_V5.0_SSS_1deg")
> > print(nc.data)
> >
> > and if you're on a suitable system with HDF4 support a gdalinfo output of
> > the file would be useful too.
> >
> > Given that you can read it with ncdf4, and if it actually is NetCDF4 (not
> > HDF4 or something else) you might help raster work with it by renaming it
> > to "Q20151522015181.L3m_MO_SCISM_V5.0_SSS_1deg.nc" since (unlike GDAL
> and
> > the NetCDF lib itself) raster uses explicit extension to dispatch to
> > different format logic code, though it ultimately sends it down to rgdal
> to
> > deal with if it can't recognize it - which is why I'm surprised you can't
> > get it to work and ( I'm guessing wildly now):
> >
> > Do you not have rgdal installed?
> >
> >
> > What system are you on? Please use sessionInfo() to share details.
> >
> > Cheers, Mike.
> > On Fri, 1 Dec 2017, 06:14 Antonio Silva, <aolinto@gmail.com> wrote:
> >
> >> Hello
> >>
> >> Some time ago I prepared scripts to extract temperature data from Modis
> >> Aqua files. It can be found at https://gist.github.com/aolinto
> >>
> >> HDF files can be downloaded at https://oceancolor.gsfc.nasa.gov/cgi/l3
> >>
> >> I got the Aquarius sea surface salinity smoothed file from June 2015.
> >>
> >> I could open and read the file:
> >>
> >> library(ncdf4)
> >> library(raster)
> >>
> >> nc.data<-nc_open("Q20151522015181.L3m_MO_SCISM_V5.0_SSS_1deg")
> >> print(nc.data)
> >> dim(ncvar_get(nc.data,"l3m_data"))
> >> ncvar_get(nc.data,"l3m_data")[c(110:160),c(110:117)]
> >>
> >> But I could not prepare a raster from it. I tryed many things as:
> >>
> >> rst.data <-
> >> raster("Q20151522015181.L3m_MO_SCISM_V5.0_SSS_1deg",varname="l3m_data")
> >> Error in .local(.Object, ...) :
> >>  `AQUARIUS/Q20151522015181.L3m_MO_SCISM_V5.0_SSS_1deg' not recognised
> as a
> >> supported file format.
> >>
> >> Error in .rasterObjectFromFile(x, band = band, objecttype =
> "RasterLayer",
> >> :
> >>  Cannot create a RasterLayer object from this file.
> >>
> >> and variations with band and layer.
> >>
> >> I would greatly appreciate a

[R-sig-Geo] make a raster from Aquarius files

2017-11-30 Thread Antonio Silva

Hello

Some time ago I prepared scripts to extract temperature data from Modis
Aqua files. It can be found at https://gist.github.com/aolinto

HDF files can be downloaded at https://oceancolor.gsfc.nasa.gov/cgi/l3

I got the Aquarius sea surface salinity smoothed file from June 2015.

I could open and read the file:

library(ncdf4)
library(raster)

nc.data<-nc_open("Q20151522015181.L3m_MO_SCISM_V5.0_SSS_1deg")
print(nc.data)
dim(ncvar_get(nc.data,"l3m_data"))
ncvar_get(nc.data,"l3m_data")[c(110:160),c(110:117)]

But I could not prepare a raster from it. I tryed many things as:

rst.data <-
raster("Q20151522015181.L3m_MO_SCISM_V5.0_SSS_1deg",varname="l3m_data")
Error in .local(.Object, ...) :
  `AQUARIUS/Q20151522015181.L3m_MO_SCISM_V5.0_SSS_1deg' not recognised as a
supported file format.

Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer",
:
  Cannot create a RasterLayer object from this file.

and variations with band and layer.

I would greatly appreciate any suggestions to solve this issue.

Thanks

-- 
Antônio Olinto Ávila da Silva
Fisheries Institute
São Paulo, Brasil

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] getting mean values for each raster in a stack

2016-05-16 Thread Antonio Silva

Dear Pebesma

No I din't, thanks. Probably I have not read enough.

extract package:raster
fun: function to summarize the values (e.g. ‘mean’). The function
  should take a single numeric vector as argument and return a
  single value (e.g. mean, min or max), and accept a ‘na.rm’
  argument ... 樂

Thanks again.

A.O.


2016-05-16 18:01 GMT-03:00 Edzer Pebesma <edzer.pebe...@uni-muenster.de>:

> Did you try:
>
> extract(stack, polygon, mean)
>
> ?
>
> On 16/05/16 22:50, Antonio Silva wrote:
> > Hi,
> >
> > After many trials and some reading I got what I need
> >
> > apply(vals[[1]],2,mean)
> >
> > Anyway I would welcome comments to improve this solution.
> >
> > Best regards
> >
> > Antonio Olinto
> >
> >
> > 2016-05-16 15:32 GMT-03:00 Antonio Silva <aolinto@gmail.com>:
> >
> >> Dear R users
> >>
> >> I have a raster stack and a polygon.
> >>
> >> I want the mean value inside this polygon for each raster in the stack.
> >>
> >> I'm doing the following:
> >>
> >> vals <- extract(stack,polygon)
> >> colMeans(r.vals[[1]]) # only for mean values
> >> # or
> >> for (i in 1:nlayers(stack)) {print(mean(vals[[1]][,i]))} # mean or other
> >> parameter
> >>
> >> I also tried to use lapply, but I got only the overall mean with
> >> lapply(vals, FUN=mean)
> >>
> >> I wonder if there is a way to get the mean (and other parameters) for
> each
> >> layer using lapply function.
> >>
> >> Thanks in advance.
> >>
> >> Antonio Olinto
> >>
> >>
> >
> >
>
> --
> Edzer Pebesma
> Institute for Geoinformatics  (ifgi),  University of Münster
> Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
> Journal of Statistical Software:   http://www.jstatsoft.org/
> Computers & Geosciences:   http://elsevier.com/locate/cageo/
> Spatial Statistics Society http://www.spatialstatistics.info
>
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>



-- 
Antônio Olinto Ávila da Silva
Biólogo / Oceanógrafo
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] getting mean values for each raster in a stack

2016-05-16 Thread Antonio Silva

Hi,

After many trials and some reading I got what I need

apply(vals[[1]],2,mean)

Anyway I would welcome comments to improve this solution.

Best regards

Antonio Olinto


2016-05-16 15:32 GMT-03:00 Antonio Silva <aolinto@gmail.com>:

> Dear R users
>
> I have a raster stack and a polygon.
>
> I want the mean value inside this polygon for each raster in the stack.
>
> I'm doing the following:
>
> vals <- extract(stack,polygon)
> colMeans(r.vals[[1]]) # only for mean values
> # or
> for (i in 1:nlayers(stack)) {print(mean(vals[[1]][,i]))} # mean or other
> parameter
>
> I also tried to use lapply, but I got only the overall mean with
> lapply(vals, FUN=mean)
>
> I wonder if there is a way to get the mean (and other parameters) for each
> layer using lapply function.
>
> Thanks in advance.
>
> Antonio Olinto
>
>


-- 
Antônio Olinto Ávila da Silva
Biólogo / Oceanógrafo
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil

[[alternative HTML version deleted]]

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

[R-sig-Geo] getting mean values for each raster in a stack

2016-05-16 Thread Antonio Silva

Dear R users

I have a raster stack and a polygon.

I want the mean value inside this polygon for each raster in the stack.

I'm doing the following:

vals <- extract(stack,polygon)
colMeans(r.vals[[1]]) # only for mean values
# or
for (i in 1:nlayers(stack)) {print(mean(vals[[1]][,i]))} # mean or other
parameter

I also tried to use lapply, but I got only the overall mean with
lapply(vals, FUN=mean)

I wonder if there is a way to get the mean (and other parameters) for each
layer using lapply function.

Thanks in advance.

Antonio Olinto

[[alternative HTML version deleted]]

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

[R-sig-Geo] converting shapes to raster

2016-02-04 Thread Antonio Silva

Dear list members



In order to run a niche overlap analysis with nicheOverlap (dismo package)
I need to convert some shapefiles to raster. I dropped a shapefile as
example at https://www.dropbox.com/s/g6o0vfdx7lr34ma/Species1Shape.rar?dl=0

Marine fish abundance data is aggregated in quadrants of 10 x 10 nautical
miles. As I have a coast line some quadrants are cut and also I have a
blank area (land area).

I tried many things like:

shp.frtcost <- readShapePoly("BL10_FRT_COSTEIRA.shp")

plot(shp.frtcost,axes=TRUE,xlim=c(-48,-45),border="grey")

text(coordinates(shp.frtcost),labels=as.character(shp.frtcost@data
$DENS),cex=0.7)

r <- raster(extent(shp.frtcost),nrows=10,ncols=15)

# or raster(extent(shp.frtcost),nrows=10,ncols=15,vals=shp.frtcost@data
$DENS)

r <- rasterize(shp.frtcost,r,fun="first")

plot(shp.frtcost,axes=TRUE,xlim=c(-48,-45),border="grey")

plot(r,add=TRUE)



Shape and raster don't overlay correctly and also I could not pass
densities values to raster.


Limits indicated at extent function do not indicate the limits of my grid.



Thanks for any help.



Best regards


Antônio Olinto Ávila da Silva
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] SST and chlorophyll data extraction in R

2015-09-30 Thread Antonio Silva

Thanks Joseph, they are already there.

https://gist.github.com/aolinto/79e184f6c156c6ab21b3

https://gist.github.com/aolinto/436d9d69798dd95e5ab0

Regards

AO

2015-09-29 17:41 GMT-03:00 Stachelek, Joseph :

> Hi Antonio,
>
> My suggestion would be turn your scripts into Github "gists" so they could
> be more widely accessible.
>
> https://gist.github.com/
>
> Joseph Stachelek
>

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] SST and chlorophyll data extraction in R

2015-09-29 Thread Antonio Silva

Hello,

Some weeks ago I wrote to this list asking about available methods to
extract sst and chlorophyll concentration data from hdf files.

After some work and much reading I wrote two scripts (see below) to do the
job in batch mode, they are very similar. I began to use nc files. I hope
it is useful for other.

Please let me know if there are any error or inconsistency. Suggestions for
improvement are welcome.

Regards

Antonio Olinto


*SST===*

# Antonio Olinto Avila-da-Silva, Instituto de Pesca, Brasil
# script to process Aqua MODIS Sea Surface Temperature
# Monthly means 9 km resolution
# files downloaded from http://oceancolor.gsfc.nasa.gov/cgi/l3
# all .L3m_MO_SST_sst_9km.nc files must be in the working directory
# the script will open each nc file, read lon, lat and sst data,
# select data from a specific area and write them into
# a single csv file named MODISA_sst.csv
# Some reference pages
# http://geog.uoregon.edu/GeogR/topics/netCDF-read-ncdf4.html
# https://scottishsnow.wordpress.com/2014/08/24/many-rastered-beast/

# load libraries
# ncdf4 needs libnetcdf-dev netcdf-bin in Linux
# install.packages(c("ncdf4","reshape2"))
library("ncdf4")
library("reshape2")

# set working directory
setwd("/mnt/Dados02/MODIS/SST")   # indicate the path to the files
file.exists("MODISA_sst.csv") # caution new data will be appended to
this file if it already exists
# file.rename("MODISA_sst.csv","MODISA_sst.old")
# file.remove("MODISA_sst.csv")

# list and remove objects
ls()
rm(list = ls())

# set the study area
lonmax<--40
lonmin<--50
latmax<--22
latmin<--29

# create a list of files and indicate its length
(f <- list.files(".", pattern="*.L3m_MO_SST_sst_9km.nc",full.names=F))
(lf<-length(f))

# variable
var<-"sst"

for (i in 1:lf) {
  # progress indicator
  print(paste("Processing file",i,"from",length(f),sep=" "))
  # open netCDF file
  data<-nc_open(f)
  # extract data
  lon<-ncvar_get(data,"lon")
  lat<-ncvar_get(data,"lat")
  value<-ncvar_get(data,var)
  unit<-ncatt_get(data,var,"units")$value
  # matrix to data.frame
  dimnames(value)<-list(lon=lon,lat=lat)
  dat.var<-melt(value,id="lon")
  # select data from the study area taking out missing data
  dat.varSAtmp<-subset(dat.var,lon<=lonmax & lon>=lonmin & lat<=latmax &
lat>=latmin & value<45)
  # extract date information
  dateini<-ncatt_get(data,0,"time_coverage_start")$value
  dateend<-ncatt_get(data,0,"time_coverage_end")$value

datemean<-mean(c(as.Date(dateend,"%Y-%m-%dT%H:%M:%OSZ"),as.Date(dateini,"%Y-%m-%dT%H:%M:%OSZ")))
  year<-substring(datemean,0,4)
  month<-substring(datemean,6,7)
  # prepare final data set

dat.varSA<-data.frame(rep(as.integer(year,nrow(dat.varSAtmp))),rep(as.integer(month,nrow(dat.varSAtmp))),
  dat.varSAtmp,rep(unit,nrow(dat.varSAtmp)),rep(var,nrow(dat.varSAtmp)))
  names(dat.varSA)<-c("year","month","lon","lat","value","unit","var")
  # save csv file
  fe<-file.exists("MODISA_sst.csv")

write.table(dat.varSA,"MODISA_sst.csv",row.names=FALSE,col.names=!fe,sep=";",dec=",",append=fe)
  # close connection
  nc_close(data)
  # clean workspace

rm(data,lon,lat,value,unit,dat.var,dat.varSAtmp,dateini,dateend,datemean,year,month,dat.varSA,fe)
}
rm(var,f,i,latmax,latmin,lf,lonmax,lonmin)


*CHL_a=*

# Antonio Olinto Avila-da-Silva, Instituto de Pesca, Brasil
# script to process Aqua MODIS Chlorophyll Concentration OCx Algorithm
# Monthly means 9 km resolution
# see http://oceancolor.gsfc.nasa.gov/WIKI/OCChlOCI.html
# files downloaded from http://oceancolor.gsfc.nasa.gov/cgi/l3
# all .L3m_MO_CHL_chlor_a_9km.nc files must be in the working directory
# the script will open each nc file, read lon, lat and chl data,
# select data from a specific area and write them into
# a single csv file named MODISA_chl.csv
# Some reference pages
# http://geog.uoregon.edu/GeogR/topics/netCDF-read-ncdf4.html
# https://scottishsnow.wordpress.com/2014/08/24/many-rastered-beast/

# load libraries
# ncdf4 needs libnetcdf-dev netcdf-bin in Linux
# install.packages(c("ncdf4","reshape2"))
library("ncdf4")
library("reshape2")

# set working directory
setwd("/mnt/Dados02/MODIS/CHL")   # indicate the path to the files
file.exists("MODISA_chl.csv") # caution new data will be appended to
this file if it already exists
# file.rename("MODISA_chl.csv","MODISA_chl.old")
# file.remove("MODISA_chl.csv")

# list and remove objects
ls()
rm(list = ls())

# set the study area
lonmax<--40
lonmin<--50
latmax<--22
latmin<--29

# create a list of files and indicate its length
(f <- list.files(".", pattern="*.L3m_MO_CHL_chlor_a_9km.nc",full.names=F))
(lf<-length(f))

# variable
var<-"chlor_a"

for (i in 1:lf) {
  # progress indicator
  print(paste("Processing file",i,"from",length(f),sep=" "))
  # open netCDF file
  data<-nc_open(f)
  # extract data
  lon<-ncvar_get(data,"lon")
  lat<-ncvar_get(data,"lat")
  value<-ncvar_get(data,var)
  unit<-ncatt_get(data,var,"units")$value
  # matrix to data.frame
  dimnames(value)<-list(lon=lon,lat=lat)

[R-sig-Geo] gdal_translate and gdalwarp question

2015-09-20 Thread Antonio Silva

Hello,

I have some doubts on the usage of some gdal tools.

After converting a hdf file to tif I want to reproject to SIRGAS2000 and
clip between lats of 22 to 29 S and lats of 40 to 50 W.

HDF file can be downloaded at
https://app.box.com/s/16cf7qv6af6gsz1v66staori2mtneu0r

Basically I'm following
https://scottishsnow.wordpress.com/2014/08/24/many-rastered-beast/

Well to convert HDF file I'm using:

gdal_translate("A20080012008031.L3m_MO_SST_4","georef.tif",sd_index=1,a_ullr=c(0,4320,8640,0),
a_srs="+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=6371007
+b=6371007 +units=m +no_defs")

Without a_ullr and a_srs options I was getting an error message when using
gdalwarp: "ERROR 1: Unable to compute a transformation between pixel/line
and georeferenced coordinates"

a_ullr and a_srs values I got with GDALinfo("georef.tif")

map <- raster("georef.tif")
plot(map)

My problem now is reproject to SIRGAS2000 and clip the image georef.tif:

gdalwarp("georef.tif", "georef2.tif",
s_srs="+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=6371007
+b=6371007 +units=m +no_defs",
t_srs="+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs")

map2 <- raster("georef2.tif")
plot(map2)

Map2 is not in SIRGAS2000 projection and clipping option
te=c(-50,-29,-40,-22) does not work.

Where is my mistake? I hope someone can tell me.

Thanks for any help.

Antonio Olinto



















-- 
Antônio Olinto Ávila da Silva
Biólogo / Oceanógrafo
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil

[[alternative HTML version deleted]]

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

[R-sig-Geo] SST and chlorophyll data extraction in R

2015-09-18 Thread Antonio Silva

Dear R-Sig-Geo fellows


This is my first post to this list. My name is Antonio Olinto and I work as
researcher at Fisheries Institute, São Paulo, Brazil.


Well, I'd like to hear suggestions on how to proceed in R to extract sea
surface temperature and chlorophyll concentration data and calculate some
statistics from MODIS/AQUA data in a given area.


I used to do this in QGIS but available statistical parameters are limited
and I had to work with one image at a time. I found a reference to a
package named satin (
https://www.r-project.org/conferences/useR-2009/slides/Villalobos+GonzalezRodriguez.pdf)
but it seems it's not available anymore.


I have several  *.L3m_MO_SST_4.bz2 and *.L3m_MO_CHL_chlor_a_4km.bz2 files
obtained from
http://podaac.jpl.nasa.gov/dataset/MODIS_AQUA_L3_SST_THERMAL_MONTHLY_4KM_DAYTIME
and from http://oceandata.sci.gsfc.nasa.gov/MODISA/Mapped/Monthly/4km/chlor/


I also have a shapefile with the boundaries of the area of interest. I'd
like to have some statistics (maximum, minimum, mean, quantiles, ect.) of
both variables in this area.


Any indication of tutorials, sites or packages will be welcome.  Thanks in
advance.


All the best,

Antônio Olinto

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] extract data from multilayer nc file

Re: [R-sig-Geo] extract data from multilayer nc file

[R-sig-Geo] extract data from multilayer nc file

[R-sig-Geo] select points out of a polygon

Re: [R-sig-Geo] how to get an attribute from a raster layer

[R-sig-Geo] how to get an attribute from a raster layer

Re: [R-sig-Geo] get data from nc file

[R-sig-Geo] get data from nc file

Re: [R-sig-Geo] subsetting a spatial polygons

[R-sig-Geo] subsetting a spatial polygons

Re: [R-sig-Geo] drawing a polygon from several lines

Re: [R-sig-Geo] drawing a polygon from several lines

[R-sig-Geo] drawing a polygon from several lines

Re: [R-sig-Geo] how to create several polygons from a list of vertices

Re: [R-sig-Geo] how to create several polygons from a list of vertices

[R-sig-Geo] how to create several polygons from a list of vertices

Re: [R-sig-Geo] make a raster from Aquarius files

[R-sig-Geo] make a raster from Aquarius files

Re: [R-sig-Geo] getting mean values for each raster in a stack

Re: [R-sig-Geo] getting mean values for each raster in a stack

[R-sig-Geo] getting mean values for each raster in a stack

[R-sig-Geo] converting shapes to raster

Re: [R-sig-Geo] SST and chlorophyll data extraction in R

Re: [R-sig-Geo] SST and chlorophyll data extraction in R

[R-sig-Geo] gdal_translate and gdalwarp question

[R-sig-Geo] SST and chlorophyll data extraction in R

26 matches

Site Navigation

Mail list logo

Footer information