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

Match European cities to Coordinates/NUTS 2

Mon, 11/05/2018 - 15:25
I have a dataframe (more than 50,000 observations), where one of the
variables is the city that unit is location. The units are all in Europe.

My goal is to assign NUTS-2 code to each of these cities. However, I am not
aware of any direct way of achieving this, so I wanted to first assign
coordinates to the cities and then use the 'over' function to match with
NUTS regions from EU shapefile.

I have tried to use 'geocode' from the ggmap package but there is a 2,500
per day limit. Is there any other solution? Any help will be greatly
appreciated.

Sincerely,

Milu

        [[alternative HTML version deleted]]

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

Re: Raster to rasterbrick does not preserve Date/Time

Mon, 11/05/2018 - 14:59
Thank you, Ben! This works well.

Sincerely,

Milu

On Sun, Nov 4, 2018 at 11:08 PM Ben Tupper <[hidden email]> wrote:

> Hi,
>
> I'm not sure about reading the datetime stamp directly form the file, but
> you have enough info to convert the z-value which is well described in your
> ncdf.  I'm not familiar with the '%.f' format specifier, but it looks like
> it is a fractional day.  Something like this perhaps will work?
>
> z      <- raster::getZ(rbick)
> z0     <- floor(z)
> zf     <- z - z0
> dt     <- as.POSIXct(as.character(z0), format = "%Y%m%d") + 24*60*60*zf
> rbrick <- raster::setZ(rbrick, z, name = 'datetime')
>
> You'll need to think through the timezone issues that might arise.
> as.POSIXct() has arguments that can help you with that.
>
> Cheers,
> Ben
>
>
> > On Nov 4, 2018, at 12:39 PM, Miluji Sb <[hidden email]> wrote:
> >
> > Dear all,
> >
> > I have a raster with multiple variables and date/time information;
> >
> >      time  Size:17164   *** is unlimited ***
> >            standard_name: time
> >            bounds: time_bnds
> >            units: day as %Y%m%d.%f
> >            calendar: proleptic_gregorian
> >        bnds  Size:2
> >
> > My goal is to convert the raster file to rasterbrick for extraction
> > purposes;
> >
> > rbrick <- brick(r,  values=TRUE, varname="var1")
> >
> > However, after conversion - the date/time information is not preserved.
> How
> > can I preserve this information?
> >
> > class       : RasterBrick
> > dimensions  : 600, 1440, 864000, 17164  (nrow, ncol, ncell, nlayers)
> > resolution  : 0.25, 0.25  (x, y)
> > extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
> > coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> > data source : /work/ncfile.nc4
> > names       : X19700101.4375, X19700102.4375, X19700103.4375,
> > X19700104.4375, X19700105.4375, X19700106.4375, X19700107.4375,
> > X19700108.4375, X19700109.4375, X19700110.4375, X19700111.4375,
> > X19700112.4375, X19700113.4375, X19700114.4375, X19700115.4375, ...
> > z-value     : 19700101.4375, 20161231.4375 (min, max)
> > varname     : var1
> >
> > Any help will be greatly appreciated!
> >
> > Sincerely,
> >
> > Milu
> >
> >       [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > [hidden email]
> > 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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Reordering geographical coordinates clockwise to make a polygon

Mon, 11/05/2018 - 13:41
alphahull or a hand-crafted "triangulate and cull long/large triangles"
might be workable options (geometry and RTriangle have the fastest and most
useful triangulation functions).

It's not generally a tractable problem afaics.  Some refs that came up in
recent discussions:

https://pdfs.semanticscholar.org/669b/1415cd64d39f1e44fd90c2b9d9453c42581b.pdf

http://erikdemaine.org/polygonization/

Cheers, Mike.

On Tue, 6 Nov 2018 at 05:26 Edzer Pebesma <[hidden email]>
wrote:

> Sorry for my previous answer, I didn't get the question fully.
>
> I'm afraid I still don't get the question fully, but functions that
> might help are sf::st_line_merge (creates a LINESTRING from line pieces)
> and sf::st_polygonize (creates a polygon from a LINESTRING that forms a
> closed ring)
>
> On 11/5/18 5:35 PM, Patrick Giraudoux wrote:
> >
> > Apologize to answer to myself: the way described below would work only
> > if there are no concave "peninsula" towards north or south inside the
> > polygon. Even thinking about an alternate solution, e.g. the Traveling
> > Salesman Problem (TSP), ie. the shortest route linking all points, one
> > could get wrong since points of  a narrowpeninsula could have points
> > closest to the opposite border than from the next point on the border...
> >
> > Suppose we are stuck, and should redraw polygons by hand
> >
> >
> > Le 05/11/2018 à 14:02, Patrick Giraudoux a écrit :
> >>
> >> Dear listers,
> >>
> >> There is an interesting post here:
> >>
> https://stackoverflow.com/questions/6989100/sort-points-in-clockwise-order
> >> dealing on the issue. However, I would like to know if a function has
> >> been already developped in a R package.
> >>
> >> I am coping with a young colleague's shapefile where borders of
> >> polygons have been drawn as lines quite inconsistently. To change this
> >> SpatialLinesDataFrame into a SpatialPoints object is easy. The idea is
> >> to  select the points  bordering each polygon (delete the others),
> >> define the point set as a Polygon, then rebuild step by step a
> >> SpatialPolygonsDataFrame with all its (~25) polygons. It would be much
> >> quicker than to redraw one by one the 5160 data points (two times on
> >> borders shared by two polygons).
> >>
> >> The problem is that the points must be reordered clockwise (the way
> >> lines making the  borders is all except clockwise) before making them
> >> a polygon.
> >>
> >> Any function already developped for that ?
> >>
> >> Cheers,
> >>
> >> Patrick
> >>
> >>
> >>
> >
> >
> >       [[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
> Heisenbergstrasse 2, 48151 Muenster, Germany
> Phone: +49 251 8333081 <+49%20251%208333081>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

        [[alternative HTML version deleted]]

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

Re: Reordering geographical coordinates clockwise to make a polygon

Mon, 11/05/2018 - 12:25
Sorry for my previous answer, I didn't get the question fully.

I'm afraid I still don't get the question fully, but functions that
might help are sf::st_line_merge (creates a LINESTRING from line pieces)
and sf::st_polygonize (creates a polygon from a LINESTRING that forms a
closed ring)

On 11/5/18 5:35 PM, Patrick Giraudoux wrote:
>
> Apologize to answer to myself: the way described below would work only
> if there are no concave "peninsula" towards north or south inside the
> polygon. Even thinking about an alternate solution, e.g. the Traveling
> Salesman Problem (TSP), ie. the shortest route linking all points, one
> could get wrong since points of  a narrowpeninsula could have points
> closest to the opposite border than from the next point on the border...
>
> Suppose we are stuck, and should redraw polygons by hand
>
>
> Le 05/11/2018 à 14:02, Patrick Giraudoux a écrit :
>>
>> Dear listers,
>>
>> There is an interesting post here:
>> https://stackoverflow.com/questions/6989100/sort-points-in-clockwise-order 
>> dealing on the issue. However, I would like to know if a function has
>> been already developped in a R package.
>>
>> I am coping with a young colleague's shapefile where borders of
>> polygons have been drawn as lines quite inconsistently. To change this
>> SpatialLinesDataFrame into a SpatialPoints object is easy. The idea is
>> to  select the points  bordering each polygon (delete the others),
>> define the point set as a Polygon, then rebuild step by step a
>> SpatialPolygonsDataFrame with all its (~25) polygons. It would be much
>> quicker than to redraw one by one the 5160 data points (two times on
>> borders shared by two polygons).
>>
>> The problem is that the points must be reordered clockwise (the way
>> lines making the  borders is all except clockwise) before making them
>> a polygon.
>>
>> Any function already developped for that ?
>>
>> Cheers,
>>
>> Patrick
>>
>>
>>
>
>
> [[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
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

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

pEpkey.asc (2K) Download Attachment

Re: Reordering geographical coordinates clockwise to make a polygon

Mon, 11/05/2018 - 10:35

Apologize to answer to myself: the way described below would work only
if there are no concave "peninsula" towards north or south inside the
polygon. Even thinking about an alternate solution, e.g. the Traveling
Salesman Problem (TSP), ie. the shortest route linking all points, one
could get wrong since points of  a narrowpeninsula could have points
closest to the opposite border than from the next point on the border...

Suppose we are stuck, and should redraw polygons by hand


Le 05/11/2018 à 14:02, Patrick Giraudoux a écrit :
>
> Dear listers,
>
> There is an interesting post here:
> https://stackoverflow.com/questions/6989100/sort-points-in-clockwise-order 
> dealing on the issue. However, I would like to know if a function has
> been already developped in a R package.
>
> I am coping with a young colleague's shapefile where borders of
> polygons have been drawn as lines quite inconsistently. To change this
> SpatialLinesDataFrame into a SpatialPoints object is easy. The idea is
> to  select the points  bordering each polygon (delete the others),
> define the point set as a Polygon, then rebuild step by step a
> SpatialPolygonsDataFrame with all its (~25) polygons. It would be much
> quicker than to redraw one by one the 5160 data points (two times on
> borders shared by two polygons).
>
> The problem is that the points must be reordered clockwise (the way
> lines making the  borders is all except clockwise) before making them
> a polygon.
>
> Any function already developped for that ?
>
> Cheers,
>
> Patrick
>
>
>

        [[alternative HTML version deleted]]

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

Re: Reordering geographical coordinates clockwise to make a polygon

Mon, 11/05/2018 - 09:31
Look for argument check_ring_dir in the documentation of sf::st_read.

On 11/5/18 2:02 PM, Patrick Giraudoux wrote:
> Dear listers,
>
> There is an interesting post here:
> https://stackoverflow.com/questions/6989100/sort-points-in-clockwise-order 
> dealing on the issue. However, I would like to know if a function has
> been already developped in a R package.
>
> I am coping with a young colleague's shapefile where borders of polygons
> have been drawn as lines quite inconsistently. To change this
> SpatialLinesDataFrame into a SpatialPoints object is easy. The idea is
> to  select the points  bordering each polygon (delete the others),
> define the point set as a Polygon, then rebuild step by step a
> SpatialPolygonsDataFrame with all its (~25) polygons. It would be much
> quicker than to redraw one by one the 5160 data points (two times on
> borders shared by two polygons).
>
> The problem is that the points must be reordered clockwise (the way
> lines making the  borders is all except clockwise) before making them a
> polygon.
>
> Any function already developped for that ?
>
> Cheers,
>
> Patrick
>
>
>
>
> [[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
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

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

pEpkey.asc (2K) Download Attachment

Reordering geographical coordinates clockwise to make a polygon

Mon, 11/05/2018 - 07:02
Dear listers,

There is an interesting post here:
https://stackoverflow.com/questions/6989100/sort-points-in-clockwise-order 
dealing on the issue. However, I would like to know if a function has
been already developped in a R package.

I am coping with a young colleague's shapefile where borders of polygons
have been drawn as lines quite inconsistently. To change this
SpatialLinesDataFrame into a SpatialPoints object is easy. The idea is
to  select the points  bordering each polygon (delete the others),
define the point set as a Polygon, then rebuild step by step a
SpatialPolygonsDataFrame with all its (~25) polygons. It would be much
quicker than to redraw one by one the 5160 data points (two times on
borders shared by two polygons).

The problem is that the points must be reordered clockwise (the way
lines making the  borders is all except clockwise) before making them a
polygon.

Any function already developped for that ?

Cheers,

Patrick




        [[alternative HTML version deleted]]

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

Re: Raster to rasterbrick does not preserve Date/Time

Sun, 11/04/2018 - 16:08
Hi,

I'm not sure about reading the datetime stamp directly form the file, but you have enough info to convert the z-value which is well described in your ncdf.  I'm not familiar with the '%.f' format specifier, but it looks like it is a fractional day.  Something like this perhaps will work?

z      <- raster::getZ(rbick)
z0     <- floor(z)
zf     <- z - z0
dt     <- as.POSIXct(as.character(z0), format = "%Y%m%d") + 24*60*60*zf
rbrick <- raster::setZ(rbrick, z, name = 'datetime')

You'll need to think through the timezone issues that might arise.  as.POSIXct() has arguments that can help you with that.

Cheers,
Ben


> On Nov 4, 2018, at 12:39 PM, Miluji Sb <[hidden email]> wrote:
>
> Dear all,
>
> I have a raster with multiple variables and date/time information;
>
>      time  Size:17164   *** is unlimited ***
>            standard_name: time
>            bounds: time_bnds
>            units: day as %Y%m%d.%f
>            calendar: proleptic_gregorian
>        bnds  Size:2
>
> My goal is to convert the raster file to rasterbrick for extraction
> purposes;
>
> rbrick <- brick(r,  values=TRUE, varname="var1")
>
> However, after conversion - the date/time information is not preserved. How
> can I preserve this information?
>
> class       : RasterBrick
> dimensions  : 600, 1440, 864000, 17164  (nrow, ncol, ncell, nlayers)
> resolution  : 0.25, 0.25  (x, y)
> extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> data source : /work/ncfile.nc4
> names       : X19700101.4375, X19700102.4375, X19700103.4375,
> X19700104.4375, X19700105.4375, X19700106.4375, X19700107.4375,
> X19700108.4375, X19700109.4375, X19700110.4375, X19700111.4375,
> X19700112.4375, X19700113.4375, X19700114.4375, X19700115.4375, ...
> z-value     : 19700101.4375, 20161231.4375 (min, max)
> varname     : var1
>
> Any help will be greatly appreciated!
>
> Sincerely,
>
> Milu
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> 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/

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

Raster to rasterbrick does not preserve Date/Time

Sun, 11/04/2018 - 11:39
Dear all,

I have a raster with multiple variables and date/time information;

      time  Size:17164   *** is unlimited ***
            standard_name: time
            bounds: time_bnds
            units: day as %Y%m%d.%f
            calendar: proleptic_gregorian
        bnds  Size:2

My goal is to convert the raster file to rasterbrick for extraction
purposes;

rbrick <- brick(r,  values=TRUE, varname="var1")

However, after conversion - the date/time information is not preserved. How
can I preserve this information?

class       : RasterBrick
dimensions  : 600, 1440, 864000, 17164  (nrow, ncol, ncell, nlayers)
resolution  : 0.25, 0.25  (x, y)
extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : /work/ncfile.nc4
names       : X19700101.4375, X19700102.4375, X19700103.4375,
X19700104.4375, X19700105.4375, X19700106.4375, X19700107.4375,
X19700108.4375, X19700109.4375, X19700110.4375, X19700111.4375,
X19700112.4375, X19700113.4375, X19700114.4375, X19700115.4375, ...
z-value     : 19700101.4375, 20161231.4375 (min, max)
varname     : var1

Any help will be greatly appreciated!

Sincerely,

Milu

        [[alternative HTML version deleted]]

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

Re: Mean every 12 layers (months)

Sat, 11/03/2018 - 18:44
Hi Jackson,

Regarding your first question, your indices vector has 6000X12 elements,
since you repeat each number 12 times.
You could maybe do:
indices <- rep(1:500,each=12), so you end up with 6000 elements

For the second question, one way you could do that is to add the zoo
package and do the following:
library(zoo)
TS20@z$date <- as.yearmon(...) # give a vector with the dates of the layers
with this way you will have the date of each layer in a 'date' format and
thus you will be able to easily select the layers you want on a later stage
(eg, choose specific months, specific years, etc).
Otherwise, you could use:
names(TS20) <- ... # give a vector with the new layer names

Hope the info can be helpful to you.

Thanks,
Nikos

*-------------------------------------------------------*
*Nikolaos Mastrantonas*

*Phone: *+44 (0) 7599 843061
*E-mail: [hidden email] <[hidden email]>*
*Linkedin*: nikolaosmastrantonas
<http://www.linkedin.com/in/nikolaosmastrantonas>


On Sun, 4 Nov 2018 at 00:14, Jackson Rodrigues <[hidden email]>
wrote:

> Hi fellows,
>
> My name is Jackson.
>
> I have 2 questions. I am gratful for any help for anyone.
>
> ----------------------
> *Question 1*
> Regarding extracting multiple mean from long raster file.
>  I want to reduce 6000 layers (1 per month) to 500 (1 per year)
>
> library(raster)
> >TS20 = brick("trace.20.10200-09701BP.cam2.h0.TS.nc", level=1, var="TS")
> #It extents from 10200 BP to 9701 BP- a total of 6000 layers (1 per month)
> >TS20 = rotate(TS20)
> >TS20
> class       : RasterBrick
> dimensions  : 48, 96, 4608, 6000  (nrow, ncol, ncell, nlayers)
> resolution  : 3.75, 3.708898  (x, y)
> extent      : -178.125, 181.875, -89.01354, 89.01354  (xmin, xmax, ymin,
> ymax)
>
> >indices<-rep(1:6000,each=12) #found somewhere 6000 layers / 12 months =
> 500 annual layers
> >s.mean<-stackApply(TS20, indices, fun = mean)
> Warning message:
> In ind[] <- indices :
>   number of items to replace is not a multiple of replacement length
>
> #Why do I get this warning message?  My raster has 6000 layers, so my
> indices should work, or not?
> ---------------------------
> *Question 2*
> how could rename my layers combining years and months ?
> For example?
> Layer 1 should be TS20$10200.Jan
> Layer 2 should be TS20$10200.Feb
> ...
> Layer 6000 should be TS20$9701.Dec
>
> In the end I will calculate seasonal mean. So, this is why I would like to
> rename them.
>
> Thank you all,
>
> Jackson
>
>         [[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

Mean every 12 layers (months)

Sat, 11/03/2018 - 18:13
Hi fellows,

My name is Jackson.

I have 2 questions. I am gratful for any help for anyone.

----------------------
*Question 1*
Regarding extracting multiple mean from long raster file.
 I want to reduce 6000 layers (1 per month) to 500 (1 per year)

library(raster)
>TS20 = brick("trace.20.10200-09701BP.cam2.h0.TS.nc", level=1, var="TS")
#It extents from 10200 BP to 9701 BP- a total of 6000 layers (1 per month)
>TS20 = rotate(TS20)
>TS20
class       : RasterBrick
dimensions  : 48, 96, 4608, 6000  (nrow, ncol, ncell, nlayers)
resolution  : 3.75, 3.708898  (x, y)
extent      : -178.125, 181.875, -89.01354, 89.01354  (xmin, xmax, ymin,
ymax)

>indices<-rep(1:6000,each=12) #found somewhere 6000 layers / 12 months =
500 annual layers
>s.mean<-stackApply(TS20, indices, fun = mean)
Warning message:
In ind[] <- indices :
  number of items to replace is not a multiple of replacement length

#Why do I get this warning message?  My raster has 6000 layers, so my
indices should work, or not?
---------------------------
*Question 2*
how could rename my layers combining years and months ?
For example?
Layer 1 should be TS20$10200.Jan
Layer 2 should be TS20$10200.Feb
...
Layer 6000 should be TS20$9701.Dec

In the end I will calculate seasonal mean. So, this is why I would like to
rename them.

Thank you all,

Jackson

        [[alternative HTML version deleted]]

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

Re: Cost distance in spatial weight matrix

Mon, 10/29/2018 - 07:48
Dear Roozbeh,

package geoRcb [1]* provides function distmatGen() [2], which computes a
cost-based distance matrix between any two sets of locations, given a
raster cost surface (such as elevation). The function uses package
gdistance to actually perform the computations, and it does not depend
on anything else in the package, so you can simply copy it to your code
base.

Is that what you were looking for?

ƒacu.-



[1] https://github.com/famuvie/geoRcb/

[2] https://github.com/famuvie/geoRcb/blob/master/R/distmatGen.R

* Disclaimer: I'm the author.


On 10/29/18 1:42 PM, Roozbeh Valavi wrote:
> Dear list,
>
> I wonder if it is possible to add some kind of cost/weighting when creating
> a spatial weight matrix to be used in a spatial lag model? I want to
> consider the effect of elevation as a cost for calculating distance to the
> neighbours.
> I am using spdep for creating the W and splm for the model.
>
> Thank you in advance,
> Roozbeh
        [[alternative HTML version deleted]]

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

Cost distance in spatial weight matrix

Mon, 10/29/2018 - 06:42
Dear list,

I wonder if it is possible to add some kind of cost/weighting when creating
a spatial weight matrix to be used in a spatial lag model? I want to
consider the effect of elevation as a cost for calculating distance to the
neighbours.
I am using spdep for creating the W and splm for the model.

Thank you in advance,
Roozbeh
--
*Roozbeh Valavi*
PhD Candidate
The Quantitative & Applied Ecology Group <http://qaeco.com/>
School of BioSciences | Faculty of Science
The University of Melbourne, VIC 3010, Australia
Mobile: +61 423 283 238

        [[alternative HTML version deleted]]

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

Re: Sf find number of parts in multigeometry

Thu, 10/25/2018 - 21:37
Dear all,
Thank you for the hints – I will test and report any surprises ☺
Martin

From: Michael Sumner <[hidden email]>
Date: Tuesday, 23 October 2018 at 5:59 pm
To: Martin Tomko <[hidden email]>
Cc: "[hidden email]" <[hidden email]>
Subject: Re: [R-sig-Geo] Sf find number of parts in multigeometry

I see perfectly good answers to this, but can't resist sharing my own approach.

 I use  gibble::gibble for a summary of parts. See how the multi-part object is number 4, and has 3 subobjects (subobject will repeat within object for holes).

library(sf)
x <- read_sf(system.file("gpkg/nc.gpkg", package = "sf"))
gibble::gibble(x)
#    A tibble: 108 x 5
#     nrow  ncol type         subobject object
#    <int> <int> <chr>            <int>  <int>
#1    27     2 MULTIPOLYGON         1      1
#2    26     2 MULTIPOLYGON         1      2
#3    28     2 MULTIPOLYGON         1      3
#4    26     2 MULTIPOLYGON         1      4
#5     7     2 MULTIPOLYGON         2      4
#6     5     2 MULTIPOLYGON         3      4
#7    34     2 MULTIPOLYGON         1      5
#...

(I use that for mapping out set-based operations on geometry data, it doesn't make a huge amount of sense on its own. I suppose a rapply scheme could be constructed to pluck out things, but you'd also want path extents and sizes and so forth for greater control).

But, if the biggest area of each multipolygon is what you want, give each a unique ID,  use st_cast to POLYGON, group by parent and slice out the largest.

library(dplyr)

x %>% mutate(ID = row_number()) %>% st_cast("POLYGON") %>%
group_by(ID) %>% arrange(desc(st_area(.))) %>% slice(1) %>% ungroup()

Note that holes within a part might reduce the area logic, but so will details of the map projection in use and so on. It's helpful to learn the structure of the underlying geometry lists to craft your own rogue solutions.

HTH

On Tue, Oct 23, 2018, 13:32 Martin Tomko <[hidden email]<mailto:[hidden email]>> wrote:
I am looking for an equivalent to Postgis ST_NumGeometries https://postgis.net/docs/ST_NumGeometries.html<https://postgis.net/docs/ST_NumGeometries.html>
I have multipolygons in an sf df, where most of them are likely to be single part. I want to identify those that are not single part, and possibly retain only the largest polygon part, and cast all into Polygon.
Any advice is appreciated, thanks!

Martin


        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]<mailto:[hidden email]>
https://stat.ethz.ch/mailman/listinfo/r-sig-geo<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
--
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

        [[alternative HTML version deleted]]

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

Re: R: logLik value of the spml spatial panel model

Thu, 10/25/2018 - 09:54
This is great response. Thank you all, Giovanni, Roger and Jeremie.
Based on the response, it seems that the LL element for fixed effect
model is only useful for optimization, but not cross model comparison?
The fixed effect spatial lag panel model, however, does produce the LL
value. Does that mean this LL value shall not be used for model
comparison purpose as well?

Thank you all.

Best,

Danlin


On 2018/10/25 8:46, Millo Giovanni wrote:
> Hello. Thanks for keeping us posted.
>
> As was aptly observed, while RE and pooled models do actually have a logLik element, FE ones don't but there's a reason. While it is not difficult to output the  LL for the FE models too, like Jeremie did (or perhaps even setting quiet=FALSE and taking the last loglikelihood value from optimization), we never managed to decide whether outputting this is appropriate or not. In fact the LL element from the procedure is relative to the "pooled" model on demeaned data, so it is sort of a crossbred likelihood; most people asking for it want it to compare to the LL of the random effects model or similar uses, which I would not advise.
>
> Gianfranco (author of the FE stuff) do please step in if appropriate.
> Best,
> Giovanni
>
> -----Messaggio originale-----
> Da: Jeremie Juste [mailto:[hidden email]]
> Inviato: giovedì 25 ottobre 2018 11:21
> A: Roger Bivand
> Cc: Danlin Yu; [hidden email]; Millo Giovanni; Gianfranco Piras
> Oggetto: Re: [R-sig-Geo] logLik value of the spml spatial panel model
>
>
> Hello,
>
> Thanks for the feedback.
>
> Here is the diff. It is my first one so I'm not sure I'm doing everything right. :-)
>
>
> modified   R/likelihoodsFE.R
> @@ -459,7 +459,7 @@ if(Hess) asyv <- NULL  else asyv <- asyv
>
>
> -       return<-list(coeff=betas, rho = rho, s2 = s2, rest.se = rest.se, rho.se = rho.se, s2.se = s2.se, asyvar1=asyvar1, residuals = r, asyv = asyv)
> +    return<-list(coeff=betas, rho = rho, s2 = s2, rest.se = rest.se,
> + rho.se = rho.se, s2.se = s2.se, ll=LL, asyvar1=asyvar1, residuals = r,
> + asyv = asyv)
>   }
>
> I have compressed splm with the modification with the .git inside. I have installed the modified package again and checked that the modification worked.
>
>> fespaterr <- spml(fm, data = Produc, listw = mat2listw(usaww),
>> model="within", spatial.error="b", Hess = FALSE) fespaterr$logLik
>
>> It may be,however, that there was a reason for omitting it, so checks
>> to verify that it actually is a logLik value on the same basis as
>> other logLik values (not from an abbreviated log-likelihood function
>> used for fitting by omitting constants) would be helpful. It would be
>> good to provide logLik methods for the ML fitted model objects.
> I have checked that it is the logLik that correspond to the estimation of the rho.
>
> Best regards,
>
> Jeremie
>
>
>
>
> Ai sensi del Reg. UE 2016/679 e normativa vigente si precisa che le informazioni contenute in questo messaggio sono riservate ed a uso esclusivo del destinatario. Qualora il messaggio in parola Le fosse pervenuto per errore, La invitiamo ad eliminarlo senza copiarlo e a non inoltrarlo a terzi, dandocene gentilmente comunicazione. Grazie.
>
> Pursuant to Italian Law and EU Reg. 2016/679, you are hereby informed that this message contains confidential information intended only for the use of the addressee. If you are not the addressee, and have received this message by mistake, please delete it and immediately notify us. You may not copy or disseminate this message to anyone. Thank you.
--
___________________________________________
Danlin Yu, Ph.D.
Professor of GIS and Urban Geography
Department of Earth & Environmental Studies
Montclair State University
Montclair, NJ, 07043
Tel: 973-655-4313
Fax: 973-655-4072
Office: CELS 314
Email: [hidden email]
webpage: csam.montclair.edu/~yu

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

Re: where is my raster stack? filename="" and inMemory=FALSE

Thu, 10/25/2018 - 06:00
For a RasterStack the internal organisation is a list of RasterLayers.
Hence you can retrieve the filenames like this, for example:

sapply(1:nlayers(raster), function(f) filename(raster[[f]]))


On 10/25/18 12:36 PM, Hugo Costa wrote:
> Dear list,
>
> I imported a raster stack and made some processing (crop, reclassify). I
> can plot the stack and all looks fine. However, I don't know any more where
> the file is. Namely
>
> fromDisk(raster) retrieves FALSE
> filename(raster) retrieves ""
> inMemory(raster) retrieves FALSE
> hasValues(raster) retrieves TRUE
>
> printing out the raster gives me this:
> class       : RasterStack
> dimensions  : 7041, 7221, 50843061, 11  (nrow, ncol, ncell, nlayers)
> resolution  : 0.00015, 0.00015  (x, y)
> extent      : -10.04175, -8.9586, 37.972, 39.02815  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
> +towgs84=0,0,0
> names       : X10, X20, X30, X40, X50, X60, X70, X80, X90, X100, layer
> min values  :   0,   0,   0,   0,   0,   0,   0,   0,   0,    0,     0
> max values  :   1,   1,   1,   1,   0,   1,   0,   1,   1,    0,     1
>
> I want to do everything in memory, but now I need to use gdalwarp, and
> don't know how to define the argument srcfile. I can write my stack to
> somewhere using writeRaster and get a filename, but that is too slow for my
> needs and I'm looking for an alternative.
>
> Thank you
> Hugo
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Benjamin Leutner

Department of Remote Sensing
University of Würzburg

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

where is my raster stack? filename="" and inMemory=FALSE

Thu, 10/25/2018 - 05:36
Dear list,

I imported a raster stack and made some processing (crop, reclassify). I
can plot the stack and all looks fine. However, I don't know any more where
the file is. Namely

fromDisk(raster) retrieves FALSE
filename(raster) retrieves ""
inMemory(raster) retrieves FALSE
hasValues(raster) retrieves TRUE

printing out the raster gives me this:
class       : RasterStack
dimensions  : 7041, 7221, 50843061, 11  (nrow, ncol, ncell, nlayers)
resolution  : 0.00015, 0.00015  (x, y)
extent      : -10.04175, -8.9586, 37.972, 39.02815  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0
names       : X10, X20, X30, X40, X50, X60, X70, X80, X90, X100, layer
min values  :   0,   0,   0,   0,   0,   0,   0,   0,   0,    0,     0
max values  :   1,   1,   1,   1,   0,   1,   0,   1,   1,    0,     1

I want to do everything in memory, but now I need to use gdalwarp, and
don't know how to define the argument srcfile. I can write my stack to
somewhere using writeRaster and get a filename, but that is too slow for my
needs and I'm looking for an alternative.

Thank you
Hugo

        [[alternative HTML version deleted]]

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

Re: logLik value of the spml spatial panel model

Thu, 10/25/2018 - 03:27
On Thu, 25 Oct 2018, Jeremie Juste wrote:

> Hello,
>
> Indeed the Loglik was not returned here by the function sarpanelerror
> here. I have provided a patch (attached). It's an ugly one but better
> than nothing. So instead of running spml you could run test..spml after
> sourcing the file to obtain the fespaterr$logLik.
>
> fespaterr <- test..spml(fm, data = Produc, listw = mat2listw(usaww),
>                   model="within", spatial.error="b", Hess = FALSE)
>
>
> The only modification I've made is in the function test..(sperrorlm)
> where I've added ll=LL in the list returned.
Good, thanks! And good that Danlin reported that your addition of the
value to the returned object appeared to work for him. I've added the splm
maintainers to this reply. If you can produce a patch file using diff, it
would be easier to add the extra returned value. It may be, however, that
there was a reason for omitting it, so checks to verify that it actually
is a logLik value on the same basis as other logLik values (not from an
abbreviated log-likelihood function used for fitting by omitting
constants) would be helpful. It would be good to provide logLik methods
for the ML fitted model objects.

Again thanks for a good example of collaboration to benifit us all.

Roger

>
>
> Hope it helps,
>
> Jeremie
>
>
>
>

--
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]
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

Re: logLik value of the spml spatial panel model

Thu, 10/25/2018 - 00:39
Dear Jeremie:

That worked! Thank you very much for the addition.

Best,

Danlin


On 2018/10/25 1:27, Jeremie Juste wrote:
> Hello,
>
> Indeed the Loglik was not returned here by the function sarpanelerror
> here. I have provided a patch (attached). It's an ugly one but better
> than nothing. So instead of running spml you could run test..spml after
> sourcing the file to obtain the fespaterr$logLik.
>
> fespaterr <- test..spml(fm, data = Produc, listw = mat2listw(usaww),
>                     model="within", spatial.error="b", Hess = FALSE)
>
>
> The only modification I've made is in the function test..(sperrorlm)
> where I've added ll=LL in the list returned.
>
>
> Hope it helps,
>
> Jeremie
>
>
>
>
>
>
>> Dear Jeremie:
>>
>> Thank you for the response. I actually only copied the scripts from
>> the spml help page, reproduced here:
>>
>> fm <- log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp
>> ## the two standard specifications (SEM and SAR) one with FE
>> ## and the other with RE:
>> ## fixed effects panel with spatial errors
>> fespaterr <- spml(fm, data = Produc, listw = mat2listw(usaww),
>>                     model="within", spatial.error="b", Hess = FALSE) #
>> doesn't matter if I set Hess = T or F here, and doesn't matter if I
>> set spatial.error = "b" or "kkp"
>>
>> And when I want to check the logLik value:
>>
>> fespaterr$logLik
>>
>> I got the "NULL" value
>>
>> I checked this with some other data sets, and found that if I am using
>> the spatial lag specification, I usually got a value for the logLik
>> parameter. If, however, I use the spatial error specification, I
>> usually got the "NULL" value for the logLik parameter. I was not able
>> to get to the bottom of the issue by digging into the codes a bit. Any
>> leads would be greatly appreciated.
>>
>> Best,
>>
>> Danlin
>>
>> On 2018/10/24 15:30, Jeremie Juste wrote:
>>> Danlin Yu <[hidden email]> writes:
>>>
>>>
>>> Hello,
>>>
>>> Could you post the exact command you have used with all the relevant
>>> options? If you dug into the code then you've noticed quite a lot of
>>> nested control structures so a clear idea of your commands might be
>>> useful here.
>>>
>>> for instance
>>>
>>>       fm <- logc ~ logp+ logpn + logy
>>>       sarpool <-  suppressWarnings(splm::spml(fm,cigar,
>>>                                               listw=lwcig,
>>>                                               model="pooling",
>>>       spatial.error="none",lag=TRUE,index=c("region","time")))
>>>
>>> Will produce the logLik
>>>       sarpool$logLik
>>>
>>>
>>>   From what I've understood at some point the function <nlminb> is always
>>> called so a likelihood value is always produced although not always
>>> returned.
>>>
>>>
>>> Best regards,
>>>
>>> Jeremie
>>>
>>>
>>>
>>>> Dear List:
>>>>
>>>> In a recent attempt to compare spatial panel models using splm
>>>> package, I intend to use the logLik value that was produced by the
>>>> spml() routine. When I ran the model, I realize that the logLik value
>>>> is not always produced, especially when the spatial.error is set to
>>>> "b", a NULL value is returned. I know I must have missed something
>>>> here, but can anyone point out why the logLik value is not returned
>>>> when the spatial.error value is set to "b" (a spatial panel lag model
>>>> always return a logLik value). I tried a bit to dig into the codes,
>>>> but didn't find much that can answer my question.
>>>>
>>>> Thank you.
>>>>
>>>> Best,
>>>>
>>>> Danlin
--
___________________________________________
Danlin Yu, Ph.D.
Professor of GIS and Urban Geography
Department of Earth & Environmental Studies
Montclair State University
Montclair, NJ, 07043
Tel: 973-655-4313
Fax: 973-655-4072
Office: CELS 314
Email: [hidden email]
webpage: csam.montclair.edu/~yu


        [[alternative HTML version deleted]]

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

Re: logLik value of the spml spatial panel model

Thu, 10/25/2018 - 00:27
Hello,

Indeed the Loglik was not returned here by the function sarpanelerror
here. I have provided a patch (attached). It's an ugly one but better
than nothing. So instead of running spml you could run test..spml after
sourcing the file to obtain the fespaterr$logLik.

fespaterr <- test..spml(fm, data = Produc, listw = mat2listw(usaww),
                   model="within", spatial.error="b", Hess = FALSE)


The only modification I've made is in the function test..(sperrorlm)
where I've added ll=LL in the list returned.


Hope it helps,

Jeremie






> Dear Jeremie:
>
> Thank you for the response. I actually only copied the scripts from
> the spml help page, reproduced here:
>
> fm <- log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp
> ## the two standard specifications (SEM and SAR) one with FE
> ## and the other with RE:
> ## fixed effects panel with spatial errors
> fespaterr <- spml(fm, data = Produc, listw = mat2listw(usaww),
>                    model="within", spatial.error="b", Hess = FALSE) #
> doesn't matter if I set Hess = T or F here, and doesn't matter if I
> set spatial.error = "b" or "kkp"
>
> And when I want to check the logLik value:
>
> fespaterr$logLik
>
> I got the "NULL" value
>
> I checked this with some other data sets, and found that if I am using
> the spatial lag specification, I usually got a value for the logLik
> parameter. If, however, I use the spatial error specification, I
> usually got the "NULL" value for the logLik parameter. I was not able
> to get to the bottom of the issue by digging into the codes a bit. Any
> leads would be greatly appreciated.
>
> Best,
>
> Danlin
>
> On 2018/10/24 15:30, Jeremie Juste wrote:
>> Danlin Yu <[hidden email]> writes:
>>
>>
>> Hello,
>>
>> Could you post the exact command you have used with all the relevant
>> options? If you dug into the code then you've noticed quite a lot of
>> nested control structures so a clear idea of your commands might be
>> useful here.
>>
>> for instance
>>
>>      fm <- logc ~ logp+ logpn + logy
>>      sarpool <-  suppressWarnings(splm::spml(fm,cigar,
>>                                              listw=lwcig,
>>                                              model="pooling",
>>      spatial.error="none",lag=TRUE,index=c("region","time")))
>>
>> Will produce the logLik
>>      sarpool$logLik
>>
>>
>>  From what I've understood at some point the function <nlminb> is always
>> called so a likelihood value is always produced although not always
>> returned.
>>
>>
>> Best regards,
>>
>> Jeremie
>>
>>
>>
>>> Dear List:
>>>
>>> In a recent attempt to compare spatial panel models using splm
>>> package, I intend to use the logLik value that was produced by the
>>> spml() routine. When I ran the model, I realize that the logLik value
>>> is not always produced, especially when the spatial.error is set to
>>> "b", a NULL value is returned. I know I must have missed something
>>> here, but can anyone point out why the logLik value is not returned
>>> when the spatial.error value is set to "b" (a spatial panel lag model
>>> always return a logLik value). I tried a bit to dig into the codes,
>>> but didn't find much that can answer my question.
>>>
>>> Thank you.
>>>
>>> Best,
>>>
>>> Danlin
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

test.R (44K) Download Attachment

Pages