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

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

Re: logLik value of the spml spatial panel model

Wed, 10/24/2018 - 23:20
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

_______________________________________________
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

Wed, 10/24/2018 - 14:30
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

Re: Sf find number of parts in multigeometry

Tue, 10/23/2018 - 01:59
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]> wrote:

> I am looking for an equivalent to Postgis ST_NumGeometries
> 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]
> 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: Sf find number of parts in multigeometry

Tue, 10/23/2018 - 00:21
As an addition to obrl soil's answer, I have come to like

lengths(st_geometry(x)) - mind the s!

Best
Tim


On 10/23/2018 06:39 AM, obrl soil wrote:
> Hi Martin,
>
> for a multipolygon geometry column, you're working with a list of
> lists of lists of matrices, so you can iterate over the geometry
> column like
>
>      sapply(st_geometry(x), length)
>
> (or purrr::map_int(), perhaps) to report number of polygon components.
> To subset the first polygon component and drop the rest,
>
>      st_geometry(x) <- st_sfc(lapply(st_geometry(x), function(y)
> st_polygon(y[[1]])), crs = 4326) # or whatever crs you're using
>
> I would first run st_buffer(x, dist = 0L) and/or
> lwgeom::st_make_valid() to ensure the geometries are properly
> structured, although this can't solve every problem. Note that AFAIK
> there's no inherent sorting to multipolygon components, although you
> could probably reorder them by area or number of vertices before
> subsetting if you were really keen.
>
> cheers
> L
> On Tue, Oct 23, 2018 at 12:32 PM Martin Tomko <[hidden email]> wrote:
>> I am looking for an equivalent to Postgis ST_NumGeometries 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]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Sf find number of parts in multigeometry

Mon, 10/22/2018 - 23:39
Hi Martin,

for a multipolygon geometry column, you're working with a list of
lists of lists of matrices, so you can iterate over the geometry
column like

    sapply(st_geometry(x), length)

(or purrr::map_int(), perhaps) to report number of polygon components.
To subset the first polygon component and drop the rest,

    st_geometry(x) <- st_sfc(lapply(st_geometry(x), function(y)
st_polygon(y[[1]])), crs = 4326) # or whatever crs you're using

I would first run st_buffer(x, dist = 0L) and/or
lwgeom::st_make_valid() to ensure the geometries are properly
structured, although this can't solve every problem. Note that AFAIK
there's no inherent sorting to multipolygon components, although you
could probably reorder them by area or number of vertices before
subsetting if you were really keen.

cheers
L
On Tue, Oct 23, 2018 at 12:32 PM Martin Tomko <[hidden email]> wrote:
>
> I am looking for an equivalent to Postgis ST_NumGeometries 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]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Sf find number of parts in multigeometry

Mon, 10/22/2018 - 21:32
I am looking for an equivalent to Postgis ST_NumGeometries 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]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

logLik value of the spml spatial panel model

Mon, 10/22/2018 - 10:37
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
email: [hidden email]
webpage: csam.montclair.edu/~yu

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

Pages