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

Re: sp::disaggregate() does not seem to add the necessary "comment" attribute to identify holes if Polygons does not contain exactly 1 polygon with 0 or 1 holes

Sat, 04/29/2017 - 09:06
On Sat, 29 Apr 2017, Anne C. Hanna wrote:

> Roger,
>
> This looks great, and I will try it out ASAP.  I do have one reservation
> though --- it seems you are using createSPComment() to reconstruct the
> comments, and I have seen some discussion that that may not be reliable
> in all cases (e.g. if the initial polygons are wonky in some way).  I
> don't know a lot about this, but is it possible that it would be
> preferable to parse and appropriately disaggregate the original comments
> strings (if they exist), so as to deal slightly more smoothly with such
> cases (i.e., the polygons would still be wonky, but at least the
> hole/polygon matching would track with whatever was in the original
> data)?
Contributions very welcome - this was code moved from the raster package,
so I really have little feeling for why it is necessary.

Please provide references and use cases. GEOS is strict in its treatments
of geometries, but does work in most cases. People just expressing
opinions doesn't help, and may well be misleading. It is possible to clean
geometries too, see the cleangeo package. createPolygonsComment probably
isn't foolproof, but specific reports help, because then it is possible to
do something about it. sp classes didn't use simple features when written
because sf was not widely used then - introduced in 2004, when sp was
being developed. Using sf helps, and rgeos::createPolygonsComment was a
work-around from 2010 when rgeos was written.

Even when you use sf in the sf package, you will still run into invalid
geometries because the geometries are in fact invalid, the sf package
also uses GEOS.

Roger

>
> I also had no success using createSPComment to fix disaggregate()'s
> output previously, even though my polygons are perfectly non-wonky, so I
> am perhaps a little more untrusting of it than I should be.  But I'll
> let you know how this version works with my data.  Thanks for addressing
> this so quickly!
>
> - Anne
>
>
> On 04/28/2017 12:11 PM, Roger Bivand wrote:
>> I've pushed the fix to my fork:
>>
>> https://github.com/rsbivand/sp
>>
>> and created a pull request:
>>
>> https://github.com/edzer/sp/pull/28
>>
>> Only one part of a complicated set if nested if() in disaggregate() was adding
>> comments, but in some settings the existing comments survived the
>> disaggregation. Now the Polygons object comment attribute is re-created for
>> all Polygons objects. This is version 1.2-6, also including code changes that
>> internally affect rgdal and rgeos - you may like to re-install them from
>> source after installing this sp version (shouldn't matter).
>>
>> Roger
>>
>> On Fri, 28 Apr 2017, Anne C. Hanna wrote:
>>
>>> Hello.  I first posted this issue report on the sp GitHub repo
>>> (https://github.com/edzer/sp/issues/27) and it was suggested that I redirect
>>> it to here.
>>>
>>> I am working with a geographic dataset with complex borders. The data is
>>> stored as a SpatialPolygonsDataFrame.  Each Polygons object in the data frame
>>> may be composed of multiple disjoint polygons, and each polygon may have
>>> multiple holes.  I want to disaggregate each of the Polygons objects into its
>>> individual disjoint polygons and construct an adjacency matrix for all the
>>> disjoint components, and I was using disaggregate() to do this.  However, when
>>> I run gTouches() on the disaggregated data, in order to compute the
>>> adjacencies, I get a number of warnings like this:
>>>
>>> Error in RGEOSBinPredFunc(spgeom1, spgeom2, byid, "rgeos_touches") : Polygons
>>> object missing comment attribute ignoring hole(s).  See function
>>> createSPComment.
>>>
>>> Looking at the Polygons "comment" attributes in the SpatialPolygonsDataFrame
>>> output by disaggregate(), I see that the only comment values are "0"
>>> (indicating a single polygon with no holes), "0 1" (indicating a single
>>> polygon with a single hole), and NULL (apparently no comment was written).
>>> Since I know my dataset contains several Polygons objects which are composed
>>> of multiple disjoint regions, and also several Polygons which contain more
>>> than one hole, this is not the expected result.  In reading the disaggregate()
>>> code in the sp GitHub repository (specifically, explodePolygons()), I also
>>> can't see anywhere the comment is being added for the cases where a Polygons
>>> object has more than two parts or more than two holes.  It actually seems like
>>> it's getting carried along almost accidentally in the few cases that do get
>>> comments, and neglected otherwise.
>>>
>>> Assuming I'm not failing to understand the code and the desired behavior
>>> (entirely possible, as I am new at working with this software!), this seems
>>> suboptimal to me.  My dataset is pretty well-behaved (despite its complexity),
>>> so I should be able to fix my issues with judicious application of
>>> createPolygonsComment.  But I had a heck of a time figuring out what was going
>>> wrong with gTouches, since Polygons comment management appears to be a pretty
>>> obscure field (and createSPComment wasn't working for me, for whatever
>>> reason).  So it seems like it might be better if disaggregate() just parses
>>> and passes along the comments from its input correctly, or, if it's absolutely
>>> necessary to not create comments, passes nothing and warns clearly in the
>>> manual that comments and associated hole information are being lost.  Passing
>>> along comments in some cases while silently dropping them in others seems like
>>> kind of the worst of both worlds.
>>>
>>> I've attached a set of tests I wrote to demonstrate desired/undesired
>>> behavior: disaggregate_comment_tests.R.  My R version is 3.4.0, my sp version
>>> is 1.2-4, my rgeos version is 0.3-23 (SVN revision 546), and my GEOS runtime
>>> version is 3.5.1-CAPI-1.9.1 r4246.  I am using Debian Release 9.0 with kernel
>>> version 4.9.0-2-amd64.  I hope this is useful; please let me know if you need
>>> more info or if there is a better place to post this issue.
>>>
>>> - Anne
>>>
>>
>
>
--
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

Using gdalUtils to reproject hdf file

Fri, 04/28/2017 - 19:22
 Dear all,

I am now working with GEFED dataset, which is a global wild file dataset.You can download one from FTP: fuoco.geog.umd.edu; Login name: fire; Password: burnt; The file is stored into gfed4/monthly/. Any one file will be example.

My task is to translate the hdf file into tif, which can be done by gdalUtils package. But there is no georeferecing information in the file. Instead, the only information I could find is from the data manual (can be found in gfed4 folder in that FTP): "Each product file contains seven data layers, each stored as a separate HDF4 Scientific Data Set (SDS). Units apply to each layer after multiplying by the specified scale factor. Each data layer has 720 rows and 1440 columns which correspond to the global 0.25◦ GFED grid. The center of the upper left grid cell is located at longitude 179.875◦W, 89.875◦N.". I think the original src can be written from this statement. And I can pass it to a_src and a_ullr variable in gdal_translate command. But I have no idea what is the string for a_src.

Any recommendation will be appreciated.

Tianyi



        [[alternative HTML version deleted]]

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

Re: sp::disaggregate() does not seem to add the necessary "comment" attribute to identify holes if Polygons does not contain exactly 1 polygon with 0 or 1 holes

Fri, 04/28/2017 - 18:25
Roger,

This looks great, and I will try it out ASAP.  I do have one reservation
though --- it seems you are using createSPComment() to reconstruct the
comments, and I have seen some discussion that that may not be reliable in all
cases (e.g. if the initial polygons are wonky in some way).  I don't know a
lot about this, but is it possible that it would be preferable to parse and
appropriately disaggregate the original comments strings (if they exist), so
as to deal slightly more smoothly with such cases (i.e., the polygons would
still be wonky, but at least the hole/polygon matching would track with
whatever was in the original data)?

I also had no success using createSPComment to fix disaggregate()'s output
previously, even though my polygons are perfectly non-wonky, so I am perhaps a
little more untrusting of it than I should be.  But I'll let you know how this
version works with my data.  Thanks for addressing this so quickly!

 - Anne


On 04/28/2017 12:11 PM, Roger Bivand wrote:
> I've pushed the fix to my fork:
>
> https://github.com/rsbivand/sp
>
> and created a pull request:
>
> https://github.com/edzer/sp/pull/28
>
> Only one part of a complicated set if nested if() in disaggregate() was adding
> comments, but in some settings the existing comments survived the
> disaggregation. Now the Polygons object comment attribute is re-created for
> all Polygons objects. This is version 1.2-6, also including code changes that
> internally affect rgdal and rgeos - you may like to re-install them from
> source after installing this sp version (shouldn't matter).
>
> Roger
>
> On Fri, 28 Apr 2017, Anne C. Hanna wrote:
>
>> Hello.  I first posted this issue report on the sp GitHub repo
>> (https://github.com/edzer/sp/issues/27) and it was suggested that I redirect
>> it to here.
>>
>> I am working with a geographic dataset with complex borders. The data is
>> stored as a SpatialPolygonsDataFrame.  Each Polygons object in the data frame
>> may be composed of multiple disjoint polygons, and each polygon may have
>> multiple holes.  I want to disaggregate each of the Polygons objects into its
>> individual disjoint polygons and construct an adjacency matrix for all the
>> disjoint components, and I was using disaggregate() to do this.  However, when
>> I run gTouches() on the disaggregated data, in order to compute the
>> adjacencies, I get a number of warnings like this:
>>
>> Error in RGEOSBinPredFunc(spgeom1, spgeom2, byid, "rgeos_touches") : Polygons
>> object missing comment attribute ignoring hole(s).  See function
>> createSPComment.
>>
>> Looking at the Polygons "comment" attributes in the SpatialPolygonsDataFrame
>> output by disaggregate(), I see that the only comment values are "0"
>> (indicating a single polygon with no holes), "0 1" (indicating a single
>> polygon with a single hole), and NULL (apparently no comment was written).
>> Since I know my dataset contains several Polygons objects which are composed
>> of multiple disjoint regions, and also several Polygons which contain more
>> than one hole, this is not the expected result.  In reading the disaggregate()
>> code in the sp GitHub repository (specifically, explodePolygons()), I also
>> can't see anywhere the comment is being added for the cases where a Polygons
>> object has more than two parts or more than two holes.  It actually seems like
>> it's getting carried along almost accidentally in the few cases that do get
>> comments, and neglected otherwise.
>>
>> Assuming I'm not failing to understand the code and the desired behavior
>> (entirely possible, as I am new at working with this software!), this seems
>> suboptimal to me.  My dataset is pretty well-behaved (despite its complexity),
>> so I should be able to fix my issues with judicious application of
>> createPolygonsComment.  But I had a heck of a time figuring out what was going
>> wrong with gTouches, since Polygons comment management appears to be a pretty
>> obscure field (and createSPComment wasn't working for me, for whatever
>> reason).  So it seems like it might be better if disaggregate() just parses
>> and passes along the comments from its input correctly, or, if it's absolutely
>> necessary to not create comments, passes nothing and warns clearly in the
>> manual that comments and associated hole information are being lost.  Passing
>> along comments in some cases while silently dropping them in others seems like
>> kind of the worst of both worlds.
>>
>> I've attached a set of tests I wrote to demonstrate desired/undesired
>> behavior: disaggregate_comment_tests.R.  My R version is 3.4.0, my sp version
>> is 1.2-4, my rgeos version is 0.3-23 (SVN revision 546), and my GEOS runtime
>> version is 3.5.1-CAPI-1.9.1 r4246.  I am using Debian Release 9.0 with kernel
>> version 4.9.0-2-amd64.  I hope this is useful; please let me know if you need
>> more info or if there is a better place to post this issue.
>>
>> - Anne
>>
>
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
signature.asc (201 bytes) Download Attachment

Re: sp::disaggregate() does not seem to add the necessary "comment" attribute to identify holes if Polygons does not contain exactly 1 polygon with 0 or 1 holes

Fri, 04/28/2017 - 11:11
I've pushed the fix to my fork:

https://github.com/rsbivand/sp

and created a pull request:

https://github.com/edzer/sp/pull/28

Only one part of a complicated set if nested if() in disaggregate() was
adding comments, but in some settings the existing comments survived the
disaggregation. Now the Polygons object comment attribute is re-created
for all Polygons objects. This is version 1.2-6, also including code
changes that internally affect rgdal and rgeos - you may like to
re-install them from source after installing this sp version (shouldn't
matter).

Roger

On Fri, 28 Apr 2017, Anne C. Hanna wrote:

> Hello.  I first posted this issue report on the sp GitHub repo
> (https://github.com/edzer/sp/issues/27) and it was suggested that I redirect
> it to here.
>
> I am working with a geographic dataset with complex borders. The data is
> stored as a SpatialPolygonsDataFrame.  Each Polygons object in the data frame
> may be composed of multiple disjoint polygons, and each polygon may have
> multiple holes.  I want to disaggregate each of the Polygons objects into its
> individual disjoint polygons and construct an adjacency matrix for all the
> disjoint components, and I was using disaggregate() to do this.  However, when
> I run gTouches() on the disaggregated data, in order to compute the
> adjacencies, I get a number of warnings like this:
>
> Error in RGEOSBinPredFunc(spgeom1, spgeom2, byid, "rgeos_touches") : Polygons
> object missing comment attribute ignoring hole(s).  See function createSPComment.
>
> Looking at the Polygons "comment" attributes in the SpatialPolygonsDataFrame
> output by disaggregate(), I see that the only comment values are "0"
> (indicating a single polygon with no holes), "0 1" (indicating a single
> polygon with a single hole), and NULL (apparently no comment was written).
> Since I know my dataset contains several Polygons objects which are composed
> of multiple disjoint regions, and also several Polygons which contain more
> than one hole, this is not the expected result.  In reading the disaggregate()
> code in the sp GitHub repository (specifically, explodePolygons()), I also
> can't see anywhere the comment is being added for the cases where a Polygons
> object has more than two parts or more than two holes.  It actually seems like
> it's getting carried along almost accidentally in the few cases that do get
> comments, and neglected otherwise.
>
> Assuming I'm not failing to understand the code and the desired behavior
> (entirely possible, as I am new at working with this software!), this seems
> suboptimal to me.  My dataset is pretty well-behaved (despite its complexity),
> so I should be able to fix my issues with judicious application of
> createPolygonsComment.  But I had a heck of a time figuring out what was going
> wrong with gTouches, since Polygons comment management appears to be a pretty
> obscure field (and createSPComment wasn't working for me, for whatever
> reason).  So it seems like it might be better if disaggregate() just parses
> and passes along the comments from its input correctly, or, if it's absolutely
> necessary to not create comments, passes nothing and warns clearly in the
> manual that comments and associated hole information are being lost.  Passing
> along comments in some cases while silently dropping them in others seems like
> kind of the worst of both worlds.
>
> I've attached a set of tests I wrote to demonstrate desired/undesired
> behavior: disaggregate_comment_tests.R.  My R version is 3.4.0, my sp version
> is 1.2-4, my rgeos version is 0.3-23 (SVN revision 546), and my GEOS runtime
> version is 3.5.1-CAPI-1.9.1 r4246.  I am using Debian Release 9.0 with kernel
> version 4.9.0-2-amd64.  I hope this is useful; please let me know if you need
> more info or if there is a better place to post this issue.
>
> - Anne
>
--
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

patial distribution of the ecological niche of cowpea,

Fri, 04/28/2017 - 09:17
I want to model the spatial distribution of the ecological niche of cowpea, I have a set of data (environmental variables and bioclimatic) and data of presence of the varieties of cowpea.I wanted to use the software Maxent for do it; But I understood through the literature that the R software is now widely used. It is with this in mind that I would like to focus on this. By the way, I need help to understand the procedures for doing this modeling using R. 1. I would like to do a Shapiro-Wilks normality test, the null hypothesis of which is to assert that The data follow a normal distribution and where the alternative hypothesis implies an asymmetric or abnormal distribution, makes it possible to conclude that the observations are of very asymmetric distribution. 2. I would like to do a logarithmic transformation of the response data is therefore performed in order to reduce the asymmetry of the data and thus improve the power of the statistical tests.3. I would like to do a canonical redundancy analysis on all the data tables. 4. I would like to make a bidirectional progressive selection in order to select the relevant variables according to Akaike's information criterion or AIC and Of the representativeness (p-value) .5. I would like to make a progressive selection by adding, according to the adjusted R-square, in order to compare the result of the selection 6. Finally, the procedures to be followed to model the distribution of the varieties.

Best regards 
 SADDA Abou-Soufianou -------------------------------------- DoctorantUniversité Dan Dicko Dankoulodo de Maradi,120 avenue Maman Koraou,ADS, MaradiBP 465 Maradi/Niger E-mail: [hidden email] : (+227)96269987
        [[alternative HTML version deleted]]

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

Re: modélisation de la distribution spatiale du niébé

Fri, 04/28/2017 - 09:12
Cher Sadda,

vous avez écrit que "J'ai voulu utilisé le logiciel Maxent pour le
faire; mais j'ai compris à travers la littérature que le logiciel R est
de nos jours largement utilisé.", mais vous pouvez utiliser MaxEnt via
le logiciel R:
install.packages("dismo")
library(dismo)
m <- maxent(x = tableau_des_prédicteurs_environnementaux, p =
points_de_présence)

Vous devez d'abord télécharger le programme à partir d'ici:
http://biodiversityinformatics.amnh.org/open_source/maxent/

J'espère que cela t'aides,
Ákos Bede-Fazekas
Académie hongroise des sciences


2017.04.28. 15:38 keltezéssel, Soufianou Abou via R-sig-Geo írta:
> Bonsoir chers, j'aimerai modéliser la distribution spatiale de la niche écologique du niébé , je dispose d'un jeu de données (variables environnementales et bioclimatiques) et des données de présence des variétés du niébé.J'ai voulu utilisé le logiciel Maxent pour le faire; mais j'ai compris à travers la littérature que le logiciel R est de nos jours largement utilisé. C'est dans cette optique que j'aimerai m'y plancher là dessus. Au fait, j'ai besoin de l'aide pour comprendre les procédures pour faire cette modélisation à l'aide de R.   1. J'aimerai faire Un test de normalité de Shapiro-Wilks, dont l’hypothèse nulle consiste à affirmer que les données suivent une distribution normale et où l'hypothèse alternative implique une distribution asymétrique ou anormale, permet de conclure que les observations sont de distribution très asymétrique. 2. J'aimerai faire une transformation logarithmique des données réponses est donc effectuée afin de réduire l’asymétrie des données et, ainsi, d’améliorer la puissance des tests statistiques.
> 3. J'aimerai faire une analyse de redondance canonique  sur l’ensemble des tableaux de données .4.J'aimerai faire  une sélection progressive bidirectionnelle  afin  de choir les variables pertinentes en fonction du critère d’information d’Akaike, ou AIC et de la représentativité (p-value) .5. J'aimerai faire  une sélection progressive par ajout, en fonction du R-carré ajusté,  afin de comparer le résultat de la sélection
> 6. En fin les procédure à suivre pour modéliser la distribution des variétés.
> merci cordialement.
> SADDA Abou-Soufianou -------------------------------------- DoctorantUniversité Dan Dicko Dankoulodo de Maradi,120 avenue Maman Koraou,ADS, MaradiBP 465 Maradi/Niger E-mail: [hidden email] : (+227)96269987
> [[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

modélisation de la distribution spatiale du niébé

Fri, 04/28/2017 - 08:38
Bonsoir chers, j'aimerai modéliser la distribution spatiale de la niche écologique du niébé , je dispose d'un jeu de données (variables environnementales et bioclimatiques) et des données de présence des variétés du niébé.J'ai voulu utilisé le logiciel Maxent pour le faire; mais j'ai compris à travers la littérature que le logiciel R est de nos jours largement utilisé. C'est dans cette optique que j'aimerai m'y plancher là dessus. Au fait, j'ai besoin de l'aide pour comprendre les procédures pour faire cette modélisation à l'aide de R.   1. J'aimerai faire Un test de normalité de Shapiro-Wilks, dont l’hypothèse nulle consiste à affirmer que les données suivent une distribution normale et où l'hypothèse alternative implique une distribution asymétrique ou anormale, permet de conclure que les observations sont de distribution très asymétrique. 2. J'aimerai faire une transformation logarithmique des données réponses est donc effectuée afin de réduire l’asymétrie des données et, ainsi, d’améliorer la puissance des tests statistiques.
3. J'aimerai faire une analyse de redondance canonique  sur l’ensemble des tableaux de données .4.J'aimerai faire  une sélection progressive bidirectionnelle  afin  de choir les variables pertinentes en fonction du critère d’information d’Akaike, ou AIC et de la représentativité (p-value) .5. J'aimerai faire  une sélection progressive par ajout, en fonction du R-carré ajusté,  afin de comparer le résultat de la sélection
6. En fin les procédure à suivre pour modéliser la distribution des variétés.
merci cordialement.
SADDA Abou-Soufianou -------------------------------------- DoctorantUniversité Dan Dicko Dankoulodo de Maradi,120 avenue Maman Koraou,ADS, MaradiBP 465 Maradi/Niger E-mail: [hidden email] : (+227)96269987
        [[alternative HTML version deleted]]

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

Extracting a single grid cell from a SpatialGrid keeping the grid structure

Fri, 04/28/2017 - 03:09
Hello,
I would like to extract a single grid cell from a SpatialGrid or a SpatialGridDataFrame object without destroying the grid structure, that is, the resulting object must be 
a SpatialGrid or a SpatialGridDataFrame consisting of a single grid cell. Here is a minimum working example:

library(sp);
data(meuse.grid);
coordinates(meuse.grid) = c("x", "y") # promote to SpatialPointsDataFrame
gridded(meuse.grid) <- TRUE # promote to SpatialPixelsDataFrame
x <- as(meuse.grid, "SpatialGridDataFrame") # creates the full grid
fullgrid(meuse.grid) <- TRUE
image(meuse.grid);

# The following works:
latInd <- 20:70; lonInd <- 10:70; # the grid cell indices along latitude and longitude, respectively
a_df <- meuse.grid[latInd, lonInd, "dist"]; # the example SpatialGridDataFrame object
 a_grid <- as(meuse.grid, "SpatialGrid"); # the example SpatialGrid object
image(a_df, add = TRUE, col = bpy.colors());
# the following prints a warning and assumes the grid cellsize of the single-cell dimension is equal to the cellsize of the 
# other dimension, which is not always appropriate:
lonInd <- 50;
a1_df <- meuse.grid[latInd, lonInd, "dist"]; image(a1_df);
# the following returns a SpatialPointsDataFrame object, which is not what is desired:
latInd <- 50;
a11_df <- meuse.grid[latInd, lonInd, "dist"];
# The following simply fails:
a11_grid <- a_grid[latInd, lonInd];

# a way round is to go through a RasterLayer object:
library(raster);
latInd <- a_grid@[hidden email][2] - latInd + 1; a11_grid <- raster(x=a_grid);
# this is for the SpatialGrid case:
a11_grid <- crop(x=a11_grid, y=extent(x=a11_grid, r1=tail(x=latInd, n=1), r2=latInd[1], c1=lonInd[1], c2=tail(x=lonInd, n=1)));
a11_grid <- as(object=a11_grid, Class="SpatialGrid");
# and this is for the SpatialGridDataFrame case:
a11_df <- raster(x=meuse.grid["dist"]);
a11_df <- crop(x=a11_df, y=extent(x=a11_df, r1=tail(x=latInd, n=1), r2=latInd[1], c1=lonInd[1], c2=tail(x=lonInd, n=1)));
a11_df <- as(object=a11_df, Class="SpatialGridDataFrame");
# This approach also works if I only want to extract a single spatial column or a single spatial row and still keep the original grid structure

I would like to ask if there is a simpler and more elegant way to avoid the RasterLayer object? Certainly, this can be achieved directly by hacking the code of the 
"[" operator for the SpatialGrid and the SpatialGridDataFrame classes. However, I do not know how to access that code ... 

Any suggestions will be appreciated.

Thank you very much for your attention and your time.

Best regards,
Martin

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

sp::disaggregate() does not seem to add the necessary "comment" attribute to identify holes if Polygons does not contain exactly 1 polygon with 0 or 1 holes

Fri, 04/28/2017 - 03:09
Hello.  I first posted this issue report on the sp GitHub repo
(https://github.com/edzer/sp/issues/27) and it was suggested that I redirect
it to here.

I am working with a geographic dataset with complex borders. The data is
stored as a SpatialPolygonsDataFrame.  Each Polygons object in the data frame
may be composed of multiple disjoint polygons, and each polygon may have
multiple holes.  I want to disaggregate each of the Polygons objects into its
individual disjoint polygons and construct an adjacency matrix for all the
disjoint components, and I was using disaggregate() to do this.  However, when
I run gTouches() on the disaggregated data, in order to compute the
adjacencies, I get a number of warnings like this:

Error in RGEOSBinPredFunc(spgeom1, spgeom2, byid, "rgeos_touches") : Polygons
object missing comment attribute ignoring hole(s).  See function createSPComment.

Looking at the Polygons "comment" attributes in the SpatialPolygonsDataFrame
output by disaggregate(), I see that the only comment values are "0"
(indicating a single polygon with no holes), "0 1" (indicating a single
polygon with a single hole), and NULL (apparently no comment was written).
Since I know my dataset contains several Polygons objects which are composed
of multiple disjoint regions, and also several Polygons which contain more
than one hole, this is not the expected result.  In reading the disaggregate()
code in the sp GitHub repository (specifically, explodePolygons()), I also
can't see anywhere the comment is being added for the cases where a Polygons
object has more than two parts or more than two holes.  It actually seems like
it's getting carried along almost accidentally in the few cases that do get
comments, and neglected otherwise.

Assuming I'm not failing to understand the code and the desired behavior
(entirely possible, as I am new at working with this software!), this seems
suboptimal to me.  My dataset is pretty well-behaved (despite its complexity),
so I should be able to fix my issues with judicious application of
createPolygonsComment.  But I had a heck of a time figuring out what was going
wrong with gTouches, since Polygons comment management appears to be a pretty
obscure field (and createSPComment wasn't working for me, for whatever
reason).  So it seems like it might be better if disaggregate() just parses
and passes along the comments from its input correctly, or, if it's absolutely
necessary to not create comments, passes nothing and warns clearly in the
manual that comments and associated hole information are being lost.  Passing
along comments in some cases while silently dropping them in others seems like
kind of the worst of both worlds.

I've attached a set of tests I wrote to demonstrate desired/undesired
behavior: disaggregate_comment_tests.R.  My R version is 3.4.0, my sp version
is 1.2-4, my rgeos version is 0.3-23 (SVN revision 546), and my GEOS runtime
version is 3.5.1-CAPI-1.9.1 r4246.  I am using Debian Release 9.0 with kernel
version 4.9.0-2-amd64.  I hope this is useful; please let me know if you need
more info or if there is a better place to post this issue.

 - Anne

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
disaggregate_comment_tests.R (5K) Download Attachment
signature.asc (201 bytes) Download Attachment

alpha hulls (always complete) and incomplete envelops (confidence intervals)

Mon, 04/24/2017 - 19:19
I have a series of noodles, which is to say line shapes, that I am trying
to come to grips.

If I generalize these lines away to point clouds in time order, then an
alpha hull will contain all points, as would convex hull, but potentially
with less wasted space. The thing about a convex hull or alpha hull is that
all samples are 'within'.  This is to say that it is converged as sample =
infinity, or, more normally, sample = total number of samples available
which is generally less than infinity.

In a very general way, I am wondering if there is an accepted method to
allow, say 30 percent of the 'alpha hulled' samples (which clearly are not
designed to allow such) to reside outside the alpha hull (essentially
creating confidence intervals upon the alpha hull (or perhaps I haven't
read enough)) . Secondarily, is there a way to compare 'confidence
interval' alpha hulls, where 70 percent of the sample points reside within,
and the rest exogenous.

I ask as my noodles may share, to an unknown amount, but perhaps
extensively, commonalities of co-existance to as much as 70 percent (it is
speculated) of an alpha hull space, and the variance that I am trying to
account for is the 30 percent that I would guess and perhaps to define as a
separate alpha hull, or some sort space, outside the alpha hull space.

Having it both ways: *if *100% of points are within the alpha hull, how
might one reduce this to 70% or some such, because hulls are always sample
complete.

I realize that this probably sounds inchoate, or am I fantabulizing, but I
think I am asking about comparing the shapes of constrained (incomplete)
alpha hulls in the context of Parzen windows (whose shape I wonder about,
boxes?).

Another way of looking at this a point process is that given a 1280x1280
grid, there are an enormous number of cells that will always be NA, a
smaller number that will be visited once. So how to proceed to compare the
noodles.

While I sense I will be killed on this question: Any thoughts or suggested
reading appreciated so I might address more intelligently anon.
Chris

        [[alternative HTML version deleted]]

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

Error in GDAL for library MODIS

Mon, 04/24/2017 - 18:02
I keep getting a GDAL error when I try to use the MODIS library to download
and process MODIS data.

PROCESSING:
_______________
GDAL           : Not available. Use 'MODIS:::checkTools('GDAL')' for more
information!

Note: I have rgdal installed, along with GDAL separately installed from
osgeo4w.

Any suggestions?


--
Data Manager,
Barbara Han Lab,
Cary Institute of Ecosystem Studies,
2801 Sharon Turnpike, Millbrook, NY 12545
Lab Site : http://www.hanlab.science/
Personal Site : http://evolecol.weebly.com/
Phone : (845)-677-7600 Ext: 241

        [[alternative HTML version deleted]]

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

Re: Pycnophylactic Interpolation

Mon, 04/24/2017 - 09:44
Breno,

Just an article comparing pycnophylatic interpolation to area-to-point
kriging
https://www.researchgate.net/publication/227727493_Reconstructing_Population_Density_Surfaces_from_Areal_Data_A_Comparison_of_Tobler's_Pycnophylactic_Interpolation_Method_and_Area-to-Point_Kriging
Chris

On Sun, Apr 23, 2017 at 4:35 PM, Breno Oliveira <[hidden email]>
wrote:

> Hi,
>
> I'm a student of Information Systems in Brazil and a junior software
> archtect. To my final work on the undergraduate program I'm implementing a
> software with the Potential model. After calculate the potential values to
> each point of analysis I would like to be able to create a map to represent
> the result. The map should looks like this one:
> https://www.r-bloggers.com/contour-and-density-layers-with-ggmap/. I'm
> new in R so I have no a big knowhow and a lot of doubts. I would like ask
> for community's help. My problem is that I have no to many points like the
> example so the map looks incorrect, actualy I have a few points (x,y and
> z). I think I must Interpolate these points to get a good result, I would
> like to use the Pycnophylactic Interpolation avaliable in pycno package.
>
> Here is a example of data I obtain (no interpolation applied):
>
>
> City    Area    Score   latitudes       longitudes      Potentials
> Belo Horizone   330,95  2513451 -19,91663263    -43,93337198    164874,1295
> Betim   346     412003  -19,96730327    -44,20119724    82629,32514
> Contagem        195,268 648766  -19,91615326    -44,08087722    140851,9626
> Sarzedo 61,892  25798   -20,06250329    -44,11504723    79307,30623
> Rio Acima       50      90      -20,08830329    -43,79104717    55226,91779
>
>
> Please let me know if you need any extra information.
> Thank you for your patience and cooperation.
> Breno Oliveira
>
>
>         [[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: R-help - Shiny and Datatable

Mon, 04/24/2017 - 09:01
On Mon, Apr 24, 2017 at 6:00 AM, <[hidden email]> wrote:
>
> Message: 4
> Date: Sun, 23 Apr 2017 18:19:59 -0400
> From: Joe Larson <[hidden email]>
>
> Hello,
>  I am getting wrapped around the axial on a shiny project, being a new R
> learner I am asking for some help.... If you want the full code, to let me
> know.


Yes, if you can post your code and data somewhere that would make it easier
to help.

Kent

        [[alternative HTML version deleted]]

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

Re: R-help - Shiny and Datatable

Sun, 04/23/2017 - 19:16
Joe,

If you look again at the CDC, USA 500 ,(
https://chronicdata.cdc.gov/500-Cities/500-Cities-Local-Data-for-Better-Health/6vp6-wxuq
),
the top 26 rows are USA crude prevalence  then age adjusted prevalence by
disease (for 2014), or the 13 (either crude prevalence, or age adjusted)
that you are looking for.

Judging by what you were describing above, it doesn't sound like you're
doing a 'things got better or things got worse,
year on year', so perhaps filtering the data to a given year (2014) and
given  health outcome (age adjusted) prior to
downloading might simplify things.
HTH
Chris

On Sun, Apr 23, 2017 at 6:19 PM, Joe Larson <[hidden email]>
wrote:

> Hello,
>  I am getting wrapped around the axial on a shiny project, being a new R
> learner I am asking for some help. I am trying to combine CDC data, Shiny,
> and leaflet together to produce a map that has the top 500 USA cities with
> adverse medical outcomes.  The goal is to allow the user to select either
> the range of outcome, state or adverse outcome.  Through Shiny examples, I
> am able to get the range function to work, but not the "state" or "health
> outcomes" to function correctly. My issue is in using input type selectable
>  "selectInput("state", "State Abbr", data_c$state, selected = NULL,
> multiple = FALSE, selectize = TRUE, width = NULL, size = NULL), "  The
> challenge I am facing is selecting the correct "Render reactive output"
>  and how to mutate the datatable so the users inputs are correctly
> displayed.  If you want the full code, to let me know.
> I have been reading shiny.rstudios.com as well the shiny-cheatsheet.pdf
> with the solution not coming to mind.
> Thanks for any suggestion you can offer.
>
> The datatable has this format (first 2 cities and allThe isisue is I can
> not find out which type of 13 outcomens):
> X year state cityname crudeprevalence populationcount measureid lat long
> healthoutcome metropop
> 1 2014 AL Birmingham 32.6 212237 ARTHRITIS 33.52756638 -86.79881747
> 0.0153 Birmingham(city),
> 212237(popluation), 0.0153(Arthritis %)
> 2 2013 AL Birmingham 45.9 212237 BPHIGH 33.52756638 -86.79881747
> 0.0216 Birmingham(city),
> 212237(popluation), 0.0216(High blood pressure %)
> 3 2014 AL Birmingham 6.1 212237 CANCER 33.52756638 -86.79881747 0.0028
> Birmingham(city),
> 212237(popluation), 0.0028(Cancer (excluding skin cancer)  %)
> 4 2014 AL Birmingham 11.4 212237 CASTHMA 33.52756638 -86.79881747
> 0.0053 Birmingham(city),
> 212237(popluation), 0.0053(Asthma %)
> 5 2014 AL Birmingham 7.6 212237 CHD 33.52756638 -86.79881747 0.0035
> Birmingham(city),
> 212237(popluation), 0.0035(Coronary heart disease %)
> 6 2014 AL Birmingham 9.4 212237 COPD 33.52756638 -86.79881747 0.0044
> Birmingham(city),
> 212237(popluation), 0.0044(Chronic obstructive pulmonary disease %)
> 7 2014 AL Birmingham 16.1 212237 DIABETES 33.52756638 -86.79881747
> 0.0075 Birmingham(city),
> 212237(popluation), 0.0075(Diabetes %)
> 8 2013 AL Birmingham 35.4 212237 HIGHCHOL 33.52756638 -86.79881747
> 0.0166 Birmingham(city),
> 212237(popluation), 0.0166(High cholesterol %)
> 9 2014 AL Birmingham 3.3 212237 KIDNEY 33.52756638 -86.79881747 0.0015
> Birmingham(city),
> 212237(popluation), 0.0015(Kidney disease %)
> 10 2014 AL Birmingham 17 212237 MHLTH 33.52756638 -86.79881747 0.008
> Birmingham(city),
> 212237(popluation), 0.008(Mental health %)
> 11 2014 AL Birmingham 18.3 212237 PHLTH 33.52756638 -86.79881747
> 0.0086 Birmingham(city),
> 212237(popluation), 0.0086(Physical health %)
> 12 2014 AL Birmingham 5 212237 STROKE 33.52756638 -86.79881747 0.0023
> Birmingham(city),
> 212237(popluation), 0.0023(Stroke %)
> 13 2014 AL Birmingham 25.9 212237 TEETHLOST 33.52756638 -86.79881747
> 0.0122 Birmingham(city),
> 212237(popluation), 0.0122(All teeth lost %)
> 14 2014 AL Hoover 25.3 81619 ARTHRITIS 33.37676027 -86.80519376 0.0309
> Hoover(city),
> 81619(popluation), 0.0309(Arthritis %)
> Thanks for you help and insight in advance
> Joe
>
>         [[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: Help for STFDF Object Creation

Sun, 04/23/2017 - 18:20
Thanks for your reply. Sorry for my second post, I tried to be more clear.

Here is link of my complete source code, data and R instance:https://drive.google.com/drive/folders/0B3tdiSZ9hBj4dVZXOXItUjdrRlE?usp=sharing


Hope, it would be easier to understand my problem.

    On Sunday, April 23, 2017 5:25 PM, Edzer Pebesma <[hidden email]> wrote:
 

 Please provide a script of an example that we can reproduce, preferably
using data that is loaded from a package (to safe other people's time),
otherwise it is impossible for anyone to help you, or unlikely somebody
will. Your script does not create count data but assumes it is there;
you don't show which package you loaded; it is unclear what your problem is.

Your second post on this list came through completely scrambled, and
within 24 hours. Don't expect that reposting, or scrambling, helps.

On 23/04/17 06:02, Ankur Sarker via R-sig-Geo wrote:
> Hi,
> I am trying to create a STFDF object and draw variogram. However, I am getting several errors. Here is my r code:
> #create count data
> c0157ts <- as.numeric(c(t(c0157)))
> c0001ts <- as.numeric(c(t(c0001)))
>
> #create coordinates
>  lat <- 0
>  lon <- 0
>  lat[1] <- 34.0793684
>  lon[1] <- -81.1301722
>  lat[2] <- 33.979127
>  lon[2] <- -81.321166398
>  sp <- data.frame(lat,lon)
>  coordinates(sp)=~lon+lat
>  projection(sp)=CRS("+init=epsg:4326")
>
> #create time object
> time_index <- seq(from = as.POSIXct("2017-01-02 01:00", tz = 'UTC'), to = as.POSIXct("2017-03-22 00:00", tz = 'UTC'), by = "hour")
>
> #combine the data
> mydata <- data.frame("count"=c(c0157ts,c0001ts))
>
> #create STFDF object
> stfdf = STFDF(ozoneSP, time_index, mydata)
> #trying to create variogram
>    var <- variogramST(count~1,data=stfdf,tunit="hours",assumeRegular=F,na.omit=T)which gives me following error:0%Error in apply(do.call(cbind, lapply(ret, function(x) x$np)), 1, sum,  :
>      dim(X) must have a positive length
> I guess my STFDF object creation is not correct. Can anyone give us any hint? Which elements in STFDF object is not correct.Here is the description of STFDF object:> str(stfdf)
> Formal class 'STFDF' [package "spacetime"] with 4 slots
>  ..@ data  :'data.frame':    3792 obs. of  1 variable:
>  .. ..$ count: num [1:3792] 40 32 64 40 55 89 71 43 69 44 ...
>  ..@ sp    :Formal class 'SpatialPoints' [package "sp"] with 3 slots
>  .. .. ..@ coords    : num [1:2, 1:2] -9031369 -9052631 4015522 4002120
>  .. .. .. ..- attr(*, "dimnames")=List of 2
>  .. .. .. .. ..$ : NULL
>  .. .. .. .. ..$ : chr [1:2] "lon" "lat"
>  .. .. ..@ bbox      : num [1:2, 1:2] -9052631 4002120 -9031369 4015522
>  .. .. .. ..- attr(*, "dimnames")=List of 2
>  .. .. .. .. ..$ : chr [1:2] "lon" "lat"
>  .. .. .. .. ..$ : chr [1:2] "min" "max"
>  .. .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
>  .. .. .. .. ..@ projargs: chr "+init=epsg:3395 +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
>  ..@ time  :An ‘xts’ object on 2017-01-02 01:00:00/2017-03-22 containing:
>  Data: int [1:1896, 1] 1 2 3 4 5 6 7 8 9 10 ...
>  - attr(*, "dimnames")=List of 2
>  ..$ : NULL
>  ..$ : chr "timeIndex"
>  Indexed by objects of class: [POSIXct,POSIXt] TZ: UTC
>  xts Attributes: 
>  NULL
> Thanks for your help!
> With Regards,Rukna
>     [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
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/
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

   
        [[alternative HTML version deleted]]

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

R-help - Shiny and Datatable

Sun, 04/23/2017 - 17:19
Hello,
 I am getting wrapped around the axial on a shiny project, being a new R
learner I am asking for some help. I am trying to combine CDC data, Shiny,
and leaflet together to produce a map that has the top 500 USA cities with
adverse medical outcomes.  The goal is to allow the user to select either
the range of outcome, state or adverse outcome.  Through Shiny examples, I
am able to get the range function to work, but not the "state" or "health
outcomes" to function correctly. My issue is in using input type selectable
 "selectInput("state", "State Abbr", data_c$state, selected = NULL,
multiple = FALSE, selectize = TRUE, width = NULL, size = NULL), "  The
challenge I am facing is selecting the correct "Render reactive output"
 and how to mutate the datatable so the users inputs are correctly
displayed.  If you want the full code, to let me know.
I have been reading shiny.rstudios.com as well the shiny-cheatsheet.pdf
with the solution not coming to mind.
Thanks for any suggestion you can offer.

The datatable has this format (first 2 cities and allThe isisue is I can
not find out which type of 13 outcomens):
X year state cityname crudeprevalence populationcount measureid lat long
healthoutcome metropop
1 2014 AL Birmingham 32.6 212237 ARTHRITIS 33.52756638 -86.79881747
0.0153 Birmingham(city),
212237(popluation), 0.0153(Arthritis %)
2 2013 AL Birmingham 45.9 212237 BPHIGH 33.52756638 -86.79881747
0.0216 Birmingham(city),
212237(popluation), 0.0216(High blood pressure %)
3 2014 AL Birmingham 6.1 212237 CANCER 33.52756638 -86.79881747 0.0028
Birmingham(city),
212237(popluation), 0.0028(Cancer (excluding skin cancer)  %)
4 2014 AL Birmingham 11.4 212237 CASTHMA 33.52756638 -86.79881747
0.0053 Birmingham(city),
212237(popluation), 0.0053(Asthma %)
5 2014 AL Birmingham 7.6 212237 CHD 33.52756638 -86.79881747 0.0035
Birmingham(city),
212237(popluation), 0.0035(Coronary heart disease %)
6 2014 AL Birmingham 9.4 212237 COPD 33.52756638 -86.79881747 0.0044
Birmingham(city),
212237(popluation), 0.0044(Chronic obstructive pulmonary disease %)
7 2014 AL Birmingham 16.1 212237 DIABETES 33.52756638 -86.79881747
0.0075 Birmingham(city),
212237(popluation), 0.0075(Diabetes %)
8 2013 AL Birmingham 35.4 212237 HIGHCHOL 33.52756638 -86.79881747
0.0166 Birmingham(city),
212237(popluation), 0.0166(High cholesterol %)
9 2014 AL Birmingham 3.3 212237 KIDNEY 33.52756638 -86.79881747 0.0015
Birmingham(city),
212237(popluation), 0.0015(Kidney disease %)
10 2014 AL Birmingham 17 212237 MHLTH 33.52756638 -86.79881747 0.008
Birmingham(city),
212237(popluation), 0.008(Mental health %)
11 2014 AL Birmingham 18.3 212237 PHLTH 33.52756638 -86.79881747
0.0086 Birmingham(city),
212237(popluation), 0.0086(Physical health %)
12 2014 AL Birmingham 5 212237 STROKE 33.52756638 -86.79881747 0.0023
Birmingham(city),
212237(popluation), 0.0023(Stroke %)
13 2014 AL Birmingham 25.9 212237 TEETHLOST 33.52756638 -86.79881747
0.0122 Birmingham(city),
212237(popluation), 0.0122(All teeth lost %)
14 2014 AL Hoover 25.3 81619 ARTHRITIS 33.37676027 -86.80519376 0.0309
Hoover(city),
81619(popluation), 0.0309(Arthritis %)
Thanks for you help and insight in advance
Joe

        [[alternative HTML version deleted]]

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

Re: Help for STFDF Object Creation

Sun, 04/23/2017 - 16:24
Please provide a script of an example that we can reproduce, preferably
using data that is loaded from a package (to safe other people's time),
otherwise it is impossible for anyone to help you, or unlikely somebody
will. Your script does not create count data but assumes it is there;
you don't show which package you loaded; it is unclear what your problem is.

Your second post on this list came through completely scrambled, and
within 24 hours. Don't expect that reposting, or scrambling, helps.

On 23/04/17 06:02, Ankur Sarker via R-sig-Geo wrote:
> Hi,
> I am trying to create a STFDF object and draw variogram. However, I am getting several errors. Here is my r code:
> #create count data
> c0157ts <- as.numeric(c(t(c0157)))
> c0001ts <- as.numeric(c(t(c0001)))
>
> #create coordinates
>   lat <- 0
>   lon <- 0
>   lat[1] <- 34.0793684
>   lon[1] <- -81.1301722
>   lat[2] <- 33.979127
>   lon[2] <- -81.321166398
>   sp <- data.frame(lat,lon)
>   coordinates(sp)=~lon+lat
>   projection(sp)=CRS("+init=epsg:4326")
>
> #create time object
> time_index <- seq(from = as.POSIXct("2017-01-02 01:00", tz = 'UTC'), to = as.POSIXct("2017-03-22 00:00", tz = 'UTC'), by = "hour")
>
> #combine the data
> mydata <- data.frame("count"=c(c0157ts,c0001ts))
>
> #create STFDF object
> stfdf = STFDF(ozoneSP, time_index, mydata)
> #trying to create variogram
>     var <- variogramST(count~1,data=stfdf,tunit="hours",assumeRegular=F,na.omit=T)which gives me following error:0%Error in apply(do.call(cbind, lapply(ret, function(x) x$np)), 1, sum,  :
>       dim(X) must have a positive length
> I guess my STFDF object creation is not correct. Can anyone give us any hint? Which elements in STFDF object is not correct.Here is the description of STFDF object:> str(stfdf)
> Formal class 'STFDF' [package "spacetime"] with 4 slots
>   ..@ data   :'data.frame':     3792 obs. of  1 variable:
>   .. ..$ count: num [1:3792] 40 32 64 40 55 89 71 43 69 44 ...
>   ..@ sp     :Formal class 'SpatialPoints' [package "sp"] with 3 slots
>   .. .. ..@ coords     : num [1:2, 1:2] -9031369 -9052631 4015522 4002120
>   .. .. .. ..- attr(*, "dimnames")=List of 2
>   .. .. .. .. ..$ : NULL
>   .. .. .. .. ..$ : chr [1:2] "lon" "lat"
>   .. .. ..@ bbox       : num [1:2, 1:2] -9052631 4002120 -9031369 4015522
>   .. .. .. ..- attr(*, "dimnames")=List of 2
>   .. .. .. .. ..$ : chr [1:2] "lon" "lat"
>   .. .. .. .. ..$ : chr [1:2] "min" "max"
>   .. .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
>   .. .. .. .. ..@ projargs: chr "+init=epsg:3395 +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
>   ..@ time   :An ‘xts’ object on 2017-01-02 01:00:00/2017-03-22 containing:
>   Data: int [1:1896, 1] 1 2 3 4 5 6 7 8 9 10 ...
>  - attr(*, "dimnames")=List of 2
>   ..$ : NULL
>   ..$ : chr "timeIndex"
>   Indexed by objects of class: [POSIXct,POSIXt] TZ: UTC
>   xts Attributes:  
>  NULL
> Thanks for your help!
> With Regards,Rukna
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
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/


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

STFDF Object Creation and Variogram

Sun, 04/23/2017 - 16:00
Hi,I am trying to create a STFDF object and draw variogram. However, I am getting several errors. 
Here is my r code:
    #create count data    c0157ts <- as.numeric(c(t(c0157)))    c0001ts <- as.numeric(c(t(c0001)))
    #create coordinates       lat <- 0      lon <- 0      lat[1] <- 34.0793684      lon[1] <- -81.1301722      lat[2] <- 33.979127      lon[2] <- -81.321166398      sp <- data.frame(lat,lon)      coordinates(sp)=~lon+lat      projection(sp)=CRS("+init=epsg:4326")
    #create time object       time_index <- seq(from = as.POSIXct("2017-01-02 01:00", tz = 'UTC'), to = as.POSIXct("2017-03-22 00:00", tz = 'UTC'), by = "hour")
    #combine the data       mydata <- data.frame("count"=c(c0157ts,c0001ts)) 
    #create STFDF object    stfdf = STFDF(ozoneSP, time_index, mydata)
    #trying to create variogram    var <- variogramST(count~1,data=stfdf,tunit="hours",assumeRegular=F,na.omit=T)

which gives me following error:    Error in apply(do.call(cbind, lapply(ret, function(x) x$np)), 1, sum,  : dim(X) must have a positive length

I guess my STFDF object creation is not correct. Can anyone give me any hint? Which elements in STFDF object is not correct. Is it because of my time object?

Here is the description of STFDF object:str(stfdf)
Output:
    Formal class 'STFDF' [package "spacetime"] with 4 slots      ..@ data  :'data.frame':    3792 obs. of  1 variable:      .. ..$ count: num [1:3792] 40 32 64 40 55 89 71 43 69 44 ...      ..@ sp    :Formal class 'SpatialPoints' [package "sp"] with 3 slots      .. .. ..@ coords    : num [1:2, 1:2] -9031369 -9052631 4015522 4002120      .. .. .. ..- attr(*, "dimnames")=List of 2      .. .. .. .. ..$ : NULL      .. .. .. .. ..$ : chr [1:2] "lon" "lat"      .. .. ..@ bbox      : num [1:2, 1:2] -9052631 4002120 -9031369 4015522      .. .. .. ..- attr(*, "dimnames")=List of 2      .. .. .. .. ..$ : chr [1:2] "lon" "lat"      .. .. .. .. ..$ : chr [1:2] "min" "max"      .. .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot      .. .. .. .. ..@ projargs: chr "+init=epsg:3395 +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"      ..@ time  :An ‘xts’ object on 2017-01-02 01:00:00/2017-03-22 containing:      Data: int [1:1896, 1] 1 2 3 4 5 6 7 8 9 10 ...    - attr(*, "dimnames")=List of 2      ..$ : NULL      ..$ : chr "timeIndex"      Indexed by objects of class: [POSIXct,POSIXt] TZ: UTC      xts Attributes:      NULL


Here, I have attached the link of Rdat for your better understanding. Thanks for your help!
STFDF.RData

 
|  
|    |  
STFDF.RData
   |  |

  |

 

        [[alternative HTML version deleted]]

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

Pycnophylactic Interpolation

Sun, 04/23/2017 - 15:35
Hi,

I'm a student of Information Systems in Brazil and a junior software archtect. To my final work on the undergraduate program I'm implementing a software with the Potential model. After calculate the potential values to each point of analysis I would like to be able to create a map to represent the result. The map should looks like this one: https://www.r-bloggers.com/contour-and-density-layers-with-ggmap/. I'm new in R so I have no a big knowhow and a lot of doubts. I would like ask for community's help. My problem is that I have no to many points like the example so the map looks incorrect, actualy I have a few points (x,y and z). I think I must Interpolate these points to get a good result, I would like to use the Pycnophylactic Interpolation avaliable in pycno package.

Here is a example of data I obtain (no interpolation applied):


City    Area    Score   latitudes       longitudes      Potentials
Belo Horizone   330,95  2513451 -19,91663263    -43,93337198    164874,1295
Betim   346     412003  -19,96730327    -44,20119724    82629,32514
Contagem        195,268 648766  -19,91615326    -44,08087722    140851,9626
Sarzedo 61,892  25798   -20,06250329    -44,11504723    79307,30623
Rio Acima       50      90      -20,08830329    -43,79104717    55226,91779


Please let me know if you need any extra information.
Thank you for your patience and cooperation.
Breno Oliveira


        [[alternative HTML version deleted]]

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

Help for STFDF Object Creation

Sat, 04/22/2017 - 23:02
Hi,
I am trying to create a STFDF object and draw variogram. However, I am getting several errors. Here is my r code:
#create count data
c0157ts <- as.numeric(c(t(c0157)))
c0001ts <- as.numeric(c(t(c0001)))

#create coordinates
  lat <- 0
  lon <- 0
  lat[1] <- 34.0793684
  lon[1] <- -81.1301722
  lat[2] <- 33.979127
  lon[2] <- -81.321166398
  sp <- data.frame(lat,lon)
  coordinates(sp)=~lon+lat
  projection(sp)=CRS("+init=epsg:4326")

#create time object
time_index <- seq(from = as.POSIXct("2017-01-02 01:00", tz = 'UTC'), to = as.POSIXct("2017-03-22 00:00", tz = 'UTC'), by = "hour")

#combine the data
mydata <- data.frame("count"=c(c0157ts,c0001ts))

#create STFDF object
stfdf = STFDF(ozoneSP, time_index, mydata)
#trying to create variogram
    var <- variogramST(count~1,data=stfdf,tunit="hours",assumeRegular=F,na.omit=T)which gives me following error:0%Error in apply(do.call(cbind, lapply(ret, function(x) x$np)), 1, sum,  :
      dim(X) must have a positive length
I guess my STFDF object creation is not correct. Can anyone give us any hint? Which elements in STFDF object is not correct.Here is the description of STFDF object:> str(stfdf)
Formal class 'STFDF' [package "spacetime"] with 4 slots
  ..@ data   :'data.frame':     3792 obs. of  1 variable:
  .. ..$ count: num [1:3792] 40 32 64 40 55 89 71 43 69 44 ...
  ..@ sp     :Formal class 'SpatialPoints' [package "sp"] with 3 slots
  .. .. ..@ coords     : num [1:2, 1:2] -9031369 -9052631 4015522 4002120
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : NULL
  .. .. .. .. ..$ : chr [1:2] "lon" "lat"
  .. .. ..@ bbox       : num [1:2, 1:2] -9052631 4002120 -9031369 4015522
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "lon" "lat"
  .. .. .. .. ..$ : chr [1:2] "min" "max"
  .. .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. .. .. ..@ projargs: chr "+init=epsg:3395 +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
  ..@ time   :An ‘xts’ object on 2017-01-02 01:00:00/2017-03-22 containing:
  Data: int [1:1896, 1] 1 2 3 4 5 6 7 8 9 10 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr "timeIndex"
  Indexed by objects of class: [POSIXct,POSIXt] TZ: UTC
  xts Attributes:  
 NULL
Thanks for your help!
With Regards,Rukna
        [[alternative HTML version deleted]]

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

Pages