|
|
相关文章推荐
I am thinking of applying spatial filetring to panel data model. Therefore
I would like to ask if you might know of some packages/codes or any other
resources for panel data that I could apply for my analysis.
It might be worth looking at the spmoran package (not panel, but under
active development). You can also look at adespatial for PCNM. In general,
the eigenvectors for the time-invariant neighbour object are also
time-invariant, unless your panel model is dynamic. Also find examples in
the literature and see whether they provide references or links to
software, for example but without direct software references:
https://research.vu.nl/ws/portalfiles/portal/73350412/06049,
https://doi.org/10.1177%2F0160017610386482.
https://scholar.google.no/scholar?oi=bibs=en=15767057313182942537_sdt=5
gives citations of Patuelli et al. (2011).
Hope this helps,
Roger
Best,
Martin Hulenyi
[[alternative HTML version deleted]]
R-sig-Geo mailing list
[email protected]
https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geodata=04%7C01%7CRoger.Bivand%40nhh.no%7C6c9307111ec74c16969508d9ced2d73b%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C637768228115992524%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C2000sdata=jqOHHKrl6SxZ441roWZNQKU5hTtS76nZFJ%2FEL56xGRI%3Dreserved=0
Roger Bivand
Emeritus Professor
Department of Economics, Norwegian School of Economics,
Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
e-mail: [email protected]
https://orcid.org/-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0J=en___
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
see if gdal understands the coordinate arrays already, you can just churn
it through the warper, otherwise construct a vrt that specifies them, which
is a bit obscure how you do it admittedly - but this is a gdal not an R
topic.
On Mon, 15 Nov 2021, 00:22 Edzer Pebesma,
wrote:
> On 14/11/2021 13:59, Edzer Pebesma wrote:
> > Gabriel shared the data with me (off-list, 1Gb) and it contains two hdf5
> > datasets, one containing the raster of the longitude and latitude, the
> > other rasters with the attributres (LST) but no coordinates. I could get
> > a (rough: factor 50 downsampling) plot with stars using:
> > library(stars)
> > lat =
> read_stars("HDF5:ECOSTRESS_L1B_GEO_11712_012_20200730T061029_0601_01.h5://Geolocation/latitude")
> > lon =
> read_stars("HDF5:ECOSTRESS_L1B_GEO_11712_012_20200730T061029_0601_01.h5://Geolocation/longitude")
> > var =
> > "HDF5:ECOSTRESS_L2_LSTE_11712_012_20200730T061029_0601_01.h5://SDS/LST"
> > lst = read_stars(var, curvilinear = list(x = lon, y = lat))
> > plot(lst, downsample = 50, axes = TRUE, reset = FALSE)
> > library(rnaturalearth)
> > ne = ne_countries(returnclass = "sf")
> > plot(st_geometry(ne), add = TRUE, border = 'yellow')
> > # california border seems to make sense
> > but I'm not sure how to make a regular raster (in some CRS) out of this
> > for exporting to GeoTIFF.
> This might be useful: https://github.com/r-spatial/stars/issues/386
> but converting to polygons will be incredibly slow and memory hungry.
> > On 14/11/2021 11:00, Edzer Pebesma wrote:
> >> This might be because the grid is not regular but possibly
> >> curvilinear, having two raster layers with the lon & lat values for
> >> each pixel. Can you share a sample data set?
> >> On 14/11/2021 10:47, Gabriel Cotlier wrote:
> >>> Dear Michael,
> >>> Following your advice I have tired :
> >>> library(terra)
> >>> library(raster)
> >>> library(dplyr)
> >>> r = terra::rast(file.choose())
> >>> plot(r$LST)
> >>> r_lst = r$LST
> >>> r_lst %>% raster()
> >>> plot(r_lst)
> >>> Then I effectively get a raster layer format from the raster package
> >>> wanted, but looking at the plot and the raster layer information--see
> >>> bellow-- the scene appears to need both latitude and longitude for
> >>> correct
> >>> geolocation, how can I do this procedure?
> r_lst %>% raster()
> >>> class : RasterLayer
> >>> dimensions : 5632, 5400, 30412800 (nrow, ncol, ncell)
> >>> resolution : 1, 1 (x, y)
> >>> extent : 0, 5400, 0, 5632 (xmin, xmax, ymin, ymax)
> >>> crs: NA
> >>> source : LST
> >>> names : LST
> >>> values : 0, 65535 (min, max)
> >>> The extent and the resolution are not the correct ones.
> >>> Maybe in the same .h5 is all the data for georeferencing it although I
> >>> could not see how to get to it and apply it to the layer.
> >>> There is also another file with complementary data related to
> >>> geolocation
> >>> but how can I query it, and get the data for geolocation and apply it
> >>> the raster layer?.
> >>> The objective is to save it as a GeoTIFF file to be opened in QGIS with
> >>> correct geolocation.
> >>> What do you suggest for converting the .h5 file to a georeferenced
> >>> Raster
> >>> Layer in GeoTIFF format ?
> >>> Thanks a lot.
> >>> Kind regards,
> >>> Gabriel
> >>> On Sat, Nov 13, 2021 at 12:22 PM Gabriel Cotlier
> >>> wrote:
> Dear Michael Summer,
> Thanks a lot for your swift reply.
> I will give a try to the procedure you mentioned.
> Kind regards,
> Gabriel
> On Sat, Nov 13, 2021 at 12:14 PM Michael Sumner
> wrote:
> > try terra::rast() or stars::read_stars() on the file, both use gdal
> > interrogate and read. What you get depends on the structure of the
> > files,
> > variables interpreted as subdatasets. stars is more general but sees
> > variables as 2D arrays with bands, terra calls these bands layers.
> > Alternatively try raster::raster() or RNetCDF or tidync which might
> > more straightforward than rhdf5 but all are quite different to each
> > other
> > and to the gdal view.
> > On Sat, 13 Nov 2021, 18:56 Gabriel Cotlier,
> > wrote:
> >> Hello,
> >> I would like to be able to visualize / plot a specific data layer
> >> within
> >> the structure a hierarchical data file format .h5 in R. I could
> >> observe
> >> the structure as follows :
> >> library(rhdf5)
> >> r= h5ls(file.choose())
> >> and also in this way identify the Group and Name of the data to be
> >> plotted,
> >> however I've tried I could not access the data itself, nor plot it
> >> raster layer with its corresponding geographic coordinates.
> >> I would like to be able to
for help with that.) But
> here are some hints:
> flipped color bar - I love to use RSeek.org which is an R-centric search
> engine - you may want to bookmark it or add it to your search engine list
> if your browser supports that. https://rseek.org/?q=leaflet+flip+color+bar
> This search hits a number of discussions on flipping the color bar.
> color intensity - oh, there is so much discussion in the world about
> colors. By default raster images have opacity of 1 (fully opaque) when
> rendered in leaflet, so really it is a matter of finding the color table
> you want and then applying colorNumeric judiciously. For me it comes down
> to trial and error. If you are looking for eye pop! then perhaps check out
> the viridis options. See ?colorNumeric
> Cheers,
> On Thu, Aug 12, 2021 at 2:24 PM Gabriel Cotlier
> wrote:
>> Dear Ben,
>> Thanks a lot for your help!
>> Actually it worked very well for me with the link you gave me.
>> However, for some reason the numbers of the values in color scale goes
>> from the lower values set at the top (blue) of the sacale to the
>> higher values set at the bottom (red ) and it would be better for me if it
>> goes from from the lower values in the bottom (blue) to higher value in the
>> top (red), since is temperature going from the lower to the higher.
>> Another issue I found a bit problematic for me to modify is that the
>> color palette is too "light", for instance not like the color palette "Jet"
>> in Matlab or Python, it seems as if the intensity of the colors of the
>> temperature colors is low or light maybe is a transparency and should
>> change opacity ?
>> Maybe you or somebody knows a possible way to improve these two issues a
>> bit ?
>> Here is the code :
>> # color palettes
>> # pal > "transparent" , reverse = TRUE)
>> pal > "transparent" , reverse = TRUE)
>> # plot map
>> leaflet() %>%
>> addTiles(urlTemplate = "
>> https://mts1.google.com/vt/lyrs=s=en=app={x}={y}={z}=G;,
>> attribution = 'Google') %>%
>> addPolygons(data = polygon,weight=5,col = 'black') %>%
>> addCircles(data = points, color = "White", radius = 500, fillOpacity =
>> 1,opacity = 9)%>%
>> addRasterImage(raster, project = FALSE, colors = pal) %>%
>> addLegend(pal = pal, values = values(raster), title = "Temperature",
>> opacity = 9) # ,labFormat = labelFormat(transform = function(x) sort(x,
>> decreasing = TRUE)))
>> Thanks a lot again for your help.
>> Kind regards,
>> Gabriel
>> On Thu, Aug 12, 2021 at 5:07 PM Ben Tupper wrote:
>>> See the "markers" and "raster images" sections here
>>> http://rstudio.github.io/leaflet/markers.html
>>> Cheers,
>>> On Thu, Aug 12, 2021 at 5:02 AM Gabriel Cotlier
>>> wrote:
Hello.
I would like to use Leaflet package to plot over a Google Satellite base
map :
1. a shapefile of polygon
2. a shapefile of points
3, a GeoTIFF image
I could use the Leaflet package to get plotted successfully only the
first
Item of the list above with the following code :
require(rgdal)
library(rgeos)
library(raster)
shapeData %
addTiles(urlTemplate = "
https://mts1.google.com/vt/lyrs=s=en=app={x}={y}={z}=G;,
attribution = 'Google') %>%
addPolygons(data=shapeData,weight=5,col = 'red')
How is it possible to find a way to complete the code above for plotting
the item 2 and 3 as well in the same Leaflet figure.
Is there any possible solution that can enable this task?
Thanks for your help.
Kind regards
Gabriel
[[alternative HTML version deleted]]
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>> Ben Tupper (he/him)
>>> Bigelow
for help with that.) But here
are some hints:
flipped color bar - I love to use RSeek.org which is an R-centric search
engine - you may want to bookmark it or add it to your search engine list
if your browser supports that. https://rseek.org/?q=leaflet+flip+color+bar
This search hits a number of discussions on flipping the color bar.
color intensity - oh, there is so much discussion in the world about
colors. By default raster images have opacity of 1 (fully opaque) when
rendered in leaflet, so really it is a matter of finding the color table
you want and then applying colorNumeric judiciously. For me it comes down
to trial and error. If you are looking for eye pop! then perhaps check out
the viridis options. See ?colorNumeric
Cheers,
On Thu, Aug 12, 2021 at 2:24 PM Gabriel Cotlier wrote:
> Dear Ben,
> Thanks a lot for your help!
> Actually it worked very well for me with the link you gave me.
> However, for some reason the numbers of the values in color scale goes
> from the lower values set at the top (blue) of the sacale to the
> higher values set at the bottom (red ) and it would be better for me if it
> goes from from the lower values in the bottom (blue) to higher value in the
> top (red), since is temperature going from the lower to the higher.
> Another issue I found a bit problematic for me to modify is that the
> color palette is too "light", for instance not like the color palette "Jet"
> in Matlab or Python, it seems as if the intensity of the colors of the
> temperature colors is low or light maybe is a transparency and should
> change opacity ?
> Maybe you or somebody knows a possible way to improve these two issues a
> bit ?
> Here is the code :
> # color palettes
> # pal "transparent" , reverse = TRUE)
> pal "transparent" , reverse = TRUE)
> # plot map
> leaflet() %>%
> addTiles(urlTemplate = "
> https://mts1.google.com/vt/lyrs=s=en=app={x}={y}={z}=G;,
> attribution = 'Google') %>%
> addPolygons(data = polygon,weight=5,col = 'black') %>%
> addCircles(data = points, color = "White", radius = 500, fillOpacity =
> 1,opacity = 9)%>%
> addRasterImage(raster, project = FALSE, colors = pal) %>%
> addLegend(pal = pal, values = values(raster), title = "Temperature",
> opacity = 9) # ,labFormat = labelFormat(transform = function(x) sort(x,
> decreasing = TRUE)))
> Thanks a lot again for your help.
> Kind regards,
> Gabriel
> On Thu, Aug 12, 2021 at 5:07 PM Ben Tupper wrote:
>> See the "markers" and "raster images" sections here
>> http://rstudio.github.io/leaflet/markers.html
>> Cheers,
>> On Thu, Aug 12, 2021 at 5:02 AM Gabriel Cotlier
>> wrote:
>>> Hello.
>>> I would like to use Leaflet package to plot over a Google Satellite base
>>> map :
>>> 1. a shapefile of polygon
>>> 2. a shapefile of points
>>> 3, a GeoTIFF image
>>> I could use the Leaflet package to get plotted successfully only the
>>> first
>>> Item of the list above with the following code :
>>> require(rgdal)
>>> library(rgeos)
>>> library(raster)
>>> shapeData >> shapeData >>
>>> leaflet() %>%
>>> addTiles(urlTemplate = "
>>> https://mts1.google.com/vt/lyrs=s=en=app={x}={y}={z}=G;,
>>> attribution = 'Google') %>%
>>> addPolygons(data=shapeData,weight=5,col = 'red')
>>> How is it possible to find a way to complete the code above for plotting
>>> the item 2 and 3 as well in the same Leaflet figure.
>>> Is there any possible solution that can enable this task?
>>> Thanks for your help.
>>> Kind regards
>>> Gabriel
>>> [[alternative HTML version deleted]]
>>> R-sig-Geo mailing list
>>> [email protected]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> Ben Tupper (he/him)
>> Bigelow Laboratory for Ocean Science
>> East Boothbay, Maine
>> http://www.bigelow.org/
>> https://eco.bigelow.org
Ben Tupper (he/him)
Bigelow Laboratory for Ocean Science
East Boothbay, Maine
http://www.bigelow.org/
https://eco.bigelow.org
[[alternative HTML version deleted]]
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Thanks a lot for your help!
Actually it worked very well for me with the link you gave me.
However, for some reason the numbers of the values in color scale goes from
the lower values set at the top (blue) of the sacale to the higher values
set at the bottom (red ) and it would be better for me if it goes from from
the lower values in the bottom (blue) to higher value in the top (red),
since is temperature going from the lower to the higher.
Another issue I found a bit problematic for me to modify is that the color
palette is too "light", for instance not like the color palette "Jet" in
Matlab or Python, it seems as if the intensity of the colors of the
temperature colors is low or light maybe is a transparency and should
change opacity ?
Maybe you or somebody knows a possible way to improve these two issues a
bit ?
Here is the code :
# color palettes
# pal %
addTiles(urlTemplate = "
https://mts1.google.com/vt/lyrs=s=en=app={x}={y}={z}=G;,
attribution = 'Google') %>%
addPolygons(data = polygon,weight=5,col = 'black') %>%
addCircles(data = points, color = "White", radius = 500, fillOpacity =
1,opacity = 9)%>%
addRasterImage(raster, project = FALSE, colors = pal) %>%
addLegend(pal = pal, values = values(raster), title = "Temperature",
opacity = 9) # ,labFormat = labelFormat(transform = function(x) sort(x,
decreasing = TRUE)))
Thanks a lot again for your help.
Kind regards,
Gabriel
On Thu, Aug 12, 2021 at 5:07 PM Ben Tupper wrote:
> See the "markers" and "raster images" sections here
> http://rstudio.github.io/leaflet/markers.html
> Cheers,
> On Thu, Aug 12, 2021 at 5:02 AM Gabriel Cotlier
> wrote:
>> Hello.
>> I would like to use Leaflet package to plot over a Google Satellite base
>> map :
>> 1. a shapefile of polygon
>> 2. a shapefile of points
>> 3, a GeoTIFF image
>> I could use the Leaflet package to get plotted successfully only the first
>> Item of the list above with the following code :
>> require(rgdal)
>> library(rgeos)
>> library(raster)
>> shapeData > shapeData >
>> leaflet() %>%
>> addTiles(urlTemplate = "
>> https://mts1.google.com/vt/lyrs=s=en=app={x}={y}={z}=G;,
>> attribution = 'Google') %>%
>> addPolygons(data=shapeData,weight=5,col = 'red')
>> How is it possible to find a way to complete the code above for plotting
>> the item 2 and 3 as well in the same Leaflet figure.
>> Is there any possible solution that can enable this task?
>> Thanks for your help.
>> Kind regards
>> Gabriel
>> [[alternative HTML version deleted]]
>> R-sig-Geo mailing list
>> [email protected]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> Ben Tupper (he/him)
> Bigelow Laboratory for Ocean Science
> East Boothbay, Maine
> http://www.bigelow.org/
> https://eco.bigelow.org
[[alternative HTML version deleted]]
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
See the "markers" and "raster images" sections here
http://rstudio.github.io/leaflet/markers.html
Cheers,
On Thu, Aug 12, 2021 at 5:02 AM Gabriel Cotlier wrote:
> Hello.
> I would like to use Leaflet package to plot over a Google Satellite base
> map :
> 1. a shapefile of polygon
> 2. a shapefile of points
> 3, a GeoTIFF image
> I could use the Leaflet package to get plotted successfully only the first
> Item of the list above with the following code :
> require(rgdal)
> library(rgeos)
> library(raster)
> shapeData shapeData
> leaflet() %>%
> addTiles(urlTemplate = "
> https://mts1.google.com/vt/lyrs=s=en=app={x}={y}={z}=G;,
> attribution = 'Google') %>%
> addPolygons(data=shapeData,weight=5,col = 'red')
> How is it possible to find a way to complete the code above for plotting
> the item 2 and 3 as well in the same Leaflet figure.
> Is there any possible solution that can enable this task?
> Thanks for your help.
> Kind regards
> Gabriel
> [[alternative HTML version deleted]]
> R-sig-Geo mailing list
> [email protected]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Ben Tupper (he/him)
Bigelow Laboratory for Ocean Science
East Boothbay, Maine
http://www.bigelow.org/
https://eco.bigelow.org
[[alternative HTML version deleted]]
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> Maybe consider st_coordinates() to get a matrix of coordinates,
> preferably from st_geometry() of the "sf" object, as hard-coding the name
> of the column containing the "sfc" object might be affected in some cases
> by a different name. Then subset the matrix as you wish.
> Hope this helps,
> Roger
> On Sun, 1 Aug 2021, Erin Hodgess wrote:
> > Hello!
> > I have the following sf object:
> >> bt.cent
> > Simple feature collection with 8 features and 1 field
> > Geometry type: POINT
> > Dimension: XY
> > Bounding box: xmin: -105.0746 ymin: 38.97641 xmax: -104.7917 ymax:
> 39.82494
> > CRS: +proj=longlat +datum=WGS84 +no_defs
> > rnorm.10. geometry
> > 1 0.3297790 POINT (-105.0746 38.97641)
> > 2 NA POINT (-104.7917 38.97641)
> > 3 -1.0734667 POINT (-105.0746 39.25925)
> > 4 NA POINT (-104.7917 39.25925)
> > 5 -0.1284898 POINT (-105.0746 39.5421)
> > 6 -0.4071287 POINT (-104.7917 39.5421)
> > 7 NA POINT (-105.0746 39.82494)
> > 8 1.4171128 POINT (-104.7917 39.82494)
> > And I am accessing the geometry (to build covariances) as:
> >> bt.cent$geometry[[1]][1]
> > [1] -105.0746
> >> bt.cent$geometry[[1]][2]
> > [1] 38.97641
> > Is that the best way to handle the geometry items, please?
> > Thanks in advance,
> > Sincerely,
> > Erin Hodgess, PhD
> > mailto: [email protected]
> > [[alternative HTML version deleted]]
> > R-sig-Geo mailing list
> > [email protected]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> Roger Bivand
> Emeritus Professor
> Department of Economics, Norwegian School of Economics,
> Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
> e-mail: [email protected]
> https://orcid.org/-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0J=en
Erin Hodgess, PhD
mailto: [email protected]
[[alternative HTML version deleted]]
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Maybe consider st_coordinates() to get a matrix of coordinates,
preferably from st_geometry() of the "sf" object, as hard-coding the name
of the column containing the "sfc" object might be affected in some cases
by a different name. Then subset the matrix as you wish.
Hope this helps,
Roger
On Sun, 1 Aug 2021, Erin Hodgess wrote:
Hello!
I have the following sf object:
bt.cent
Simple feature collection with 8 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -105.0746 ymin: 38.97641 xmax: -104.7917 ymax: 39.82494
CRS: +proj=longlat +datum=WGS84 +no_defs
rnorm.10. geometry
1 0.3297790 POINT (-105.0746 38.97641)
2 NA POINT (-104.7917 38.97641)
3 -1.0734667 POINT (-105.0746 39.25925)
4 NA POINT (-104.7917 39.25925)
5 -0.1284898 POINT (-105.0746 39.5421)
6 -0.4071287 POINT (-104.7917 39.5421)
7 NA POINT (-105.0746 39.82494)
8 1.4171128 POINT (-104.7917 39.82494)
And I am accessing the geometry (to build covariances) as:
bt.cent$geometry[[1]][1]
[1] -105.0746
bt.cent$geometry[[1]][2]
[1] 38.97641
Emeritus Professor
Department of Economics, Norwegian School of Economics,
Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
e-mail: [email protected]
https://orcid.org/-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0J=en
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
It suggests you do not have the required packages for those methods on your
machine. You can first simply run the following code and then try yours:
sdm::installAll()
Hope this helps,
Best wishes,
Babak
On Tue, Jun 29, 2021 at 1:00 PM wrote:
> Date: Mon, 28 Jun 2021 19:18:01 +0100
> From: Francisco Borges
> To: [email protected]
> Subject: [R-sig-Geo] Question concerning the SDM package in R
> Message-ID:
> ca+btngownxrvp_o47_a5p1u_2zcjhsok1dv8qzemdhi2hf7...@mail.gmail.com>
> Content-Type: text/plain; charset="utf-8"
> Hello all!
> I have encountered an issue with the SDM package by Babak Naimi and Miguel
> Araújo in R. Specifically, when re-running a model which I had successfully
> run a couple of months ago, I get the following error:
> «Error in sdmSetting(test.percent = 30, replication = c("sub","boot", :
> methods rf and brt do not exist!»
> This is the code I'm inputting:
> model sdm(
> occ ~.,
> data = dalter,
> methods = c('glm','rf','brt'),
> test.p = 30,
> replication = c("sub", "boot", "cv"),
> n.folds = 4,
> n = 5
> This is weird. It is stating that the sdm package does not have boosted
> regression tree nor random forest as an available method. I have
> reinstalled the package through CRAN in R today.
> Any suggestions?
> Thank you in advance,
> *Francisco Borges*
> MARE – Marine and Environmental Sciences Centre
> Laboratório Marítimo da Guia - Faculdade de Ciências da Universidade de
> Lisboa
> Av. Nossa Senhora do Cabo, 939
> 2750-374 Cascais, Portugal
> Tel: +351 214 869 211
> orcid.org/-0002-3911-0421
> http://www.mare-centre.pt/pt/user/8232
> https://www.cienciavitae.pt//pt/B814-9E8E-F556
> [[alternative HTML version deleted]]
> R-sig-Geo mailing list
> [email protected]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> End of R-sig-Geo Digest, Vol 214, Issue 14
[[alternative HTML version deleted]]
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> The general syntax is:
> writeRaster(NameRaster, file="PathOfFolderSTinDesktop/NameFile.tif",
> format= ""GTiff", overwrite=TRUE)
> If you want to add options with bands:
> writeRaster(NameRaster, file="FolderPath/NameFile.tif", format= ""GTiff",
> overwrite=TRUE, options= c("INTERLEAVE=BAND","COMPRESS=LZW"))
> NameRaster is a RasterLayer or RasterBrick t in your global environment.
> Hope this works,
> Yasna
>> I'm trying to save a raster on workspace as GeoTIFF file into a folder
>> named ST in Desktop
>> writeRaster(x = ST , file = paste(n,"_Rescaled_band",sep=" "),
>> format="GTiff", overwrite = T)
>> and got the following error message :
>> Error in .local(x, filename, ...) :
>> Attempting to write a file to a path that does not exist:
>> Thanks for your help.
>> Regards.
>> [[alternative HTML version deleted]]
>> R-sig-Geo mailing list
>> [email protected]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> *Yasna K Palmeiro Silva.*
> *Tildes han sido omitidos intencionalmente
[[alternative HTML version deleted]]
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
The general syntax is:
writeRaster(NameRaster, file="PathOfFolderSTinDesktop/NameFile.tif",
format= ""GTiff", overwrite=TRUE)
If you want to add options with bands:
writeRaster(NameRaster, file="FolderPath/NameFile.tif", format= ""GTiff",
overwrite=TRUE, options= c("INTERLEAVE=BAND","COMPRESS=LZW"))
NameRaster is a RasterLayer or RasterBrick t in your global environment.
Hope this works,
Yasna
> I'm trying to save a raster on workspace as GeoTIFF file into a folder
> named ST in Desktop
> writeRaster(x = ST , file = paste(n,"_Rescaled_band",sep=" "),
> format="GTiff", overwrite = T)
> and got the following error message :
> Error in .local(x, filename, ...) :
> Attempting to write a file to a path that does not exist:
> Thanks for your help.
> Regards.
> [[alternative HTML version deleted]]
> R-sig-Geo mailing list
> [email protected]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
*Yasna K Palmeiro Silva.*
*Tildes han sido omitidos intencionalmente
[[alternative HTML version deleted]]
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Shape files are notoriously bad, and you just found one reason why.
I would suggest to look int e.g. GRASS (https://grass.osgeo.org
) or QGIS (https://www.qgis.org/en/site/) for doing
the cleaning and joining. Afterwards, you can go back to R to do further
analysis.
Cheers,
Rainer
> On 18 Mar 2021, at 14:02, Sean Trende wrote:
> Hello,
> First, thank you for allowing me to join this group. To be able to (quietly)
> observe discussions from so many great minds in an area so important for my
> day-to-day job is something I am greatly looking forward to.
> I do have a question that I hope someone can help with. For my dissertation,
> I am building an areal model off of Ohio precincts. The shapefile
> (https://github.com/mggg/ohio-precincts) is notoriously filled with slivers,
> which complicates the adjacency matrix.
> I've come up with a wonky solution (code at bottom) where, after loading the
> shapefile and filtering by county I manually increase the size of the snap by
> .01 and look to see whether the number of nonzero links increases. If it
> does, I can then plot the new and old adjacency matrices, and eyeball the new
> link(s) to determine whether it is "real" or if the snap has grown so large I
> am linking precincts that are not actually adjacent. I am hoping I can find
> a snap level that links most of the truly adjacent precincts, without linking
> too many precincts that are not actually adjacent. I've gotten through ten
> counties, and it looking like the best snap is ~30.
> This is obviously time-consuming, and in cities it can be tough to see
> whether precincts are actually adjacent given their small size. My questions
> 1) Is there a way to extract the number nonzero links from an 'nb' class of
> files and assign it to a variable? This would allow me to automate the
> inquiry and examine only those snaps that increased the number of links,
> saving time.
> 2) Am I missing a much more straightforward solution to cleaning the
> shapefile entirely?
> Thank you for your time,
> df vect j
> selected_county df_county
> ia ib
> nb_a nb_b coords
> nb_b #This allows you to look and see if there are more links with the
> larger snap than with the smaller snap
> plot(df_county)
> plot(nb_b, coords, col = "blue", add = T, lwd = 2)
> plot(nb_a, coords, col = "grey", add = T, lwd = 2) #Because of the
> overlay, the new connection will be blue
> ia ib the comparison by a similar amount
> #so if this is the first time you run it, your next comparison will be
> snaps of 0.01 and 0.02
> #To keep it from resetting I just highlight up to the assignment of nb_a,
> and then hit cntrl_enter until the number of links increases
> R-sig-Geo mailing list
> [email protected]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Hello Roger, thanks for the tips.In fact, I believe that I will extract the coordinates by the gCentroid function, from the rgeos package.The results look better.Also, thank you very much for the GWR tip, in my case, it is just a complementary result, but in fact, is good to know about it.A hug.Your tips are helping me so much!Em 4 de nov de 2020 17:32, Roger Bivand escreveu:On Wed, 4 Nov 2020, Pietro Andre Telatin Paschoalino wrote:
> Hello Roger,
> Yes, I'm using spdep package to construct the matrices.
Do not use matrices, use neighbour or weights objects. Dense matrices with
mostly zero values are a waste of space.
> I used the knearneigh function (in which I entered the coordinates), as
> well as the poly2nb function.
Yes, however, you change the support of the aggregated data units by doing
this, so knn is less well supported than contiguity in general.
> I think that I only have one of these rings but the size of the polygons
> is relatively large (they are political divisions) of a country. In that
> case, can I keep extracting coordinates by the coordinates function?
The centroids are more dependent (and so the distances between point
representations) on the boundary geometries than are simple contiguities.
Change of support is always something of a problem. Practically, yes, you
can do it, but whether it makes sense or not in terms of the data is a
different question.
> I didn't understand why you suggested avoiding GWR. Are you commenting
> on the methodology itself?
Yes, it may be used for exploring to detect possible omitted covariates
but should not be used for modelling unless that is definitely the data
generating process. Random data generates patterned (insignificant)
coefficient maps.
Roger
> Pietro Andre Telatin Paschoalino
> Doutorando em Ciências Econômicas da Universidade Estadual de Maringá - PCE.
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [email protected]
https://orcid.org/-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0J=en
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Do not use matrices, use neighbour or weights objects. Dense matrices with
mostly zero values are a waste of space.
I used the knearneigh function (in which I entered the coordinates), as
well as the poly2nb function.
Yes, however, you change the support of the aggregated data units by doing
this, so knn is less well supported than contiguity in general.
I think that I only have one of these rings but the size of the polygons
is relatively large (they are political divisions) of a country. In that
case, can I keep extracting coordinates by the coordinates function?
The centroids are more dependent (and so the distances between point
representations) on the boundary geometries than are simple contiguities.
Change of support is always something of a problem. Practically, yes, you
can do it, but whether it makes sense or not in terms of the data is a
different question.
I didn't understand why you suggested avoiding GWR. Are you commenting
on the methodology itself?
Yes, it may be used for exploring to detect possible omitted covariates
but should not be used for modelling unless that is definitely the data
generating process. Random data generates patterned (insignificant)
coefficient maps.
Roger
Pietro Andre Telatin Paschoalino
Doutorando em Ciências Econômicas da Universidade Estadual de Maringá - PCE.
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [email protected]
https://orcid.org/-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0J=en___
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Yes, I'm using spdep package to construct the matrices.
I used the knearneigh function (in which I entered the coordinates), as well as
the poly2nb function.
I think that I only have one of these rings but the size of the polygons is
relatively large (they are political divisions) of a country. In that case, can
I keep extracting coordinates by the coordinates function?
I didn't understand why you suggested avoiding GWR. Are you commenting on the
methodology itself?
Pietro Andre Telatin Paschoalino
Doutorando em Ci�ncias Econ�micas da Universidade Estadual de Maring� - PCE.
[[alternative HTML version deleted]]
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
What do you mean by hole? Non-contiguous regions? If so, I don't think
it's a problem, because I don't have any non-contiguous regions. I have
a blank polygon that is not computed because it is a lake (it is blank).
An interior ring within a polygon (like a lake, or another region
fully within another, which then has its own exterior ring).
I am using coordinates to build spatial weight matrices to use in
spatial regression / geographically weighted regressions.
Please do not try to do this by hand. You can use topological predicates
(in rgeos or sf), or simply use spdep::poly2nb(), which does this. Avoid
GWR whenever possible, please.
Roger
Pietro Andre Telatin Paschoalino
Doutorando em Ciências Econômicas da Universidade Estadual de Maringá - PCE.
De: Roger Bivand
Enviado: quarta-feira, 4 de novembro de 2020 10:38
Para: Pietro Andre Telatin Paschoalino
Cc: [email protected]
Assunto: Re: [R-sig-Geo] Question about function COORDINATES in SP PACKAGE
On Wed, 4 Nov 2020, Pietro Andre Telatin Paschoalino wrote:
Hello everyone.
Could anyone help with a small question about the coordinates in sp
package function?
I want to retrieve the geographic coordinates of polygons (shape file),
and reading the documentation, I did not find out if these coordinates
refer to the centroid of polygon.
So, in fact, do such coordinates refer to the centroid of each polygon?
Yes, they report the centroid of the largest polygon, but may not handle
holes correctly. Why do you need to retreive the complete set of
coordinates? That would be much easier using the sf package, but is seldom
needed.
Roger
Thanks.
Pietro Andre Telatin Paschoalino
Doutorando em Ci�ncias Econ�micas da Universidade Estadual de Maring� - PCE.
[[alternative HTML version deleted]]
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [email protected]
https://orcid.org/-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0J=en
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [email protected]
https://orcid.org/-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0J=en___
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
What do you mean by hole? Non-contiguous regions? If so, I don't think it's a
problem, because I don't have any non-contiguous regions. I have a blank
polygon that is not computed because it is a lake (it is blank).
I am using coordinates to build spatial weight matrices to use in spatial
regression / geographically weighted regressions.
Pietro Andre Telatin Paschoalino
Doutorando em Ciências Econômicas da Universidade Estadual de Maringá - PCE.
De: Roger Bivand
Enviado: quarta-feira, 4 de novembro de 2020 10:38
Para: Pietro Andre Telatin Paschoalino
Cc: [email protected]
Assunto: Re: [R-sig-Geo] Question about function COORDINATES in SP PACKAGE
On Wed, 4 Nov 2020, Pietro Andre Telatin Paschoalino wrote:
> Hello everyone.
> Could anyone help with a small question about the coordinates in sp
> package function?
> I want to retrieve the geographic coordinates of polygons (shape file),
> and reading the documentation, I did not find out if these coordinates
> refer to the centroid of polygon.
> So, in fact, do such coordinates refer to the centroid of each polygon?
Yes, they report the centroid of the largest polygon, but may not handle
holes correctly. Why do you need to retreive the complete set of
coordinates? That would be much easier using the sf package, but is seldom
needed.
Roger
> Thanks.
> Pietro Andre Telatin Paschoalino
> Doutorando em Ci�ncias Econ�micas da Universidade Estadual de Maring� - PCE.
>[[alternative HTML version deleted]]
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [email protected]
https://orcid.org/-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0J=en
[[alternative HTML version deleted]]
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
I want to retrieve the geographic coordinates of polygons (shape file),
and reading the documentation, I did not find out if these coordinates
refer to the centroid of polygon.
So, in fact, do such coordinates refer to the centroid of each polygon?
Yes, they report the centroid of the largest polygon, but may not handle
holes correctly. Why do you need to retreive the complete set of
coordinates? That would be much easier using the sf package, but is seldom
needed.
Roger
Thanks.
Pietro Andre Telatin Paschoalino
Doutorando em Ci�ncias Econ�micas da Universidade Estadual de Maring� - PCE.
[[alternative HTML version deleted]]
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [email protected]
https://orcid.org/-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0J=en___
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Dear Prof. Bivand, thank you very much for the super quick reply and fix.
The approach using *spatialreg* and as(x, "CsparseMatrix") works perfectly.
This step makes little sense - why is it included?
I included the Cholesky decomposition since I found this problem while I
was estimating the Cholesky decomposition of a precision matrix in a
multivariate CAR model (as you probably guessed) and the function *chol *was
complaining that D - alpha * W was not symmetric (where D is a diagonal
matrix with the total number of neighbours for each spatial unit and alpha
is the parameter controlling the strength in the proper CAR model). D is
clearly symmetric by construction so the problem was with W and I couldn't
understand why. I did not include all the details about how or why I found
the error because I think they were not important but, unfortunately, the
example with *chol(W)* is also meaningless.
Anyway thank you very much for the solution and the other suggestions.
Kind regards
Andrea
> On Wed, 29 Jul 2020, Andrea Gilardi wrote:
> > Dear all, I think the following message should work. Excuse me for the
> > previous message.
> > library(spdep)
> > library(rgdal)
> > library(Matrix)
> > nc_sids
> Do things in steps:
> adj is.symmetric.nb(adj, force=TRUE)
> # [1] TRUE
> W0 all(W0 == t(W0))
> # [1] TRUE
> > Matrix::isSymmetric(W)
> > chol(W)
> W isSymmetric(W)
> #[1] TRUE
> res #Error in asMethod(object) : not a positive definite matrix
> This step makes little sense - why is it included?
> I don't know why coercing to Matrix does not see the matrix as symmetric.
> lw library(spatialreg)
> W #> isSymmetric(W)
> #[1] TRUE
> This works and is what you might need:
> or see ?Matrix::USCounties.
> Roger
> > If it's still not working I added a reprex here:
> > https://gist.github.com/agila5/59d3a173df9b1efabe13800f748a2d48
> Not useful, not in thread nor do any answers stay in thread.
> > Kind regards
> > Andrea
> > Il giorno mer 29 lug 2020 alle ore 17:19 Roger Bivand <
> [email protected]>
> > ha scritto:
> >> Please follow-up in this thread posting in plain text (not HTML). As you
> >> can see below, HTML garbles your message.
> >> On Wed, 29 Jul 2020, Andrea Gilardi wrote:
> >>> Dear all, probably I'm missing something obvious but I have a question
> >>> the matrix built using nb2mat with argument style = "B". More
> precisely,
> >>> don't understand why the following code says that the matrix *W* (which
> >>> should be a binary first-order adjacency matrix) is not symmetric:
> >>> *library(spdep)library(rgdal)library(Matrix)# datanc.sids >>> readOGR(system.file("shapes/sids.shp", package = "spData"))adj >>> poly2nb(nc.sids)W >>> "Matrix")Matrix::isSymmetric(W)#> [1] FALSEchol(W)#> Error in
> >>> asMethod(object): not a symmetric matrix; consider forceSymmetric() or
> >>> symmpart()*
> >>> Kind regards
> >>> Andrea
> >>> [[alternative HTML version deleted]]
> >>> R-sig-Geo mailing list
> >>> [email protected]
> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >> Roger Bivand
> >> Department of Economics, Norwegian School of Economics,
> >> Helleveien 30, N-5045 Bergen, Norway.
> >> voice: +47 55 95 93 55; e-mail: [email protected]
> >> https://orcid.org/-0003-2392-6140
> >> https://scholar.google.no/citations?user=AWeghB0J=en
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: [email protected]
> https://orcid.org/-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0J=en
[[alternative HTML version deleted]]
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
library(Matrix)
nc_sids Please follow-up in this thread posting in plain text (not HTML). As you
> can see below, HTML garbles your message.
> On Wed, 29 Jul 2020, Andrea Gilardi wrote:
> > Dear all, probably I'm missing something obvious but I have a question on
> > the matrix built using nb2mat with argument style = "B". More precisely,
> > don't understand why the following code says that the matrix *W* (which
> > should be a binary first-order adjacency matrix) is not symmetric:
> > *library(spdep)library(rgdal)library(Matrix)# datanc.sids > readOGR(system.file("shapes/sids.shp", package = "spData"))adj > poly2nb(nc.sids)W > "Matrix")Matrix::isSymmetric(W)#> [1] FALSEchol(W)#> Error in
> > asMethod(object): not a symmetric matrix; consider forceSymmetric() or
> > symmpart()*
> > Kind regards
> > Andrea
> > [[alternative HTML version deleted]]
> > R-sig-Geo mailing list
> > [email protected]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: [email protected]
> https://orcid.org/-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0J=en
[[alternative HTML version deleted]]
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Please follow-up in this thread posting in plain text (not HTML). As you
can see below, HTML garbles your message.
On Wed, 29 Jul 2020, Andrea Gilardi wrote:
Dear all, probably I'm missing something obvious but I have a question on
the matrix built using nb2mat with argument style = "B". More precisely, I
don't understand why the following code says that the matrix *W* (which
should be a binary first-order adjacency matrix) is not symmetric:
*library(spdep)library(rgdal)library(Matrix)# datanc.sids [1] FALSEchol(W)#> Error in
asMethod(object): not a symmetric matrix; consider forceSymmetric() or
symmpart()*
Kind regards
Andrea
[[alternative HTML version deleted]]
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [email protected]
https://orcid.org/-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0J=en
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Dear users,
I would like to download OSM data and convert it to shapefile and change
coordinates to same as Landsat data.
I'm following the code bellow.
Is this the correct way for saving/writing to swapfile the OSM data?
Since you are using the sf package, why not use it for the other
operations? There is no need to involve sp or rgdal unless you want to for
some other reason. Maybe you are copying from an old blog or SO source?
Have you visited: https://www.r-spatial.org/ with plenty of posts on
migrating to sf?
With regard to coordinate operations, you (and everybody else) will need
to find out how to define CRS without using Proj.4 strings; this change is
coming soon! Yours looks like EPSG:32636. Both sf and sp will be moving
very soon to WKT2-2019 from Proj.4 strings (Proj.4 strings will still be
available, but datums and other specification components may simply drop
out without warning because of changes in PROJ and GDAL).
Finally, stop using shapefiles and change to GeoPackage straight away,
there is no reason to continue to use a format past its end-of-life.
Roger
Thanks a lot for your response.
Best,
library(osmar)
library(osmdata)
library(OpenStreetMap)
library(tmap)
# Paris administrative borders ID
opq_string () %>%
osmdata_sf ()
qtm(dat2$osm_multipolygons)
library(sf)
st_write(dat2$osm_multipolygons, "s.shp")
# change coordinate system to Landsat data
landst_crs= "+proj=utm +zone=36 +datum=WGS84 +units=m +no_defs +ellps=WGS84
+towgs84=0,0,0"
s=readOGR('/Users/gcotlier/Desktop/s.shp')
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55 ; e-mail:
> [email protected]
> https://orcid.org/-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0J=en
[[alternative HTML version deleted]]
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Thank you very much for your reply.
My thesis supervisor told me to use GWR to explore small-scale differences and
test the consistency of the SAR models.
So, I tried to fit the GWR model with all my observations. However, if I
understood you correctly, I would have to choose smaller sections of the grid
and try to fit a GWR there. What I am not sure about how to do this in order to
get informative results.
As for the SAR model, the dependent variable consists n blocks where the m-th
entry of each block corresponds to the m-th cell of my grid. So, if I have to
reduce the grid to about 5K entries, I would need to take a subset of roughly
250 cells and then take all the 20 observations for each of these cells? And
then try to fit a GWR there?
Best,
Raphael
> Am 17.05.2019 um 14:20 schrieb Roger Bivand :
> On Fri, 17 May 2019, Raphael Mesaric via R-sig-Geo wrote:
>> Dear all,
>> Is there an option to shorten the running time for the gwr() function,
>> similar to the ‚LU‘ method for lagsarlm() in the spdep package? Because I
>> have a model with roughly 500’000 observations, and the running time at the
>> moment is quite long, respectively it has not yet terminated.
> What are you actually doing? Why did you choose GWR? Are you fitting a GWR
> with 500K observations, or have you fitted a GWR with many fewer
> observations, and are now rendering that fitted model with 500K fit points?
> GWR is only for detecting possible non-stationarity or similar
> mis-specification in moderately sized data sets. Trying to fit with 500K
> gives a dense hat matrix of 500K x 500K, which is imposssible (or were it
> possible would be uninformative). Think of 5K as a sensible maximum if GWR is
> condidered sensible at all. I would think that finding a bandwidth is
> impossible too.
> Roger
>> Thank you for your help in advance.
>> Best regards,
>> Raphael Mesaric
>> R-sig-Geo mailing list
>> [email protected]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: [email protected]
> https://orcid.org/-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0J=en
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
I am working on a spatial regression model which is built in such a way
that the dependent variable consists of n blocks of multiple entries,
and each block is referring to the same spatial grid. In order to run
the lagsarlm() function or the errsarlm() function, I now need to extend
my initial nb object from one grid to n grids. In other words, I need
the nb object to be repeated n times. However, so far I did not find a
way to do this properly.
I tried to just replicate the nb object by repeating the entries of the
neighbours list and the weights list, respectively.
dnn800_3w$neighbours
> To: "[email protected]"
> Subject: [R-sig-Geo] Question niche based model whith R
> Message-ID: <[email protected]>
> Content-Type: text/plain; charset="UTF-8"
> ?Hello dears,I work on the modeling of the ecological niche of striga, I
> have presence data (data) and bioclimatic variables and other variables
> around.Indeed, my question how should I proceed to model the ecological
> niche with software R??Pending receipt of my request, please accept my best
> regards.SADDA Abou-Soufianou --
> DoctorantUniversit? Dan Dicko Dankoulodo de Maradi,120 avenue Maman
> Koraou,ADS, MaradiBP 465 Maradi/Niger E-mail:[email protected]?:
> (+227)96269987
> [[alternative HTML version deleted]]
Hi Abou,
it is a quite broad question...a simple and easy answer would be "you
can just model yous data.frame with ?lm function or ?glm but you would
lose many of the iteresting stuff!! Then, in my opinion, first of all
you must decide between "presence-only" and "presence-absence" models
(my suggestion is for the second and in your case you must learn how
to generate pseudo-absences or background points. So, in my
experience, there are two main packages which are the most interesting
and suited for this purpose and which may help you: dismo and biomod2
(https://cran.r-project.org/web/packages/dismo/dismo.pdf -
https://cran.r-project.org/web/packages/biomod2/biomod2.pdf). There
are also many tutorials which can help you, for instance
https://cran.r-project.org/web/packages/dismo/vignettes/sdm.pdf
Then once you will handle the subject you will also be able to create
your own SDM using the algorithm you prefer (e.g. RandomForest,
MaxEnt, GLM, MARS, GAM and so on...) and without any additional
package. Anyway, my advice is to use one of the two above cited
packages because they are object-oriented and optimized for Ecological
Modelling.
Just to cover possible next question (which model is the best
solution) you will see that it is impossible to define the BEST method
which can perform correctly in all the circumstances. Many Researchers
are used to work with a single-algorithm (for example Generalized
Linear Models, Maximum Entropy, Generalized Additive Models, Neural
network, Random Forest etc.) while others are more willing to work
with ensemble projections and consensus models (i.e. mean or weighted
mean of more than one algorithm). Some examples of single-model and
ensemble model are here:
http://dx.doi.org/10.5424/fs/2016253-09476
http://dx.doi.org/10./gcb.12604
http://dx.doi.org/10.1007/s13595-014-0439-4
http://dx.doi.org/10./gcb.12476
Cheers,
Maurizio Marchi
Skype ID: maurizioxyz
linux user 552742
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> One more question, I assume that if I create independent lag variables
>> with: data$w_xxx > can get the impacts of these lag variables; however, if I add them into
>> model, do I need to change the type="lag" instead of "mixed" when running
>> lagsarlm command?
> Don't even think about this, it is senseless. The impact of x_r is as
> given in the formula, the only contrast is between the impacts in the
> model="lag" case and the model="mixed" case. Because of the interaction
> between \rho and the \betas (and \gammas) in the y = \rho W y + ... models,
> you only get one set of impacts per x_r per model. You can interprete the
> differences between the three impact components between models, but avoid
> playing with stuff without doing all of the maths first.
> Roger
>> Many thanks for your attention.
>> Jorge
>> *Ing. Jorge Alfredo Cárcamo, M. Sc., Ph. D. (c)*
>> Agriculture economics
>> Georg-August-Universität Göttingen
>> On Sun, Jun 5, 2016 at 2:17 PM, Roger Bivand
wrote:
>> On Sat, 4 Jun 2016, Jorge Cárcamo wrote:
>>> Dear all,
I am working with a lagsarlm "mixed" model I executed:
library(spdep)
library(coda)
dsts15 |z|)
(Intercept) 1.4255019 0.4824828 2.9545 0.003132
SD46 0.0149057 0.0116146 1.2834 0.199365
Totaland -0.0217556 0.0072475 -3.0018 0.002684
yearrain 0.0557373 0.0485612 1.1478 0.251062
yearfdi -0.0109979 0.0069536 -1.5816 0.113741
lag.SD46 0.0262273 0.0278045 0.9433 0.345539
lag.Totaland 0.0060992 0.0141515 0.4310 0.666474
lag.yearrain -0.2083118 0.0772751 -2.6957 0.007024
lag.yearfdi 0.0111460 0.0081562 1.3666 0.171761
Rho: -0.54702, LR test value: 10.052, p-value: 0.0015215
Immediatly, I exectued impacts(mod.sdm.15, R=1000), to get the impacts:
Direct Indirect Total
SD46 0.013136761 0.013451739 0.0265884995
Totaland -0.023517332 0.013397014 -0.0101203177
yearrain 0.079110119 -0.177734631 -0.0986245121
yearfdi -0.012677636 0.012773412 0.957761
Through simulation I manage to get credible intervals for direct,
indirect
and total impacts (HPDinterval(impacts, choice="XXX")).
>>> As you must know from the references on the help page for impacts
>>> methods,
>>> these are the combined impacts of the variables:
>>> S_r(W) = (I - \rho W)^{-1} (\beta_r I - \gamma_r W)
>>> where the direct impacts are sum(S_r(W))/n, etc. The \gamma_r are the
>>> coefficients on W x_r.
>>> However, how can I get the impacts of the lagged variables? I have seen
some publications such as: Läpple, D., & Kelley, H. (2014). Spatial
dependence in the adoption of organic drystock farming in Ireland. That
report a posterior mean and credible intervals for the lagged variables.
>>> Given the above, either you are misreading Läpple & Kelley (I do not have
>>> access), or both you and they are wrong. There are by definition on
>>> separable impacts for the lagged X variables.
>>> Hope this clarifies,
>>> Roger
>>> All suggestions are welcome.
Best regards,
*Ing. Jorge Alfredo Cárcamo, M. Sc., Ph. D. (c)*
Agriculture economics
Georg-August-Universität Göttingen
[[alternative HTML version deleted]]
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>> Roger Bivand
>>> Department of Economics, Norwegian School of Economics,
>>> Helleveien 30, N-5045 Bergen, Norway.
>>> voice: +47 55 95 93 55; fax +47 55 95 91 00
>>> e-mail: [email protected]
>>> http://orcid.org/-0003-2392-6140
>>> https://scholar.google.no/citations?user=AWeghB0J=en
>>> http://depsy.org/person/434412
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; fax +47 55 95 91 00
Re: [R-sig-Geo] Question on hierarchical data format 5 .h5 file format
2021-11-14 Thread Edzer Pebesma On 14/11/2021 13:59, Edzer Pebesma wrote: Gabriel shared the data with me (off-list, 1Gb) and it contains two hdf5 datasets, one containing the raster of the longitude and latitude, the other rasters with the attributres (LST) but no coordinates. I could get a (rough: factor 50 downsampling) plot with stars using: library(stars) lat = read_stars("HDF5:ECOSTRESS_L1B_GEO_11712_012_20200730T061029_0601_01.h5://Geolocation/latitude") lon = read_stars("HDF5:ECOSTRESS_L1B_GEO_11712_012_20200730T061029_0601_01.h5://Geolocation/longitude") var = "HDF5:ECOSTRESS_L2_LSTE_11712_012_20200730T061029_0601_01.h5://SDS/LST" lst = read_stars(var, curvilinear = list(x = lon, y = lat)) plot(lst, downsample = 50, axes = TRUE, reset = FALSE) library(rnaturalearth) ne = ne_countries(returnclass = "sf") plot(st_geometry(ne), add = TRUE, border = 'yellow') # california border seems to make sense but I'm not sure how to make a regular raster (in some CRS) out of this for exporting to GeoTIFF. This might be useful: https://github.com/r-spatial/stars/issues/386 but converting to polygons will be incredibly slow and memory hungry. On 14/11/2021 11:00, Edzer Pebesma wrote: This might be because the grid is not regular but possibly curvilinear, having two raster layers with the lon & lat values for each pixel. Can you share a sample data set? On 14/11/2021 10:47, Gabriel Cotlier wrote: Dear Michael, Following your advice I have tired : library(terra) library(raster) library(dplyr) r = terra::rast(file.choose()) plot(r$LST) r_lst = r$LST r_lst %>% raster() plot(r_lst) Then I effectively get a raster layer format from the raster package wanted, but looking at the plot and the raster layer information--see bellow-- the scene appears to need both latitude and longitude for correct geolocation, how can I do this procedure? r_lst %>% raster() class : RasterLayer dimensions : 5632, 5400, 30412800 (nrow, ncol, ncell) resolution : 1, 1 (x, y) extent : 0, 5400, 0, 5632 (xmin, xmax, ymin, ymax) crs : NA source : LST names : LST values : 0, 65535 (min, max) The extent and the resolution are not the correct ones. Maybe in the same .h5 is all the data for georeferencing it although I could not see how to get to it and apply it to the layer. There is also another file with complementary data related to geolocation but how can I query it, and get the data for geolocation and apply it to the raster layer?. The objective is to save it as a GeoTIFF file to be opened in QGIS with correct geolocation. What do you suggest for converting the .h5 file to a georeferenced Raster Layer in GeoTIFF format ? Thanks a lot. Kind regards, Gabriel On Sat, Nov 13, 2021 at 12:22 PM Gabriel Cotlier wrote: Dear Michael Summer, Thanks a lot for your swift reply. I will give a try to the procedure you mentioned. Kind regards, Gabriel On Sat, Nov 13, 2021 at 12:14 PM Michael Sumner wrote: try terra::rast() or stars::read_stars() on the file, both use gdal to interrogate and read. What you get depends on the structure of the files, variables interpreted as subdatasets. stars is more general but sees variables as 2D arrays with bands, terra calls these bands layers. Alternatively try raster::raster() or RNetCDF or tidync which might be more straightforward than rhdf5 but all are quite different to each other and to the gdal view. On Sat, 13 Nov 2021, 18:56 Gabriel Cotlier, wrote: Hello, I would like to be able to visualize / plot a specific data layer within the structure a hierarchical data file format .h5 in R. I could observe the structure as follows : library(rhdf5) r= h5ls(file.choose()) and also in this way identify the Group and Name of the data to be plotted, however I've tried I could not access the data itself, nor plot it raster layer with its corresponding geographic coordinates. I would like to be able to select such data, plot it as a raster layer if possible save it as GeoTIFF file format to be further processed raster package as a common raster layer. If any small guidance example, reprex, of this procedure is possible it would be appreciated. Thanks a lot in advance Kind regards, Gabriel Cotlier Gabriel Cotlier, PhD Haifa Research Center for Theoretical Physics and Astrophysics (HCTPA) University of Haifa Israel [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 R-sig-Geo mailing list [email protected]Re: [R-sig-Geo] Question on hierarchical data format 5 .h5 file format
2021-11-14 Thread Gabriel Cotlier Hello Edzer, Thanks a lot. The scaling factor to convert to LST is 0.02 Regards, Gabriel On Sun, Nov 14, 2021 at 2:59 PM Edzer Pebesma wrote: > Gabriel shared the data with me (off-list, 1Gb) and it contains two hdf5 > datasets, one containing the raster of the longitude and latitude, the > other rasters with the attributres (LST) but no coordinates. I could get > a (rough: factor 50 downsampling) plot with stars using: > library(stars) > lat = > read_stars("HDF5:ECOSTRESS_L1B_GEO_11712_012_20200730T061029_0601_01.h5://Geolocation/latitude") > lon = > read_stars("HDF5:ECOSTRESS_L1B_GEO_11712_012_20200730T061029_0601_01.h5://Geolocation/longitude") > var = > "HDF5:ECOSTRESS_L2_LSTE_11712_012_20200730T061029_0601_01.h5://SDS/LST" > lst = read_stars(var, curvilinear = list(x = lon, y = lat)) > plot(lst, downsample = 50, axes = TRUE, reset = FALSE) > library(rnaturalearth) > ne = ne_countries(returnclass = "sf") > plot(st_geometry(ne), add = TRUE, border = 'yellow') > # california border seems to make sense > but I'm not sure how to make a regular raster (in some CRS) out of this > for exporting to GeoTIFF. > On 14/11/2021 11:00, Edzer Pebesma wrote: > > This might be because the grid is not regular but possibly curvilinear, > > having two raster layers with the lon & lat values for each pixel. Can > > you share a sample data set? > > On 14/11/2021 10:47, Gabriel Cotlier wrote: > >> Dear Michael, > >> Following your advice I have tired : > >> library(terra) > >> library(raster) > >> library(dplyr) > >> r = terra::rast(file.choose()) > >> plot(r$LST) > >> r_lst = r$LST > >> r_lst %>% raster() > >> plot(r_lst) > >> Then I effectively get a raster layer format from the raster package as > >> wanted, but looking at the plot and the raster layer information--see > >> bellow-- the scene appears to need both latitude and longitude for > >> correct > >> geolocation, how can I do this procedure? > >>> r_lst %>% raster() > >> class : RasterLayer > >> dimensions : 5632, 5400, 30412800 (nrow, ncol, ncell) > >> resolution : 1, 1 (x, y) > >> extent : 0, 5400, 0, 5632 (xmin, xmax, ymin, ymax) > >> crs: NA > >> source : LST > >> names : LST > >> values : 0, 65535 (min, max) > >> The extent and the resolution are not the correct ones. > >> Maybe in the same .h5 is all the data for georeferencing it although I > >> could not see how to get to it and apply it to the layer. > >> There is also another file with complementary data related to > geolocation > >> but how can I query it, and get the data for geolocation and apply it to > >> the raster layer?. > >> The objective is to save it as a GeoTIFF file to be opened in QGIS with > >> correct geolocation. > >> What do you suggest for converting the .h5 file to a georeferenced > Raster > >> Layer in GeoTIFF format ? > >> Thanks a lot. > >> Kind regards, > >> Gabriel > >> On Sat, Nov 13, 2021 at 12:22 PM Gabriel Cotlier > >> wrote: > >>> Dear Michael Summer, > >>> Thanks a lot for your swift reply. > >>> I will give a try to the procedure you mentioned. > >>> Kind regards, > >>> Gabriel > >>> On Sat, Nov 13, 2021 at 12:14 PM Michael Sumner > >>> wrote: > try terra::rast() or stars::read_stars() on the file, both use gdal to > interrogate and read. What you get depends on the structure of the > files, > variables interpreted as subdatasets. stars is more general but sees > variables as 2D arrays with bands, terra calls these bands layers. > Alternatively try raster::raster() or RNetCDF or tidync which might be > more straightforward than rhdf5 but all are quite different to each > other > and to the gdal view. > On Sat, 13 Nov 2021, 18:56 Gabriel Cotlier, > wrote: > > Hello, > > I would like to be able to visualize / plot a specific data layer > > within > > the structure a hierarchical data file format .h5 in R. I could > > observe > > the structure as follows : > > library(rhdf5) > > r= h5ls(file.choose()) > > and also in this way identify the Group and Name of the data to be > > plotted, > > however I've tried I could not access the data itself, nor plot it > > raster layer with its corresponding geographic coordinates. > > I would like to be able to select such data, plot it as a raster > layer > > if possible save it as GeoTIFF file format to be further processed > > raster package as a common raster layer. > > If any small guidance example, reprex, of this procedure is possible > > it would be appreciated. > > Thanks a lot in advance > > Kind regards, > > Gabriel Cotlier > > Gabriel Cotlier, PhD > > Haifa Research Center for Theoretical Physics andRe: [R-sig-Geo] Question on hierarchical data format 5 .h5 file format
2021-11-14 Thread Edzer Pebesma Gabriel shared the data with me (off-list, 1Gb) and it contains two hdf5 datasets, one containing the raster of the longitude and latitude, the other rasters with the attributres (LST) but no coordinates. I could get a (rough: factor 50 downsampling) plot with stars using: library(stars) lat = read_stars("HDF5:ECOSTRESS_L1B_GEO_11712_012_20200730T061029_0601_01.h5://Geolocation/latitude") lon = read_stars("HDF5:ECOSTRESS_L1B_GEO_11712_012_20200730T061029_0601_01.h5://Geolocation/longitude") var = "HDF5:ECOSTRESS_L2_LSTE_11712_012_20200730T061029_0601_01.h5://SDS/LST" lst = read_stars(var, curvilinear = list(x = lon, y = lat)) plot(lst, downsample = 50, axes = TRUE, reset = FALSE) library(rnaturalearth) ne = ne_countries(returnclass = "sf") plot(st_geometry(ne), add = TRUE, border = 'yellow') # california border seems to make sense but I'm not sure how to make a regular raster (in some CRS) out of this for exporting to GeoTIFF. On 14/11/2021 11:00, Edzer Pebesma wrote: This might be because the grid is not regular but possibly curvilinear, having two raster layers with the lon & lat values for each pixel. Can you share a sample data set? On 14/11/2021 10:47, Gabriel Cotlier wrote: Dear Michael, Following your advice I have tired : library(terra) library(raster) library(dplyr) r = terra::rast(file.choose()) plot(r$LST) r_lst = r$LST r_lst %>% raster() plot(r_lst) Then I effectively get a raster layer format from the raster package as I wanted, but looking at the plot and the raster layer information--see bellow-- the scene appears to need both latitude and longitude for correct geolocation, how can I do this procedure? r_lst %>% raster() class : RasterLayer dimensions : 5632, 5400, 30412800 (nrow, ncol, ncell) resolution : 1, 1 (x, y) extent : 0, 5400, 0, 5632 (xmin, xmax, ymin, ymax) crs : NA source : LST names : LST values : 0, 65535 (min, max) The extent and the resolution are not the correct ones. Maybe in the same .h5 is all the data for georeferencing it although I could not see how to get to it and apply it to the layer. There is also another file with complementary data related to geolocation but how can I query it, and get the data for geolocation and apply it to the raster layer?. The objective is to save it as a GeoTIFF file to be opened in QGIS with correct geolocation. What do you suggest for converting the .h5 file to a georeferenced Raster Layer in GeoTIFF format ? Thanks a lot. Kind regards, Gabriel On Sat, Nov 13, 2021 at 12:22 PM Gabriel Cotlier wrote: Dear Michael Summer, Thanks a lot for your swift reply. I will give a try to the procedure you mentioned. Kind regards, Gabriel On Sat, Nov 13, 2021 at 12:14 PM Michael Sumner wrote: try terra::rast() or stars::read_stars() on the file, both use gdal to interrogate and read. What you get depends on the structure of the files, variables interpreted as subdatasets. stars is more general but sees variables as 2D arrays with bands, terra calls these bands layers. Alternatively try raster::raster() or RNetCDF or tidync which might be more straightforward than rhdf5 but all are quite different to each other and to the gdal view. On Sat, 13 Nov 2021, 18:56 Gabriel Cotlier, wrote: Hello, I would like to be able to visualize / plot a specific data layer within the structure a hierarchical data file format .h5 in R. I could observe the structure as follows : library(rhdf5) r= h5ls(file.choose()) and also in this way identify the Group and Name of the data to be plotted, however I've tried I could not access the data itself, nor plot it raster layer with its corresponding geographic coordinates. I would like to be able to select such data, plot it as a raster layer if possible save it as GeoTIFF file format to be further processed raster package as a common raster layer. If any small guidance example, reprex, of this procedure is possible it would be appreciated. Thanks a lot in advance Kind regards, Gabriel Cotlier Gabriel Cotlier, PhD Haifa Research Center for Theoretical Physics and Astrophysics (HCTPA) University of Haifa Israel [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo This might be because the grid is not regular but possibly curvilinear, having two raster layers with the lon & lat values for each pixel. Can you share a sample data set? On 14/11/2021 10:47, Gabriel Cotlier wrote: Dear Michael, Following your advice I have tired : library(terra) library(raster) library(dplyr) r = terra::rast(file.choose()) plot(r$LST) r_lst = r$LST r_lst %>% raster() plot(r_lst) Then I effectively get a raster layer format from the raster package as I wanted, but looking at the plot and the raster layer information--see bellow-- the scene appears to need both latitude and longitude for correct geolocation, how can I do this procedure? r_lst %>% raster() class : RasterLayer dimensions : 5632, 5400, 30412800 (nrow, ncol, ncell) resolution : 1, 1 (x, y) extent : 0, 5400, 0, 5632 (xmin, xmax, ymin, ymax) crs: NA source : LST names : LST values : 0, 65535 (min, max) The extent and the resolution are not the correct ones. Maybe in the same .h5 is all the data for georeferencing it although I could not see how to get to it and apply it to the layer. There is also another file with complementary data related to geolocation but how can I query it, and get the data for geolocation and apply it to the raster layer?. The objective is to save it as a GeoTIFF file to be opened in QGIS with correct geolocation. What do you suggest for converting the .h5 file to a georeferenced Raster Layer in GeoTIFF format ? Thanks a lot. Kind regards, Gabriel On Sat, Nov 13, 2021 at 12:22 PM Gabriel Cotlier wrote: Dear Michael Summer, Thanks a lot for your swift reply. I will give a try to the procedure you mentioned. Kind regards, Gabriel On Sat, Nov 13, 2021 at 12:14 PM Michael Sumner wrote: try terra::rast() or stars::read_stars() on the file, both use gdal to interrogate and read. What you get depends on the structure of the files, variables interpreted as subdatasets. stars is more general but sees variables as 2D arrays with bands, terra calls these bands layers. Alternatively try raster::raster() or RNetCDF or tidync which might be more straightforward than rhdf5 but all are quite different to each other and to the gdal view. On Sat, 13 Nov 2021, 18:56 Gabriel Cotlier, wrote: Hello, I would like to be able to visualize / plot a specific data layer within the structure a hierarchical data file format .h5 in R. I could observe the structure as follows : library(rhdf5) r= h5ls(file.choose()) and also in this way identify the Group and Name of the data to be plotted, however I've tried I could not access the data itself, nor plot it as a raster layer with its corresponding geographic coordinates. I would like to be able to select such data, plot it as a raster layer if possible save it as GeoTIFF file format to be further processed with raster package as a common raster layer. If any small guidance example, reprex, of this procedure is possible it would be appreciated. Thanks a lot in advance Kind regards, Gabriel Cotlier Gabriel Cotlier, PhD Haifa Research Center for Theoretical Physics and Astrophysics (HCTPA) University of Haifa Israel [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Then I effectively get a raster layer format from the raster package as I wanted, but looking at the plot and the raster layer information--see bellow-- the scene appears to need both latitude and longitude for correct geolocation, how can I do this procedure? > r_lst %>% raster() class : RasterLayer dimensions : 5632, 5400, 30412800 (nrow, ncol, ncell) resolution : 1, 1 (x, y) extent : 0, 5400, 0, 5632 (xmin, xmax, ymin, ymax) crs: NA source : LST names : LST values : 0, 65535 (min, max) The extent and the resolution are not the correct ones. Maybe in the same .h5 is all the data for georeferencing it although I could not see how to get to it and apply it to the layer. There is also another file with complementary data related to geolocation but how can I query it, and get the data for geolocation and apply it to the raster layer?. The objective is to save it as a GeoTIFF file to be opened in QGIS with correct geolocation. What do you suggest for converting the .h5 file to a georeferenced Raster Layer in GeoTIFF format ? Thanks a lot. Kind regards, Gabriel On Sat, Nov 13, 2021 at 12:22 PM Gabriel Cotlier wrote: > Dear Michael Summer, > Thanks a lot for your swift reply. > I will give a try to the procedure you mentioned. > Kind regards, > Gabriel > On Sat, Nov 13, 2021 at 12:14 PM Michael Sumner > wrote: >> try terra::rast() or stars::read_stars() on the file, both use gdal to >> interrogate and read. What you get depends on the structure of the files, >> variables interpreted as subdatasets. stars is more general but sees >> variables as 2D arrays with bands, terra calls these bands layers. >> Alternatively try raster::raster() or RNetCDF or tidync which might be >> more straightforward than rhdf5 but all are quite different to each other >> and to the gdal view. >> On Sat, 13 Nov 2021, 18:56 Gabriel Cotlier, wrote: >>> Hello, >>> I would like to be able to visualize / plot a specific data layer within >>> the structure a hierarchical data file format .h5 in R. I could observe >>> the structure as follows : >>> library(rhdf5) >>> r= h5ls(file.choose()) >>> and also in this way identify the Group and Name of the data to be >>> plotted, >>> however I've tried I could not access the data itself, nor plot it as a >>> raster layer with its corresponding geographic coordinates. >>> I would like to be able to select such data, plot it as a raster layer >>> if possible save it as GeoTIFF file format to be further processed with >>> raster package as a common raster layer. >>> If any small guidance example, reprex, of this procedure is possible >>> it would be appreciated. >>> Thanks a lot in advance >>> Kind regards, >>> Gabriel Cotlier >>> Gabriel Cotlier, PhD >>> Haifa Research Center for Theoretical Physics and Astrophysics (HCTPA) >>> University of Haifa >>> Israel >>> [[alternative HTML version deleted]] >>> R-sig-Geo mailing list >>> [email protected] >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Dear Michael Summer, Thanks a lot for your swift reply. I will give a try to the procedure you mentioned. Kind regards, Gabriel On Sat, Nov 13, 2021 at 12:14 PM Michael Sumner wrote: > try terra::rast() or stars::read_stars() on the file, both use gdal to > interrogate and read. What you get depends on the structure of the files, > variables interpreted as subdatasets. stars is more general but sees > variables as 2D arrays with bands, terra calls these bands layers. > Alternatively try raster::raster() or RNetCDF or tidync which might be > more straightforward than rhdf5 but all are quite different to each other > and to the gdal view. > On Sat, 13 Nov 2021, 18:56 Gabriel Cotlier, wrote: >> Hello, >> I would like to be able to visualize / plot a specific data layer within >> the structure a hierarchical data file format .h5 in R. I could observe >> the structure as follows : >> library(rhdf5) >> r= h5ls(file.choose()) >> and also in this way identify the Group and Name of the data to be >> plotted, >> however I've tried I could not access the data itself, nor plot it as a >> raster layer with its corresponding geographic coordinates. >> I would like to be able to select such data, plot it as a raster layer and >> if possible save it as GeoTIFF file format to be further processed with >> raster package as a common raster layer. >> If any small guidance example, reprex, of this procedure is possible >> it would be appreciated. >> Thanks a lot in advance >> Kind regards, >> Gabriel Cotlier >> Gabriel Cotlier, PhD >> Haifa Research Center for Theoretical Physics and Astrophysics (HCTPA) >> University of Haifa >> Israel >> [[alternative HTML version deleted]] >> R-sig-Geo mailing list >> [email protected] >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo try terra::rast() or stars::read_stars() on the file, both use gdal to interrogate and read. What you get depends on the structure of the files, variables interpreted as subdatasets. stars is more general but sees variables as 2D arrays with bands, terra calls these bands layers. Alternatively try raster::raster() or RNetCDF or tidync which might be more straightforward than rhdf5 but all are quite different to each other and to the gdal view. On Sat, 13 Nov 2021, 18:56 Gabriel Cotlier, wrote: > Hello, > I would like to be able to visualize / plot a specific data layer within > the structure a hierarchical data file format .h5 in R. I could observe > the structure as follows : > library(rhdf5) > r= h5ls(file.choose()) > and also in this way identify the Group and Name of the data to be plotted, > however I've tried I could not access the data itself, nor plot it as a > raster layer with its corresponding geographic coordinates. > I would like to be able to select such data, plot it as a raster layer and > if possible save it as GeoTIFF file format to be further processed with the > raster package as a common raster layer. > If any small guidance example, reprex, of this procedure is possible > it would be appreciated. > Thanks a lot in advance > Kind regards, > Gabriel Cotlier > Gabriel Cotlier, PhD > Haifa Research Center for Theoretical Physics and Astrophysics (HCTPA) > University of Haifa > Israel > [[alternative HTML version deleted]] > R-sig-Geo mailing list > [email protected] > https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Thanks a lot. I will try spdep::n.comp.nb() and read the link you suggested to me on "Proximity and aerial data" which seems very interesting. Maybe could it be due to the fact that the data layer of the city is traversed by a river, which I set to na since it was not relevant for the analysis, thus maybe from there the message "reginos with no links"? However when I do : Neighbour list object: Number of regions: 636410 Number of nonzero links: 5060866 Percentage nonzero weights: 0.001249542 Average number of links: 7.95221 8 regions with no links: 103556 104030 104502 104970 170752 290616 291628 520687 Then, even with regions with no links, I have no error when setting zero.policy = TRUE in nb2listw() > lw.l moran.test(l$dta, lw.l,zero.policy = TRUE) # default randomisation=TRUE Moran I test under randomisation data: l$dta weights: lw.l n reduced by no-neighbour observations Moran I statistic standard deviate = 1559.1, p-value < 2.2e-16 alternative hypothesis: greater sample estimates: Moran I statistic Expectation Variance 9.815120e-01 -1.571336e-06 3.963186e-07 Why could it be running apparently ok in this way ? Thanks a lot again for your help and guidance. Kind regards, Gabriel Cotlier On Mon, Oct 25, 2021 at 2:12 PM Roger Bivand wrote: > On Mon, 25 Oct 2021, Gabriel Cotlier wrote: > > Hello, > > I have tried the code for most of the dataset but only few make the > > same error: > > Neighbour list object: > > Number of regions: 975434 > > Number of nonzero links: 3869290 > > Percentage nonzero weights: 0.0004066638 > > Average number of links: 3.966737 > > 79 regions with no links: > > 3206 34632 47267 65155 81013 90934 91361 91427 98318 > > 98765 100054 100472 100884 100896 102508 102514 105207 > > 105937 106630 106965 114449 115479 122329 123614 126880 > > 128631 314459 350405 381283 397379 401275 412489 420698 > > 426384 428632 429752 432058 433230 440631 451195 451196 > > 452270 464327 467009 503242 504459 507625 510412 517865 > > 517916 524413 526766 535269 537320 537712 544576 549519 > > 551629 553413 565169 566421 567703 567791 568587 579209 > > 585464 625809 642179 648122 651096 652581 654065 708908 > > 777500 845322 846238 847153 848065 866391 > >> lw.l > Error in nb2listw(nb.l, style = "W") : Empty neighbour sets found > > Could it be something related to the input data layer? > Yes. Please also look at the output of spdep::n.comp.nb() if it works for > a dataset this large, see ?n.comp.nb and examine the nc component (number > of subgraphs). The singleton graphs have no neighbours, but there may be > other subgraphs, you find their sizes by table(table( component>)). This may well have come from the na.rm=TRUE in > rasterToPolygons(), but is unavoidable. Maybe see also: > https://keen-swartz-3146c4.netlify.app/area.html# > and expand all code (button at top right) and search for n.comp.nc. > Since you have singleton observations, zero.policy=TRUE is your only > option unless you further subset by omitting all the observations for > which card(nb.l) == 0L - since they have no neighbours, they do not inform > the spatial process. > Roger > > How could it be possible to avoid this error for running the code? > > Thanks a lot for your guidance and help. > > Kind regards, > > Gabriel Cotlier > > On Mon, Oct 25, 2021 at 12:21 PM Gabriel Cotlier > > wrote: > >> Thank you very much for your response. > >> Yes indeed it clarifies and helps a lot. > >> Kind regards, > >> Gabriel Cotlier > >> On Mon, Oct 25, 2021 at 12:00 PM Roger Bivand > wrote: > >>> On Sat, 23 Oct 2021, Gabriel Cotlier wrote: > Hello, > I would like to estimate the pseudo p-value of Moran's I from raster > layer, by converting it to polygon data in order to be able to > directly > the function spdep::mc.moran() and getting in one step not only value > Moran'sI but its statistical significance without using other > >>> simulations > but only that in the funcion spdep::mc.moran() as an alternative I > used the code below. However, since with the initial version of the > >>> code I > got several errors such as : > Error in moran. mc(My_raster, lw.l, 999) : r is not a numeric vector > So I converted "r" to a vector using as.vector() in : > >>> First, if you have the spatial weights object lw.l, > spdep::moran.test(), > >>> using analytical rather than bootstrap inference, is almost certainly > >>> efficient, and yields almost always the same outcome with > >>> randomisation=TRUE. > However then appeared another the error: > Error in na.fail.default(x) : missing values inRe: [R-sig-Geo] question on raster Moran's I statistical significance
2021-10-25 Thread Roger Bivand On Mon, 25 Oct 2021, Roger Bivand wrote: On Mon, 25 Oct 2021, Gabriel Cotlier wrote: Hello, I have tried the code for most of the dataset but only few make the same error: Neighbour list object: Number of regions: 975434 Number of nonzero links: 3869290 Percentage nonzero weights: 0.0004066638 Average number of links: 3.966737 79 regions with no links: 3206 34632 47267 65155 81013 90934 91361 91427 98318 98765 100054 100472 100884 100896 102508 102514 105207 105937 106630 106965 114449 115479 122329 123614 126880 128631 314459 350405 381283 397379 401275 412489 420698 426384 428632 429752 432058 433230 440631 451195 451196 452270 464327 467009 503242 504459 507625 510412 517865 517916 524413 526766 535269 537320 537712 544576 549519 551629 553413 565169 566421 567703 567791 568587 579209 585464 625809 642179 648122 651096 652581 654065 708908 777500 845322 846238 847153 848065 866391 lw.l )) . This may well have come from the na.rm=TRUE in rasterToPolygons(), but is unavoidable. Maybe see also: https://keen-swartz-3146c4.netlify.app/area.html# and expand all code (button at top right) and search for n.comp.nc. Since you have singleton observations, zero.policy=TRUE is your only option unless you further subset by omitting all the observations for which card(nb.l) == 0L - since they have no neighbours, they do not inform the spatial process. Roger How could it be possible to avoid this error for running the code? Thanks a lot for your guidance and help. Kind regards, Gabriel Cotlier On Mon, Oct 25, 2021 at 12:21 PM Gabriel Cotlier wrote: Thank you very much for your response. Yes indeed it clarifies and helps a lot. Kind regards, Gabriel Cotlier On Mon, Oct 25, 2021 at 12:00 PM Roger Bivand wrote: On Sat, 23 Oct 2021, Gabriel Cotlier wrote: Hello, I would like to estimate the pseudo p-value of Moran's I from raster layer, by converting it to polygon data in order to be able to directly the function spdep::mc.moran() and getting in one step not only value Moran'sI but its statistical significance without using other simulations but only that in the funcion spdep::mc.moran() as an alternative I have used the code below. However, since with the initial version of the code I got several errors such as : Error in moran. mc(My_raster, lw.l, 999) : r is not a numeric vector So I converted "r" to a vector using as.vector() in :Re: [R-sig-Geo] question on raster Moran's I statistical significance
2021-10-25 Thread Roger Bivand On Mon, 25 Oct 2021, Gabriel Cotlier wrote: Hello, I have tried the code for most of the dataset but only few make the same error: Neighbour list object: Number of regions: 975434 Number of nonzero links: 3869290 Percentage nonzero weights: 0.0004066638 Average number of links: 3.966737 79 regions with no links: 3206 34632 47267 65155 81013 90934 91361 91427 98318 98765 100054 100472 100884 100896 102508 102514 105207 105937 106630 106965 114449 115479 122329 123614 126880 128631 314459 350405 381283 397379 401275 412489 420698 426384 428632 429752 432058 433230 440631 451195 451196 452270 464327 467009 503242 504459 507625 510412 517865 517916 524413 526766 535269 537320 537712 544576 549519 551629 553413 565169 566421 567703 567791 568587 579209 585464 625809 642179 648122 651096 652581 654065 708908 777500 845322 846238 847153 848065 866391 lw.l )). This may well have come from the na.rm=TRUE in rasterToPolygons(), but is unavoidable. Maybe see also: https://keen-swartz-3146c4.netlify.app/area.html# and expand all code (button at top right) and search for n.comp.nc. Since you have singleton observations, zero.policy=TRUE is your only option unless you further subset by omitting all the observations for which card(nb.l) == 0L - since they have no neighbours, they do not inform the spatial process. Roger How could it be possible to avoid this error for running the code? Thanks a lot for your guidance and help. Kind regards, Gabriel Cotlier On Mon, Oct 25, 2021 at 12:21 PM Gabriel Cotlier wrote: Thank you very much for your response. Yes indeed it clarifies and helps a lot. Kind regards, Gabriel Cotlier On Mon, Oct 25, 2021 at 12:00 PM Roger Bivand wrote: On Sat, 23 Oct 2021, Gabriel Cotlier wrote: Hello, I would like to estimate the pseudo p-value of Moran's I from raster layer, by converting it to polygon data in order to be able to directly the function spdep::mc.moran() and getting in one step not only value Moran'sI but its statistical significance without using other simulations but only that in the funcion spdep::mc.moran() as an alternative I have used the code below. However, since with the initial version of the code I got several errors such as : Error in moran. mc(My_raster, lw.l, 999) : r is not a numeric vector So I converted "r" to a vector using as.vector() in :Re: [R-sig-Geo] question on raster Moran's I statistical significance
2021-10-25 Thread Gabriel Cotlier Hello, I have tried the code for most of the dataset but only few make the same error: Neighbour list object: Number of regions: 975434 Number of nonzero links: 3869290 Percentage nonzero weights: 0.0004066638 Average number of links: 3.966737 79 regions with no links: 3206 34632 47267 65155 81013 90934 91361 91427 98318 98765 100054 100472 100884 100896 102508 102514 105207 105937 106630 106965 114449 115479 122329 123614 126880 128631 314459 350405 381283 397379 401275 412489 420698 426384 428632 429752 432058 433230 440631 451195 451196 452270 464327 467009 503242 504459 507625 510412 517865 517916 524413 526766 535269 537320 537712 544576 549519 551629 553413 565169 566421 567703 567791 568587 579209 585464 625809 642179 648122 651096 652581 654065 708908 777500 845322 846238 847153 848065 866391 > lw.l Thank you very much for your response. > Yes indeed it clarifies and helps a lot. > Kind regards, > Gabriel Cotlier > On Mon, Oct 25, 2021 at 12:00 PM Roger Bivand wrote: >> On Sat, 23 Oct 2021, Gabriel Cotlier wrote: >> > Hello, >> > I would like to estimate the pseudo p-value of Moran's I from raster >> > layer, by converting it to polygon data in order to be able to directly >> > the function spdep::mc.moran() and getting in one step not only value >> > Moran'sI but its statistical significance without using other >> simulations >> > but only that in the funcion spdep::mc.moran() as an alternative I have >> > used the code below. However, since with the initial version of the >> code I >> > got several errors such as : >> > Error in moran. mc(My_raster, lw.l, 999) : r is not a numeric vector >> > So I converted "r" to a vector using as.vector() in : >> First, if you have the spatial weights object lw.l, spdep::moran.test(), >> using analytical rather than bootstrap inference, is almost certainly >> efficient, and yields almost always the same outcome with >> randomisation=TRUE. >> > However then appeared another the error: >> > Error in na.fail.default(x) : missing values in object >> Please see the na.action= argument to moran.mc() and moran.test(); the >> default is na.fail, but can be na.exclude, which then also subsets the >> listw spatial weights object. You can also subset l before creating the >> neighbour object: >> library(raster) >> library(spdep) >> r > values(r) > l > dim(l) # [1] 99 1 >> nb.l > nb.l # here no no-neighbour observations, no need for zero.policy= >> lw.l > moran.test(l$layer, lw.l) # default randomisation=TRUE >> > Finally these errors disappeared with the following version of the >> code-- >> > see below -- the version of the code below runs OK--unless for the input >> > raster I have used-- with no more errors or problems. >> Because rasterToPolygons() removes NA observations by default. >> Hope this clarifies, >> Roger >> > Therefore my question is : if instead of using arguments of the function >> > spdep::mc.moran() such as na.action = na.omit or other --which I have >> tried >> > and could not make work-- is it a valid way to solve the above errors >> > the way as presented in the code below? >> > I would appreciate it a lot if there is a possibility of knowing what >> could >> > be wrong in the code below? >> > I suspect it could be possible the problem is related to the way in >> which >> > the raster data layer is converted to polygon data im the code to be >> > input to the function spdep::mc.moran() , but will appreciate any >> guidance >> > and assistance in the issue to get working in the correct manner. >> > Thanks a lot. >> > Kind regards, >> > ### Estimate of Moran's I and its p-value from a raster >> > layer ## >> > library(raster) >> > library(spdep) >> > library(easycsv) >> > rm(list = ls()) >> > nb.l > > lw.l > > >> > v> > M > > M >> > plot(M, main=sprintf("Original, Moran's I (p = %0.3f)", M$p.value), >> > xlab="Monte-Carlo Simulation of Moran I", sub =NA, cex.lab=1.5, >> > cex.axis=1.5, cex.main=1.5, cex.sub=1.5) >> > abline(v=M$statistic, col="red", lw=2) >> > dev.off()Re: [R-sig-Geo] question on raster Moran's I statistical significance
2021-10-25 Thread Gabriel Cotlier Thank you very much for your response. Yes indeed it clarifies and helps a lot. Kind regards, Gabriel Cotlier On Mon, Oct 25, 2021 at 12:00 PM Roger Bivand wrote: > On Sat, 23 Oct 2021, Gabriel Cotlier wrote: > > Hello, > > I would like to estimate the pseudo p-value of Moran's I from raster data > > layer, by converting it to polygon data in order to be able to directly > > the function spdep::mc.moran() and getting in one step not only value of > > Moran'sI but its statistical significance without using other simulations > > but only that in the funcion spdep::mc.moran() as an alternative I have > > used the code below. However, since with the initial version of the code > > got several errors such as : > > Error in moran. mc(My_raster, lw.l, 999) : r is not a numeric vector > > So I converted "r" to a vector using as.vector() in : > First, if you have the spatial weights object lw.l, spdep::moran.test(), > using analytical rather than bootstrap inference, is almost certainly more > efficient, and yields almost always the same outcome with > randomisation=TRUE. > > However then appeared another the error: > > Error in na.fail.default(x) : missing values in object > Please see the na.action= argument to moran.mc() and moran.test(); the > default is na.fail, but can be na.exclude, which then also subsets the > listw spatial weights object. You can also subset l before creating the > neighbour object: > library(raster) > library(spdep) > r values(r) l dim(l) # [1] 99 1 > nb.l nb.l # here no no-neighbour observations, no need for zero.policy= > lw.l moran.test(l$layer, lw.l) # default randomisation=TRUE > > Finally these errors disappeared with the following version of the code-- > > see below -- the version of the code below runs OK--unless for the input > > raster I have used-- with no more errors or problems. > Because rasterToPolygons() removes NA observations by default. > Hope this clarifies, > Roger > > Therefore my question is : if instead of using arguments of the function > > spdep::mc.moran() such as na.action = na.omit or other --which I have > tried > > and could not make work-- is it a valid way to solve the above errors in > > the way as presented in the code below? > > I would appreciate it a lot if there is a possibility of knowing what > could > > be wrong in the code below? > > I suspect it could be possible the problem is related to the way in which > > the raster data layer is converted to polygon data im the code to be > > input to the function spdep::mc.moran() , but will appreciate any > guidance > > and assistance in the issue to get working in the correct manner. > > Thanks a lot. > > Kind regards, > > ### Estimate of Moran's I and its p-value from a raster > > layer ## > > library(raster) > > library(spdep) > > library(easycsv) > > rm(list = ls()) > > nb.l > lw.l > > > v > M > M > > plot(M, main=sprintf("Original, Moran's I (p = %0.3f)", M$p.value), > > xlab="Monte-Carlo Simulation of Moran I", sub =NA, cex.lab=1.5, > > cex.axis=1.5, cex.main=1.5, cex.sub=1.5) > > abline(v=M$statistic, col="red", lw=2) > > dev.off() > Emeritus Professor > Department of Economics, Norwegian School of Economics, > Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway. > e-mail: [email protected] > https://orcid.org/-0003-2392-6140 > https://scholar.google.no/citations?user=AWeghB0J=en Gabriel Cotlier, PhD Haifa Research Center for Theoretical Physics and Astrophysics (HCTPA) University of Haifa Israel [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo I would like to estimate the pseudo p-value of Moran's I from raster data layer, by converting it to polygon data in order to be able to directly use the function spdep::mc.moran() and getting in one step not only value of Moran'sI but its statistical significance without using other simulations but only that in the funcion spdep::mc.moran() as an alternative I have used the code below. However, since with the initial version of the code I got several errors such as : Error in moran. mc(My_raster, lw.l, 999) : r is not a numeric vector So I converted "r" to a vector using as.vector() in :Re: [R-sig-Geo] question on raster Moran's I statistical significance
2021-10-22 Thread Gabriel Cotlier Hello, With regards to my intention of solely estimating the pseudo p-value of Moran's I from raster data layer, by converting it to polygon data in order to be able to directly use the function spdep::mc.moran() and getting in one step not only value of Moran'sI but its statistical significance without using other simulations but only that in the funcion spdep::mc.moran() as an alternative I have used the code below. However, since with the initial version of the code I got several errors such as : Error in moran. mc(My_raster, lw.l, 999) : r is not a numeric vector So I converted "r" to a vector using as.vector() in : M Hello > Thanks a lot. > Indeed your explanation makes it much more clear. > Kind regards, > On Thu, Oct 21, 2021 at 11:40 AM Roger Bivand wrote: >> On Thu, 21 Oct 2021, Gabriel Cotlier wrote: >> > Hello >> > Thank you very much. >> > I have large raster layers and would like to ask, in order to reduce the >> > processing time of the simulation, choosing a smaller nsim value could >> > ? and if so, what could be the minimum nsim value recommended? >> If you examine the code, you see that raster::Moran() always performs >> multiple raster::focal() calls, effectively constructing the spatial >> weights on each run. The problem is not nsim, it is the way that >> raster::Moran() is constructed. >> Indeed, using a Hope-type test gives the same inferential outcome as >> using >> asymptotic methods for global tests for spatial autocorrelation, and is >> superfluous, but no other approach is feasible for raster::Moran() (which >> by the way only seems to use 4 neighbours although claiming it uses 8). >> 1. How large is large? >> It is possible to generate nb neighbour objects for large data sets, >> making the use of spdep::moran.test() feasible, certainly for tens of >> thousands of observations (all the census tracts in coterminous US, for >> example). >> This uses distances between raster cell centres to find neighbours: >> > library(spdep) >> Loading required package: sp >> Loading required package: spData >> Loading required package: sf >> Linking to GEOS 3.10.0, GDAL 3.3.2, PROJ 8.1.1 >> > crds > > dim(crds) >> [1] 96 2 >> > grd > > grd$z > > system.time(dnbr > user system elapsed >> 30.065 0.235 30.381 >> Neighbour list object: >> Number of regions: 96 >> Number of nonzero links: 3836000 >> Percentage nonzero weights: 0.0004162326 >> Average number of links: 3.995833 >> > system.time(dnbq > user system elapsed >> 54.502 0.080 54.699 >> Neighbour list object: >> Number of regions: 96 >> Number of nonzero links: 7668004 >> Percentage nonzero weights: 0.0008320317 >> Average number of links: 7.987504 >> Once the neighbour objects are ready, conversion to spatial weights >> objects takes some time, and computing I with a constant mean model >> depends on the numbers of neighbours: >> > system.time(lwr > user system elapsed >>6.694 0.000 6.722 >> > system.time(Ir > user system elapsed >> 24.371 0.000 24.470 >> Moran I test under randomisation >> data: grd$z >> weights: lwr >> Moran I statistic standard deviate = -0.06337, p-value = 0.5253 >> alternative hypothesis: greater >> sample estimates: >> Moran I statistic Expectation VarianceRe: [R-sig-Geo] question on raster Moran's I statistical significance
2021-10-21 Thread Gabriel Cotlier Hello Thanks a lot. Indeed your explanation makes it much more clear. Kind regards, On Thu, Oct 21, 2021 at 11:40 AM Roger Bivand wrote: > On Thu, 21 Oct 2021, Gabriel Cotlier wrote: > > Hello > > Thank you very much. > > I have large raster layers and would like to ask, in order to reduce the > > processing time of the simulation, choosing a smaller nsim value could > > ? and if so, what could be the minimum nsim value recommended? > If you examine the code, you see that raster::Moran() always performs > multiple raster::focal() calls, effectively constructing the spatial > weights on each run. The problem is not nsim, it is the way that > raster::Moran() is constructed. > Indeed, using a Hope-type test gives the same inferential outcome as using > asymptotic methods for global tests for spatial autocorrelation, and is > superfluous, but no other approach is feasible for raster::Moran() (which > by the way only seems to use 4 neighbours although claiming it uses 8). > 1. How large is large? > It is possible to generate nb neighbour objects for large data sets, > making the use of spdep::moran.test() feasible, certainly for tens of > thousands of observations (all the census tracts in coterminous US, for > example). > This uses distances between raster cell centres to find neighbours: > > library(spdep) > Loading required package: sp > Loading required package: spData > Loading required package: sf > Linking to GEOS 3.10.0, GDAL 3.3.2, PROJ 8.1.1 > > crds > dim(crds) > [1] 96 2 > > grd > grd$z > system.time(dnbr user system elapsed > 30.065 0.235 30.381 > Neighbour list object: > Number of regions: 96 > Number of nonzero links: 3836000 > Percentage nonzero weights: 0.0004162326 > Average number of links: 3.995833 > > system.time(dnbq user system elapsed > 54.502 0.080 54.699 > Neighbour list object: > Number of regions: 96 > Number of nonzero links: 7668004 > Percentage nonzero weights: 0.0008320317 > Average number of links: 7.987504 > Once the neighbour objects are ready, conversion to spatial weights > objects takes some time, and computing I with a constant mean model > depends on the numbers of neighbours: > > system.time(lwr user system elapsed >6.694 0.000 6.722 > > system.time(Ir user system elapsed > 24.371 0.000 24.470 > Moran I test under randomisation > data: grd$z > weights: lwr > Moran I statistic standard deviate = -0.06337, p-value = 0.5253 > alternative hypothesis: greater > sample estimates: > Moran I statistic Expectation Variance > -4.679899e-05 -1.041668e-06 5.213749e-07 > > system.time(lwq user system elapsed >6.804 0.000 6.828 > > system.time(Iq user system elapsed > 46.703 0.012 46.843 > Moran I test under randomisation > data: grd$z > weights: lwq > Moran I statistic standard deviate = -0.70373, p-value = 0.7592 > alternative hypothesis: greater > sample estimates: > Moran I statistic Expectation Variance > -3.604417e-04 -1.041668e-06 2.608222e-07 > 2. The larger your N, the less likely that the test means anything at all, > because the assumption is that the observed entities are not simply > arbitrary products of, say, resolution. > If you think of global Moran's I as a specification test of a regression > of the variable of interest on the constant (the mean model is just the > constant), for raster data the resolution controls the outcome > (downscaling/upscaling will shift Moran's I). If you include covariates, > patterning in the residuals of a richer model may well abate. > Hope this clarifies, > Roger > > Thanks a lot again. > > Kind regards, > > On Wed, Oct 20, 2021 at 8:45 PM Roger Bivand > wrote: > >> On Wed, 20 Oct 2021, Gabriel Cotlier wrote: > >>> Hello, > >>> I would like to estimate the Moran's *I* coefficient for raster data > >>> together with the statical significance of the spatial autocorrelation > >>> obtained. > >>> I found that the raster package function Moran() although calculates > >>> spatial autocorrelation index it apparently does not give directly the > >>> statical significance of the results obtained : > >>> https://search.r-project.org/CRAN/refmans/raster/html/autocor.html > >>> Could it be be possible to obtain the statistical significance of the > >>> results with either raster package or similar one? > >> fortunes::fortune("This is R") > >> library(raster) > >> r >> values(r) >> fRe: [R-sig-Geo] question on raster Moran's I statistical significance
2021-10-21 Thread Roger Bivand On Thu, 21 Oct 2021, Gabriel Cotlier wrote: Hello Thank you very much. I have large raster layers and would like to ask, in order to reduce the processing time of the simulation, choosing a smaller nsim value could help ? and if so, what could be the minimum nsim value recommended? If you examine the code, you see that raster::Moran() always performs multiple raster::focal() calls, effectively constructing the spatial weights on each run. The problem is not nsim, it is the way that raster::Moran() is constructed. Indeed, using a Hope-type test gives the same inferential outcome as using asymptotic methods for global tests for spatial autocorrelation, and is superfluous, but no other approach is feasible for raster::Moran() (which by the way only seems to use 4 neighbours although claiming it uses 8). 1. How large is large? It is possible to generate nb neighbour objects for large data sets, making the use of spdep::moran.test() feasible, certainly for tens of thousands of observations (all the census tracts in coterminous US, for example). This uses distances between raster cell centres to find neighbours: library(spdep) Loading required package: sp Loading required package: spData Loading required package: sf Linking to GEOS 3.10.0, GDAL 3.3.2, PROJ 8.1.1Re: [R-sig-Geo] question on raster Moran's I statistical significance
2021-10-20 Thread Gabriel Cotlier Hello Thank you very much. I have large raster layers and would like to ask, in order to reduce the processing time of the simulation, choosing a smaller nsim value could help ? and if so, what could be the minimum nsim value recommended? Thanks a lot again. Kind regards, On Wed, Oct 20, 2021 at 8:45 PM Roger Bivand wrote: > On Wed, 20 Oct 2021, Gabriel Cotlier wrote: > > Hello, > > I would like to estimate the Moran's *I* coefficient for raster data and > > together with the statical significance of the spatial autocorrelation > > obtained. > > I found that the raster package function Moran() although calculates the > > spatial autocorrelation index it apparently does not give directly the > > statical significance of the results obtained : > > https://search.r-project.org/CRAN/refmans/raster/html/autocor.html > > Could it be be possible to obtain the statistical significance of the > > results with either raster package or similar one? > fortunes::fortune("This is R") > library(raster) > r values(r) f (rI r1 nsim res set.seed(1) > for (i in 1:nsim) { >values(r1) res[i] } > Hope-type tests date back to Cliff and Ord; they are permutation > bootstraps. > r_g library(spdep) > nb set.seed(1) > o return_boot=TRUE) > x_a plot(density(o$t), xlim=x_a) > abline(v=o$t0) > lines(density(res), lty=2) > abline(v=rI, lty=2) > It is not immediately obvious from the code of raster::Moran() why it is > different, possibly because of padding the edges of the raster and thus > increasing the cell count. > For added speed, the bootstrap can be parallelized in both cases; polygon > boundaries are perhaps not ideal. > Hope this clarifies. Always provide a reproducible example, never post > mail. > Roger Bivand > > Thanks a lot. > > Kind regards, > > Gabriel > > [[alternative HTML version deleted]] > > R-sig-Geo mailing list > > [email protected] > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > Roger Bivand > Emeritus Professor > Department of Economics, Norwegian School of Economics, > Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway. > e-mail: [email protected] > https://orcid.org/-0003-2392-6140 > https://scholar.google.no/citations?user=AWeghB0J=en Gabriel Cotlier, PhD Haifa Research Center for Theoretical Physics and Astrophysics (HCTPA) University of Haifa Israel [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo I would like to estimate the Moran's *I* coefficient for raster data and together with the statical significance of the spatial autocorrelation obtained. I found that the raster package function Moran() although calculates the spatial autocorrelation index it apparently does not give directly the statical significance of the results obtained : https://search.r-project.org/CRAN/refmans/raster/html/autocor.html Could it be be possible to obtain the statistical significance of the results with either raster package or similar one? fortunes::fortune("This is R") library(raster)Re: [R-sig-Geo] question on raster Moran's I statistical significance
2021-10-20 Thread Roger Bivand On Wed, 20 Oct 2021, Gabriel Cotlier wrote: Hello, I would like to estimate the Moran's *I* coefficient for raster data and together with the statical significance of the spatial autocorrelation obtained. I found that the raster package function Moran() although calculates the spatial autocorrelation index it apparently does not give directly the statical significance of the results obtained : https://search.r-project.org/CRAN/refmans/raster/html/autocor.html Could it be be possible to obtain the statistical significance of the results with either raster package or similar one? fortunes::fortune("This is R") library(raster)Re: [R-sig-Geo] question on the use of Leaflet for plotting points and rater GeoTIFF
2021-08-13 Thread Gabriel Cotlier Dear Ben, Thanks a lot for your help with the useful and appropriate link you sent me. There was the solution for the flipping of the colorbar. I could finally flip the color bar and it worked out correctly. I also found a little bit less opaque color palette, that is the matlab "jet" color in the R "matlab" package, but it is still a bit opaque for my personal taste. I will give it some trials and errors with many different opacities to see what happened. Thanks a lot again for your helpful advice and guidance, Kind regards, Gabriel > I don't think we can run the code since it isn't reproducible (see > https://CRAN.R-project.org/package=reprex >Re: [R-sig-Geo] question on the use of Leaflet for plotting points and rater GeoTIFF
2021-08-13 Thread Ben Tupper I don't think we can run the code since it isn't reproducible (see https://CRAN.R-project.org/package=reprexRe: [R-sig-Geo] Question spdep package - aliased variables and NA z-values
2019-05-28 Thread Roger Bivand On Mon, 27 May 2019, Raphael Mesaric via R-sig-Geo wrote: Dear all, My model consists of observations for all four seasons, and each of the seasons is analyzed with respect to five different time slots. All of these observations are referring to the same grid, which means that my neighborhood matrix is repeated twenty times (I took the neighborhood matrix of the original grid and modified the entries in such a way that I have a diagonal matrix where each diagonal block entry consists of the original neighborhood object, with adjusted values). I recently run my model with the errorsarlm() and lagsarlm() functions. It worked fine, and the results coincide well with the linear regression results, which is great. However, there are a few other strange things: These are not really appropriate models, maybe look at the splm package for panel/longitudinal equivalents. To debug the aliasing in the model fitting functions in spatialreg (no longer spdep), please provide a reproducible example, preferably from a built-in data set. Roger In the lag-model, I have only estimates for some variables, but the z-value is indicated as NA (for the error model it works fine). After running the SAR models, R tells me that the two dummy variables for the seasons are aliased. However, in the linear regression I did not have any problems, and the variables are not correlated. I think the second issue might be related to the fact that the number of observations seems to correspond to the number of observations for one season, instead of all four seasons. In other words, it looks like R has realized that the neighborhood repeats itself for each season. Interestingly enough, R has not realized that the pattern repeats itself inside the seasons as well for the time slots. Or is there another possible reason for the fact that the number of observations corresponds to only one fourth of all the observations actually included? Has anybody experienced these issues as well? Any idea on how to resolve these problems would be greatly appreciated. Thank you in advance! Best regards, Raphael Mesaric [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Roger Bivand Department of Economics, Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; e-mail: [email protected] https://orcid.org/-0003-2392-6140 https://scholar.google.no/citations?user=AWeghB0J=en R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Thank you very much for your reply. My thesis supervisor told me to use GWR to explore small-scale differences and test the consistency of the SAR models. It happens. From your problem description, I'm unsure whether GWR or SAR are appropriate methods, and would suggest the mixed-model and multilevel literatures. Mapable IID and spatially structured random effects may be more appropriate. Roger So, I tried to fit the GWR model with all my observations. However, if I understood you correctly, I would have to choose smaller sections of the grid and try to fit a GWR there. What I am not sure about how to do this in order to get informative results. As for the SAR model, the dependent variable consists n blocks where the m-th entry of each block corresponds to the m-th cell of my grid. So, if I have to reduce the grid to about 5K entries, I would need to take a subset of roughly 250 cells and then take all the 20 observations for each of these cells? And then try to fit a GWR there? Best, Raphael Am 17.05.2019 um 14:20 schrieb Roger Bivand : On Fri, 17 May 2019, Raphael Mesaric via R-sig-Geo wrote: Dear all, Is there an option to shorten the running time for the gwr() function, similar to the ‚LU‘ method for lagsarlm() in the spdep package? Because I have a model with roughly 500’000 observations, and the running time at the moment is quite long, respectively it has not yet terminated. What are you actually doing? Why did you choose GWR? Are you fitting a GWR with 500K observations, or have you fitted a GWR with many fewer observations, and are now rendering that fitted model with 500K fit points? GWR is only for detecting possible non-stationarity or similar mis-specification in moderately sized data sets. Trying to fit with 500K gives a dense hat matrix of 500K x 500K, which is imposssible (or were it possible would be uninformative). Think of 5K as a sensible maximum if GWR is condidered sensible at all. I would think that finding a bandwidth is impossible too. Roger Thank you for your help in advance. Best regards, Raphael Mesaric R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Roger Bivand Department of Economics, Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; e-mail: [email protected] https://orcid.org/-0003-2392-6140 https://scholar.google.no/citations?user=AWeghB0J=en Roger Bivand Department of Economics, Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; e-mail: [email protected] https://orcid.org/-0003-2392-6140 https://scholar.google.no/citations?user=AWeghB0J=en___ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Look at multilevel models, fitting iid and/or mrf to the three regular dimensions. Look at the mixed models literature. Some references in https://doi.org/10.1016/j.spasta.2017.01.002, but there is much more. Hope this helps, Roger Roger Bivand Norwegian School of Economics Bergen, Norway Fra: Raphael Mesaric Sendt: fredag 17. mai, 14.58 Emne: Re: [R-sig-Geo] Question spdep package - extending nb object Til: Roger Bivand Kopi: [email protected] Dear Roger, Thank you very much for your reply. It is kind of a spatial panel model, but I have only one time series. To describe it differently: The dependent variable consists of twenty blocks with averaged data. Each block is related to one season and one out of five time slots (that is why I have 20 blocks - 4 seasons times 5 time slots). And each block of data refers to the very same grid, i.e. the first entry of the first block corresponds to the same call as the first entry of all the other blocks, the second entries correspond to the second cell etc. The grid is of hectare resolution, i.e. 100 x 100m (one could see this as similar to post codes). Therefore, nb2blocknb() sounds interesting, even though I do not really have „location-less observations“, but just observations with different attributes (season and time). Would it work if my ‚ID‘ argument is built in such a way that it repeats 19 times the identification tags of the first block? For example: rep(seq(from = 1, to = m, by = 1),19) , where m is the number of cells in my grid? Best, Raphael Am 17.05.2019 um 14:13 schrieb Roger Bivand mailto:[email protected]>>: On Thu, 16 May 2019, Raphael Mesaric via R-sig-Geo wrote: Dear all, I am working on a spatial regression model which is built in such a way that the dependent variable consists of n blocks of multiple entries, and each block is referring to the same spatial grid. In order to run the lagsarlm() function or the errsarlm() function, I now need to extend my initial nb object from one grid to n grids. In other words, I need the nb object to be repeated n times. However, so far I did not find a way to do this properly. I tried to just replicate the nb object by repeating the entries of the neighbours list and the weights list, respectively. dnn800_3w$neighbours <- rep(dnn800_3w$neighbours,20) dnn800_3w$weights <- rep(dnn800_3w$weights,20) When I do this, I get the following error message: Error in is.symmetric.nb(listw$neighbours, FALSE) : Not neighbours list Which makes sense, as the matrix is not symmetric any more due to the steps described above. However, I cannot think of a correct implementation at the moment. Not even worth trying, as you found. Do you mean a Kronecker product, as found in spatial panel models? Or do you mean that you have blocks of observations without position that belong to upper-level objects with known position (postcodes or similar)? If the latter, look at ?spdep::nb2blocknb. Hope this clarifies, Roger Can you help me with this issue? Any suggestions are highly appreciated. Thank you in advance. Best regards, Raphael Mesaric [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected]<mailto:[email protected]> https://stat.ethz.ch/mailman/listinfo/r-sig-geo Roger Bivand Department of Economics, Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; e-mail: [email protected]<mailto:[email protected]> https://orcid.org/-0003-2392-6140 https://scholar.google.no/citations?user=AWeghB0J=en Thank you very much for your reply. It is kind of a spatial panel model, but I have only one time series. To describe it differently: The dependent variable consists of twenty blocks with averaged data. Each block is related to one season and one out of five time slots (that is why I have 20 blocks - 4 seasons times 5 time slots). And each block of data refers to the very same grid, i.e. the first entry of the first block corresponds to the same call as the first entry of all the other blocks, the second entries correspond to the second cell etc. The grid is of hectare resolution, i.e. 100 x 100m (one could see this as similar to post codes). Therefore, nb2blocknb() sounds interesting, even though I do not really have „location-less observations“, but just observations with different attributes (season and time). Would it work if my ‚ID‘ argument is built in such a way that it repeats 19 times the identification tags of the first block? For example: rep(seq(from = 1, to = m, by = 1),19) , where m is the number of cells in my grid? Best, Raphael > Am 17.05.2019 um 14:13 schrieb Roger Bivand : > On Thu, 16 May 2019, Raphael Mesaric via R-sig-Geo wrote: >> Dear all, >> I am working on a spatial regression model which is built in such a way that >> the dependent variable consists of n blocks of multiple entries, and each >> block is referring to the same spatial grid. In order to run the lagsarlm() >> function or the errsarlm() function, I now need to extend my initial nb >> object from one grid to n grids. In other words, I need the nb object to be >> repeated n times. However, so far I did not find a way to do this properly. >> I tried to just replicate the nb object by repeating the entries of the >> neighbours list and the weights list, respectively. >> dnn800_3w$neighbours > dnn800_3w$weights > >> When I do this, I get the following error message: >> Error in is.symmetric.nb(listw$neighbours, FALSE) : Not neighbours list >> Which makes sense, as the matrix is not symmetric any more due to the steps >> described above. However, I cannot think of a correct implementation at the >> moment. > Not even worth trying, as you found. Do you mean a Kronecker product, as > found in spatial panel models? Or do you mean that you have blocks of > observations without position that belong to upper-level objects with known > position (postcodes or similar)? If the latter, look at ?spdep::nb2blocknb. > Hope this clarifies, > Roger >> Can you help me with this issue? Any suggestions are highly appreciated. >> Thank you in advance. >> Best regards, >> Raphael Mesaric >> [[alternative HTML version deleted]] >> R-sig-Geo mailing list >> [email protected]Re: [R-sig-Geo] Question spdep package - lagsarlm not terminating
2019-05-06 Thread Raphael Mesaric via R-sig-Geo Dear Roger, Thank you again for your fast response and the helpful answers and clarifications. I will try to resolve all the issues with your recommendations in the coming days. Just a few comments here: I created the listw object by nb2listw() (without any error messages) and set „zero.policy=TRUE“ while calling the function. However, I used the function subset() instead of subset.nb(). I will play around with it a bit and see which works best for me. But it is good to know that observations with missing values are dropped by default and the matrix adjusted. With regards to your question related to GWR: Yes, the error occurs while printing a fitted object. And I also thought that it might be linked to the collinearity in the covariates, just as it is the case for the error with lagsarlm(). However, if you do not recommend it for any other purposes I will have to discuss its use with my thesis supervisor again. Have a nice day! Best regards, Raphael > Am 06.05.2019 um 17:20 schrieb Roger Bivand : > On Mon, 6 May 2019, Raphael Mesaric via R-sig-Geo wrote: >>> Dear Roger, >>> Thank you very much for your fast response. >>> My weight matrix stems from a rectangular (but not square) grid. So, I >>> think I will have to use the „LU“ method. > No, this is a misunderstanding. All weights matrices used for fitting models > are square nxn matrices by definition. If your grid had r rows and c columns, > n = r * c, and the weights matrix has n rows and n columns. >>> However, by doing that several other error messages popped up. The first >>> one was the following: >>> Error in spatialreg::errorsarlm(formula = formula, data = data, listw = >>> listw, : >>> NAs in lagged dependent variable >>> In addition: Warning messages: >>> 1: Function errorsarlm moved to the spatialreg package >>> 2: In lag.listw(listw, y, zero.policy = zero.policy) : >>> NAs in lagged values >>> And this happened even though I do actually not have NaN values in my >>> dataset (I removed them in advance). When I completely reload all >>> variables, it works sometimes for a few tries but after a while the error >>> messages pops up again. Do you know why this might happen? And could I >>> theoretically use a dataset with NaN values and the function omits them? > By default, observations with missing values are dropped and the weights > object is subsetted to match. This may produce no-neighbour observations. You > can choose to set their spatially lagged values to zero by using > zero.policy=TRUE, which is the probable cause of your error. It could very > well be that your weights object itself contains no-neighbour observations if > not created using spdep::nb2listw() - which will alert you to the problem. >>> The second error message was the following: >>> Error in solve.default(-(mat), tol.solve = tol.solve) : >>> system is computationally singular: reciprocal condition number = >>> 2.71727e-26 >>> I found a reference online (related to glm) that this might be linked to >>> explanatory variables which have a high correlation. Would eliminating some >>> of the variables do the job here as well? > Please do not use references online (especially if you do not link to them). > The information you need is in the help page, referring to the tol.solve= > argument. Indeed, your covariates are scaled such that their are either > colinear (unlikely), or that the numerical values of the coefficient > variances are very different from the spatial coefficient. You may need to > re-scale the response or covariates in order to invert the matrix needed to > yield the variances of the coefficients. >>> Then I have two other questions, where I couldn’t find any answers online: >>> Is there an option to eliminate specific rows of a nb object? I tried the >>> following command: >>> rook234x259_2 > spdep::subset.nb() only. Your suggestion only creates a big mess, as > subsetting reduces the number of rows, and neighbours are indexed between 1 > and n. Consequently, at least some of the remaining vectors will point to > neighbours with ids > n, and many of the others will point to the wrong > neighbours. You can convert to a sparse matrix, subset using "[", and then > convert back, but the subset method should work. >>> where index is a vector with the row numbers I would like to keep. But this >>> converts the nb to a list object, and there is no list2nb function (as much >>> as I am aware of). And the option „subset“ does not only remove the rows, >>> but also the corresponding cells which leads to a different neighborhood >>> matrix. > I do not know what you mean, of course it has to remove links in both > directions. >>> The last question refers to the package spgwr. I would like to run this >>> model as well. The gwr function runs successfully, but when callingRe: [R-sig-Geo] Question spdep package - lagsarlm not terminating
2019-05-05 Thread Roger Bivand On Sat, 4 May 2019, Raphael Mesaric via R-sig-Geo wrote: Dear all, I have a question with regards to the function „lagsarlm" from the package spdep. My problem is that the function is not terminating. Of course, I have quite a big grid (depending on the selection either 34000 or 60600 entries) and I also have a lot of explanatory variables (about 40). But I am still wondering whether there is something wrong. If you read the help page (the function is in the spatialreg package and will be dropped from spdep shortly), and look at the references (Bivand et al. 2013), you will see that the default value if the method= argument is "eigen". For small numbers of observations, solving the eigenproblem of a dense weights object is not a problem, but becomes demanding on memory as n increases. You are using virual memory (and may run out of that too) which makes your machine unresponsive. If you choose an alternative method, typically "Matrix" for symmetric or similar to symmetric sparse weights, or "LU" for asymmetric sparse weights. data(house, package="spData") dim(house) [1] 2535724 LO_nb Neighbour list object: Number of regions: 25357 Number of nonzero links: 74874 Percentage nonzero weights: 0.01164489 Average number of links: 2.952794Re: [R-sig-Geo] Question about HSAR package
2018-09-27 Thread Justin Schon Ok, it's helpful to know that I need to zoom in on those three things. I created the Random Effects Matrix by hand in Excel, so I read that into R and put the matrix into the format recommended here: https://cran.r-project.org/web/packages/HSAR/vignettes/hsar.html Below, I show how I made W, M, and Delta.mat, before I try to estimate the model. Hopefully this helps. constit HSAR.model1 This code tells nothing, the problem is in your construction of W, M > and/or Delta. Pleaseng show this code too, best as a reproducible example. > Tip: sometimes running traceback() after an error shows where it happens. > Roger Bivand > Norwegian School of Economics > Bergen, Norway > Fra: Justin Schon > Sendt: torsdag 27. september, 21.36 > Emne: [R-sig-Geo] Question about HSAR package > Til: [email protected] > Dear all, I am receiving the error "not an S4 object" when I attempt to > estimate the hierarchal spatial auto-regressive model from the HSAR > package. I have attempted several ways of creating the lower level matrix > and higher level matrix. Rather than asking if members of this list can > help with the code, I am first wondering if anyone can explain why this > error would appear. I am including the code that estimates the model, as > well as the error, below: > HSAR.model1 Justin Schon Post-Doctoral Researcher on Environmental Change and Migration MURI Migration Research Team: http://murimigration.org/ University of Florida Fellow, Initiative for Sustainable Energy Policy (ISEP) [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo This code tells nothing, the problem is in your construction of W, M and/or Delta. Pleaseng show this code too, best as a reproducible example. Tip: sometimes running traceback() after an error shows where it happens. Roger Bivand Norwegian School of Economics Bergen, Norway Fra: Justin Schon Sendt: torsdag 27. september, 21.36 Emne: [R-sig-Geo] Question about HSAR package Til: [email protected] Dear all, I am receiving the error "not an S4 object" when I attempt to estimate the hierarchal spatial auto-regressive model from the HSAR package. I have attempted several ways of creating the lower level matrix and higher level matrix. Rather than asking if members of this list can help with the code, I am first wondering if anyone can explain why this error would appear. I am including the code that estimates the model, as well as the error, below: > HSAR.model1 [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo I've posted this to the mixed models forum as well but thought people here might know: I want to use the the correlation setting with corSpher in nlme to account for potential spatial autocorrelation in my data. My data include observations from across the globe with locations in latitude and longitude (decimal degrees). From the R help for corSpher, the syntax for their wheat example would be something like: Nothing in the documentation nor the underlying book suggests that the coordinates are anything other than planar, so geographical coordinates in decimal degrees should not be used. Unfortunately, I do not have access to the underlying article in crop science, but there is no reason to think that the agricultural experiment extended from France to Turkey and France to south of Ghana. The documentation in nlme::corSpher() refers to stats::dist, and I can't see any way to insert fields::rdist.earth() in its place. So you can safely accept that the coefficients must be projected (i.e. planar). This is all accessible in the documentation and Pinheiro & Bates (2000). It is pretty certain that they did not envisage global scope in geographical coordinates, so maybe you should be using fields or gstat which are more likely to be able to handle Great Circle distances. Hope this clarifies, Roger fm1Wheat2Re: [R-sig-Geo] question about rgeos "gDistance" parameters
2017-11-10 Thread Roger Bivand Please post plain text only - since you didn't provide an example, your code wasn't deformed by HTML formatting, but had you provided an example it would have been. Did you look at the examples on the help page? library(rgeos) pt1 = readWKT("POINT(0.5 0.5)") pt2 = readWKT("POINT(2 2)") p1 = readWKT("POLYGON((0 0,1 0,1 1,0 1,0 0))") p2 = readWKT("POLYGON((2 0,3 1,4 0,2 0))") gDistance(pt1,pt2) gDistance(p1,pt1) gDistance(p1,pt2) gDistance(p1,p2) Modifying this:Re: [R-sig-Geo] Question about the anisotropy factor in the sum metric spatiotemporalcovariance gstat
2017-10-04 Thread Dr . Benedikt Gräler Dear Sara, the anisotropy in the metric spatio-temporal covariance models describes how spatial distances relate to temporal distance (i.e. to answer the question: At which temporal separation (x hours/days/weeks/...) do have two locations 1 m/km/... apart the same correlation?). Every model that includes the metric model has to have a spatio-temporal anisotropy parameter. For more details, please refer to the vignette of gstat [1]. https://cran.r-project.org/web/packages/gstat/vignettes/spatio-temporal-kriging.pdf Dear all. I have a question about the sum metric spatiotemporalcovariance function in gstat. I need to understand the theoretical background and the meaning of the anisotropy (Anis) in it. If anyone can help me in this I will be deeply appreciated Best regards [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Dr. Benedikt Gräler 52°North Initiative for Geospatial Open Source Software GmbH Martin-Luther-King-Weg 24 48155 Muenster, Germany E-Mail: [email protected] Fon: +49-(0)-251/396371-39 Fax: +49-(0)-251/396371-11 http://52north.org/ Twitter: @FiveTwoN General Managers: Dr. Albert Remke, Dr. Andreas Wytzisk Local Court Muenster HRB 10849 R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geoRe: [R-sig-Geo] Question niche based model whith R
2017-04-14 Thread Mathieu Basille You can also have a look at exploratory niche analyses, such as the Ecological niche factor analysis. A few of them are available in the 'adehabitatHS' package: https://cran.r-project.org/package=adehabitatHS Detailed explanations are available in the vignette (section 3: Design I studies): https://cran.r-project.org/web/packages/adehabitatHS/vignettes/adehabitatHS.pdf Hope this helps, Mathieu. On 04/10/2017 05:12 AM, Soufianou Abou via R-sig-Geo wrote: > Hello dears,I work on the modeling of the ecological niche of striga, I have > presence data (data) and bioclimatic variables and other variables > around.Indeed, my question how should I proceed to model the ecological niche > with software R? Pending receipt of my request, please accept my best > regards.SADDA Abou-Soufianou -- > DoctorantUniversité Dan Dicko Dankoulodo de Maradi,120 avenue Maman > Koraou,ADS, MaradiBP 465 Maradi/Niger E-mail: [email protected] : > (+227)96269987 > [[alternative HTML version deleted]] > R-sig-Geo mailing list > [email protected] > https://stat.ethz.ch/mailman/listinfo/r-sig-geo Mathieu Basille [email protected] | http://ase-research.org/basille +1 954-577-6314 | University of Florida FLREC « Le tout est de tout dire, et je manque de mots Et je manque de temps, et je manque d'audace. » — Paul Éluard R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geoRe: [R-sig-Geo] Question niche based model whith R
2017-04-10 Thread Babak Naimi There are several packages available in R to do ecological niche modelling. One easy-to-use option would be the sdm package. It is a comprehensive modelling and simulation framework to fit species distribution models supports several algorithms, and also is user-friendly. You can find a tutorial for a quick start at: http://biogeoinformatics.org Good luck, Babak > Date: Mon, 10 Apr 2017 09:12:54 + (UTC) > From: Soufianou AbouRe: [R-sig-Geo] Question impacts lagged independent variables lagsarlm
2016-06-05 Thread Jorge Cárcamo Prof. Bivand, Many thanks for your clarification. I think that I should read Läpple & Kelley again. One more question, I assume that if I create independent lag variables with: data$w_xxx wrote: > On Sat, 4 Jun 2016, Jorge Cárcamo wrote: > Dear all, >> I am working with a lagsarlm "mixed" model I executed: >> library(spdep) >> library(coda) >> dsts15 > idw15 > sy15 > mod.sdm.15> saa >> + owue + yearrain + yearfdi, data=data, listw=sy15, type="mixed", >> tol.solve=1.0e-12) >> summary(mod.sdm.15) >> And got following results (abridged table): >> Estimate Std. Error z value Pr(>|z|) >> (Intercept) 1.4255019 0.4824828 2.9545 0.003132 >> SD46 0.0149057 0.0116146 1.2834 0.199365 >> Totaland -0.0217556 0.0072475 -3.0018 0.002684 >> yearrain 0.0557373 0.0485612 1.1478 0.251062 >> yearfdi -0.0109979 0.0069536 -1.5816 0.113741 >> lag.SD46 0.0262273 0.0278045 0.9433 0.345539 >> lag.Totaland 0.0060992 0.0141515 0.4310 0.666474 >> lag.yearrain -0.2083118 0.0772751 -2.6957 0.007024 >> lag.yearfdi 0.0111460 0.0081562 1.3666 0.171761 >> Rho: -0.54702, LR test value: 10.052, p-value: 0.0015215 >> Immediatly, I exectued impacts(mod.sdm.15, R=1000), to get the impacts: >> Direct Indirect Total >> SD46 0.013136761 0.013451739 0.0265884995 >> Totaland -0.023517332 0.013397014 -0.0101203177 >> yearrain 0.079110119 -0.177734631 -0.0986245121 >> yearfdi -0.012677636 0.012773412 0.957761 >> Through simulation I manage to get credible intervals for direct, indirect >> and total impacts (HPDinterval(impacts, choice="XXX")). > As you must know from the references on the help page for impacts methods, > these are the combined impacts of the variables: > S_r(W) = (I - \rho W)^{-1} (\beta_r I - \gamma_r W) > where the direct impacts are sum(S_r(W))/n, etc. The \gamma_r are the > coefficients on W x_r. >> However, how can I get the impacts of the lagged variables? I have seen >> some publications such as: Läpple, D., & Kelley, H. (2014). Spatial >> dependence in the adoption of organic drystock farming in Ireland. That >> report a posterior mean and credible intervals for the lagged variables. > Given the above, either you are misreading Läpple & Kelley (I do not have > access), or both you and they are wrong. There are by definition on > separable impacts for the lagged X variables. > Hope this clarifies, > Roger >> All suggestions are welcome. >> Best regards, >> *Ing. Jorge Alfredo Cárcamo, M. Sc., Ph. D. (c)* >> Agriculture economics >> Georg-August-Universität Göttingen >> [[alternative HTML version deleted]] >> R-sig-Geo mailing list >> [email protected] >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > Roger Bivand > Department of Economics, Norwegian School of Economics, > Helleveien 30, N-5045 Bergen, Norway. > voice: +47 55 95 93 55; fax +47 55 95 91 00 > e-mail: [email protected] > http://orcid.org/-0003-2392-6140 > https://scholar.google.no/citations?user=AWeghB0J=en > http://depsy.org/person/434412 [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geoRe: [R-sig-Geo] Question impacts lagged independent variables lagsarlm
2016-06-05 Thread Roger Bivand On Sat, 4 Jun 2016, Jorge Cárcamo wrote: Dear all, I am working with a lagsarlm "mixed" model I executed: library(spdep) library(coda) dsts15 |z|) (Intercept) 1.4255019 0.4824828 2.9545 0.003132 SD46 0.0149057 0.0116146 1.2834 0.199365 Totaland -0.0217556 0.0072475 -3.0018 0.002684 yearrain 0.0557373 0.0485612 1.1478 0.251062 yearfdi -0.0109979 0.0069536 -1.5816 0.113741 lag.SD46 0.0262273 0.0278045 0.9433 0.345539 lag.Totaland 0.0060992 0.0141515 0.4310 0.666474 lag.yearrain -0.2083118 0.0772751 -2.6957 0.007024 lag.yearfdi 0.0111460 0.0081562 1.3666 0.171761 Rho: -0.54702, LR test value: 10.052, p-value: 0.0015215 Immediatly, I exectued impacts(mod.sdm.15, R=1000), to get the impacts: Direct Indirect Total SD46 0.013136761 0.013451739 0.0265884995 Totaland -0.023517332 0.013397014 -0.0101203177 yearrain 0.079110119 -0.177734631 -0.0986245121 yearfdi -0.012677636 0.012773412 0.957761 Through simulation I manage to get credible intervals for direct, indirect and total impacts (HPDinterval(impacts, choice="XXX")). As you must know from the references on the help page for impacts methods, these are the combined impacts of the variables: S_r(W) = (I - \rho W)^{-1} (\beta_r I - \gamma_r W) where the direct impacts are sum(S_r(W))/n, etc. The \gamma_r are the coefficients on W x_r. However, how can I get the impacts of the lagged variables? I have seen some publications such as: Läpple, D., & Kelley, H. (2014). Spatial dependence in the adoption of organic drystock farming in Ireland. That report a posterior mean and credible intervals for the lagged variables. Given the above, either you are misreading Läpple & Kelley (I do not have access), or both you and they are wrong. There are by definition on separable impacts for the lagged X variables. Hope this clarifies, Roger All suggestions are welcome. Best regards, *Ing. Jorge Alfredo Cárcamo, M. Sc., Ph. D. (c)* Agriculture economics Georg-August-Universität Göttingen [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Roger Bivand Department of Economics, Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 91 00 e-mail: [email protected] http://orcid.org/-0003-2392-6140 https://scholar.google.no/citations?user=AWeghB0J=en http://depsy.org/person/434412___ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geoRe: [R-sig-Geo] question about nb2listw row standardization
2015-08-04 Thread Roger Bivand On Tue, 4 Aug 2015, Mueller,Drew wrote: I was hoping for a clarification on the row standardization of nb2listw. I have implemented a nb list based on distance, and a function dlist - nbdists(col.gal.nb, coords) dlist - lapply(dlist, function(x) 1/x) my question is, if I call Wa - nb2listw(x, glist = dlist, style = W) Wb - nb2listw(x, glist = dlist, stlye = B) is Wa then row standardized based on inverse distance or does the glist override this and I get actual undstandardized inverse distance weights? Also, would Wb actually be binary, or does the glist argument override this to make the weights 1/x instead of 1's and 0's? If my goal is to maintain the actual inverse distance weight without standardization (for better economic interpretation), does Wb give me an unstandardized W, or do I need to call a different function to avoid row standardization? If my goal IS to row standardize a weights matrix with weights based on inverse distance, then is Wa the correct function call? Mainly I was having trouble interpreting the help file for nb2listw as to whether style = B always creates 1's and 0's regardless of what glist is set to. Thanks Please do check: library(spdep) example(columbus) coords - coordinates(columbus) dlist - nbdists(col.gal.nb, coords) dlist - lapply(dlist, function(x) 1/x) Wa - nb2listw(col.gal.nb, glist = dlist, style = W) Wb - nb2listw(col.gal.nb, glist = dlist, style = B) summary(sapply(Wa$weights, sum)) summary(sapply(Wb$weights, sum)) W row standardises the general weights, B passes through general weights. IDW are perverse, and without - contrary to your assertion - economic meaning, because the scaling depends closely on the metric of the coordinates: summary(unlist(dlist)) are all 1. Roger [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Roger Bivand Department of Economics, Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 91 00 e-mail: [email protected] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo # read external grid with rgdal; for example this data on package sp x - readGDAL(system.file(external/test.ag, package=sp)[1]) #transform to im format of spatstat x.im- as(x, im) plot(x.im) # you should transform the values to correspond to the desired density of hospitals with arithmetic operations, for example: x.im2- x.im/10e7 plot(x.im2) # And now you rpoispp to simulate random locations of hospitals an use this to compute envelopes rpoispp(x.im2) points(rpoispp(x.im2)) MArcelino El 05/05/2015 a las 19:37, eugeneby escribió: Good afternoon everyone! I'm new to this forum and was hoping to get some clarity about running k-function analysis in R. I would like to see whether hospitals in a certain state are clustered, dispersed or randomly distributed. I could just use the regular Ripley's k-function analysis to examine the spatial distribution of the hospitals and simulate confidence envelopes using the envelope command. However, it won't be surprising to find that hospitals are clustered, as more will be located in areas where the population (or population density) is higher. So my question is whether there any way for R to take into consideration a raster file showing the different population densities in various parts of my study region when simulating confidence envelopes? Is this what the Kinhom function does? And if so, how do I get it to read my raster file? Thank you! View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Question-about-inhomogeneous-k-functions-tp7588180.html Sent from the R-sig-geo mailing list archive at Nabble.com. R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Thank you so much for your quick and helpful reply! If I may, a couple follow up questions: 1. Instead of the rpoispp(x.im2), should I use rpoint(260,x.im2,win=BoundaryPolygonsOW) to generate random point patterns which contain the same number of points as in the original pattern (i.e., 260)? Because the rpoispp command generates a ridiculously large number of points and it doesn't seem that there's a simple way to indicate the number of points that should be generated. 2. Should I then manually simply use the Kest command to calculate the k-function for each of the, say, 99 simulations and compare my original pattern with the min and max of these k-functions at different distances? Or is there code which does that automatically? E.g., Envelopes - spatstat::envelope(pp, Kinhom, lambda=x.im2, nsim=99, correction=Ripley)? (where pp is my original point pattern). Thanks so much! All the best, ~Eugene View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Question-about-inhomogeneous-k-functions-tp7588180p7588183.html Sent from the R-sig-geo mailing list archive at Nabble.com. R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geoRe: [R-sig-Geo] Question about derivative work - what is the license for map derived using e.g. spatial predict function?
2014-11-27 Thread Jochen Albrecht What makes Tomislav's question really hard to answer is that there is no *one* copyright law. Just look at the long list of exceptions in the Google legal notices that he provided. My experience is only with the United States and here the matter is relatively clear (in spite of Hollywood's extreme provisions). 1. As mentioned, facts are not copyrightable; so if you have just plain data without any create packaging/mapping/analysis then these are not protected - regardless of what the vendor may claim (and some vendor are very creative with their claims). 2. Your derivation is your work and does therefore not inherit any copyright either- even if the original was copyrighted, say because of the packaging. The Google Earth digitization is a case here. Possibly not in Israel, but here in the U.S. you may freely digitize , i.e., create new data - the boundaries that you are digitizing are not part of the source material but your (creative) interpretation. You may, however, not copy; that is, you may not re-assemble the tiles. 3. The notions of original vs. derivative are every once in a while contested in court but so far, the courts have usually sided with the defense. There are only two authors in our field that have published on this topic that I am aware off. Check out either Onsrud or Cho for a long list of cases. Unless you find a way to make millions of your work (which would interest me immensely) you are unlikely to be prosecuted across national boundaries because the claimant's effort would not be worth it. Cheers, Jochen Dr. Jochen Albrecht Computational and Theoretical Geography HN 1030 Hunter College CUNY 695 Park Avenue New York, NY 10065 On Thu, Nov 27, 2014 at 4:50 PM, Tomislav Hengl [email protected] wrote: Hi Barry, Thanks for reminding me about the you can not copyright facts principle ( http://www.newmediarights.org/business_models/artist/are_facts_copyrighted). So the point measured values can be considered to be, in fact, - facts. So here is one more time the question to everyone (I am sure some of you have experienced these issues in practice): 1. A point data set has a restrictive license. 2. From this point data you can derive predictions - these predictions can not be used to re-produce point locations and values, hence the source is untraceable. 3. Does this derivative work (spatial prediction) has to follow the same data license as used for the point data, or this a new data sets for which you can set license freely. Some basic principles to consider: 1. One can not copyright facts (http://www.newmediarights. org/business_models/artist/are_facts_copyrighted) but facts can be protected using e.g. the trade secret laws. 2. If one would digitize content from Google Earth / maps imagery i.e. digitize polygon maps representing geomorphology without a permission from Google this would mean breaking of the Legal Notices ( http://www.google.com/help/legalnotices_maps.html) - has anyone bee prosecuted for this ever? 3. The transformation, modification or adaptation of the work must be substantial and bear its author's personality to be original and thus protected by copyright (http://en.wikipedia.org/wiki/Derivative_work) hence it is OK to create derivative works as long as one can prove that there are substantial new elements. Do spatial predictions from points fall into this category? thank you! T. Hengl On 27-11-2014 16:38, Barry Rowlingson wrote: On Thu, Nov 27, 2014 at 2:44 PM, Edzer Pebesma [email protected] wrote: Tom, in your example below, x contains the kriging variance; points with zero kriging variance must be observation locations, with the predicted value equal to the observation. In case the nugget in m would have been replaced by a measurement error component (Err in m and m1 below), you would not have this effect, and also have no discontinuity in the interpolated surface at observation locations: Now all that is going to be fun to explain to a judge and jury. Suppose you took all the notes of a Bach fugue as X=time, Y=pitch, and interpolated them in time, then created a new piece using the interpolation prediction at time points between the notes, would this be a derivative work? Yes, I think: A “derivative work” is a work based upon one or more preexisting works, such as a translation, musical arrangement, dramatization, fictionalization, motion picture version, sound recording, art reproduction, abridgment, condensation, or any other form in which a work may be recast, transformed, or adapted. [Wikipedia, where I get all my legal advice from] - they key word being transformed. But does a dataset count as a work here? It is supposedly a set of measurements of the truth, rather than something that is a creative work. A random australian web site that comes up tops in a google search says: A dataset may attract copyright protection (as aRe: [R-sig-Geo] Question on histograms from Raster and RasterVIS packages
2014-11-19 Thread Pascal Oettli Here is an example, with the mean in dotted red line and the median in dotted blue line: library(rasterVis) f - system.file(external/test.grd, package=raster) r - raster(f) mn - cellStats(r, mean) md - cellStats(r, median) histogram(r) + layer(panel.abline(v = mn, col = red, lty = 2)) + layer(panel.abline(v = md, col = blue, lty = 2)) On Wed, Nov 19, 2014 at 2:56 AM, Nuno Sá [email protected] wrote: Hello! My aim is to add a Median or a Mean line to an histogram plot in R. The problem is the following: The hist function from the raster package uses a maximum of 100 000 values for generating the histogram but allows me to edit the plot, so I can easily add a line within the plot using abline. The histogram function in rasterVIS uses all my dataset but does not allow editing of the plot area, so I cannot use the abline function to add this line (or I do not know how to do it). If you have an alternative or a solution to work around this, I am all hears Thank you for any help in advance! Ciao! Nuno César de Sá +351 91 961 90 37 [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Pascal Oettli Project Scientist JAMSTEC Yokohama, Japan R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo The hist function from the raster package uses a maximum of 100 000 values for generating the histogram That is the default, but you do not have to use that. Did you read the help file?? ?raster::hist Shows that you can set the maxpixels to any number you want. E.g. maxpixels=Inf if you want them all. With very large raster you may run out of memory. If that happens, then you can create the data for a histogram yourself by using the function 'freq', perhaps after 'reclassify' (and then probably 'barplot'). Robert My aim is to add a Median or a Mean line to an histogram plot in R. The problem is the following: The hist function from the raster package uses a maximum of 100 000 values for generating the histogram but allows me to edit the plot, so I can easily add a line within the plot using abline. The histogram function in rasterVIS uses all my dataset but does not allow editing of the plot area, so I cannot use the abline function to add this line (or I do not know how to do it). If you have an alternative or a solution to work around this, I am all hears Thank you for any help in advance! Ciao! Nuno César de Sá +351 91 961 90 37 [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Mr. Pascal, thank you, your suggestion works quite well. First, I was having error calculating the median using the cellStats because I was using brackets (median as opposed to median) and the error message mislead me to think it was not implemented (which obviously I obviously found very strange but moved quickly onwards after Thomas Adams suggestion). (This was the error) Error in .local(x, stat, ...) : invalid 'stat'. Should be sum, min, max, sd, mean, or 'countNA' The + operator worked 100% as you've said. So great feedback, thank you soo much.I am now free to make some more beautiful looking plots hehe Mr. Hijmanns, I did read the package help but failed to find a mention that it was not limited to a maximum of 100, so my mistake there. It works perfectly as you've said as well, so thank you so much for the feedback. As usual, computers are always right and the user was the error. Good thing we have these forums and such a quick access to much better users. Since I have your attention Mr. Hijmanns I ask a quick extra question - in the cellStats you mention that it will fail gracefully for very large rasters (not my case here since it gracefully succeeded), do you know what is the limit of cells for that failure? Once again, thank you all for the help! Best of luck on everything for all! On 19 November 2014 15:14, Robert J. Hijmans [email protected] wrote: Nuno, The hist function from the raster package uses a maximum of 100 000 values for generating the histogram That is the default, but you do not have to use that. Did you read the help file?? ?raster::hist Shows that you can set the maxpixels to any number you want. E.g. maxpixels=Inf if you want them all. With very large raster you may run out of memory. If that happens, then you can create the data for a histogram yourself by using the function 'freq', perhaps after 'reclassify' (and then probably 'barplot'). Robert My aim is to add a Median or a Mean line to an histogram plot in R. The problem is the following: The hist function from the raster package uses a maximum of 100 000 values for generating the histogram but allows me to edit the plot, so I can easily add a line within the plot using abline. The histogram function in rasterVIS uses all my dataset but does not allow editing of the plot area, so I cannot use the abline function to this line (or I do not know how to do it). If you have an alternative or a solution to work around this, I am all hears Thank you for any help in advance! Ciao! Nuno César de Sá +351 91 961 90 37 [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo cellStats will fail gracefully for very large rasters do you know what is the limit of cells for that failure? That depends on your computer hardware (RAM), software (32 or 64 bits OS) and perhaps memory management (will RAM expand to hard disk if needed?). On Wed, Nov 19, 2014 at 11:39 AM, Nuno Sá [email protected] wrote: Hello! Mr. Pascal, thank you, your suggestion works quite well. First, I was having error calculating the median using the cellStats because I was using brackets (median as opposed to median) and the error message mislead me to think it was not implemented (which obviously I obviously found very strange but moved quickly onwards after Thomas Adams suggestion). (This was the error) Error in .local(x, stat, ...) : invalid 'stat'. Should be sum, min, max, sd, mean, or 'countNA' The + operator worked 100% as you've said. So great feedback, thank you soo much.I am now free to make some more beautiful looking plots hehe Mr. Hijmanns, I did read the package help but failed to find a mention that it was not limited to a maximum of 100, so my mistake there. It works perfectly as you've said as well, so thank you so much for the feedback. As usual, computers are always right and the user was the error. Good thing we have these forums and such a quick access to much better users. Since I have your attention Mr. Hijmanns I ask a quick extra question - in the cellStats you mention that it will fail gracefully for very large rasters (not my case here since it gracefully succeeded), do you know what is the limit of cells for that failure? Once again, thank you all for the help! Best of luck on everything for all! On 19 November 2014 15:14, Robert J. Hijmans [email protected] wrote: Nuno, The hist function from the raster package uses a maximum of 100 000 values for generating the histogram That is the default, but you do not have to use that. Did you read the help file?? ?raster::hist Shows that you can set the maxpixels to any number you want. E.g. maxpixels=Inf if you want them all. With very large raster you may run out of memory. If that happens, then you can create the data for a histogram yourself by using the function 'freq', perhaps after 'reclassify' (and then probably 'barplot'). Robert My aim is to add a Median or a Mean line to an histogram plot in R. The problem is the following: The hist function from the raster package uses a maximum of 100 000 values for generating the histogram but allows me to edit the plot, so I can easily add a line within the plot using abline. The histogram function in rasterVIS uses all my dataset but does not allow editing of the plot area, so I cannot use the abline function to this line (or I do not know how to do it). If you have an alternative or a solution to work around this, I am all hears Thank you for any help in advance! Ciao! Nuno César de Sá +351 91 961 90 37 [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo My aim is to add a Median or a Mean line to an histogram plot in R. The problem is the following: The hist function from the raster package uses a maximum of 100 000 values for generating the histogram but allows me to edit the plot, so I can easily add a line within the plot using abline. The histogram function in rasterVIS uses all my dataset but does not allow editing of the plot area, so I cannot use the abline function to add this line (or I do not know how to do it). If you have an alternative or a solution to work around this, I am all hears Thank you for any help in advance! Ciao! Nuno César de Sá +351 91 961 90 37 [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Yes thanks, I know it works for the normal hist function which uses 3% of all data (mind I am talking about rasters) which I am not sure that is a good enough sample, I leave that question in the air The problem is it does not work in rasterVIS (which does use all the raster). The error it gives me is the following: abline(v=m.rst,col=blue,lwd=2) Error in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...) : plot.new has not been called yet Regards!! Thank you once again! On 18 November 2014 18:04, Thomas Adams [email protected] wrote: Nuno, This works for me: require(stats) set.seed(14) x - rchisq(100, df = 4) hist(x) meanx-mean(x) meanx [1] 3.734162 abline(v=meanx, col=red) Cheers! On Tue, Nov 18, 2014 at 10:56 AM, Nuno Sá [email protected] wrote: Hello! My aim is to add a Median or a Mean line to an histogram plot in R. The problem is the following: The hist function from the raster package uses a maximum of 100 000 values for generating the histogram but allows me to edit the plot, so I can easily add a line within the plot using abline. The histogram function in rasterVIS uses all my dataset but does not allow editing of the plot area, so I cannot use the abline function to add this line (or I do not know how to do it). If you have an alternative or a solution to work around this, I am all hears Thank you for any help in advance! Ciao! Nuno César de Sá +351 91 961 90 37 [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo I thought you might say that. I have used R with GRASS GIS, which interfaces using the spgrass6 contributed package. I have analyzed some pretty large GIS raster files, looking at radar precipitation data, and the approach I suggested works well. The data imported into R using sp (the basis for spgrass6) is an sp raster object. So, from the spgrass6 manual: After loading R from the GRASS GIS terminal prompt, for example, spear - readRAST6(c(geology, elevation.dem), cat=c(TRUE, FALSE)) Try looking at this: summary(spear) you can do this: boxplot(spear$elevation.dem ~ spear$geology) or this, for instance : hist(spear$elevation.dem) meanx-mean(spear$elevation.dem) abline(v=meanx, col=red) But, maybe I don't understand the problem?? Regards, On Tue, Nov 18, 2014 at 11:21 AM, Nuno Sá [email protected] wrote: Hey Thomas! Yes thanks, I know it works for the normal hist function which uses 3% of all data (mind I am talking about rasters) which I am not sure that is a good enough sample, I leave that question in the air The problem is it does not work in rasterVIS (which does use all the raster). The error it gives me is the following: abline(v=m.rst,col=blue,lwd=2) Error in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...) : plot.new has not been called yet Regards!! Thank you once again! On 18 November 2014 18:04, Thomas Adams [email protected] wrote: Nuno, This works for me: require(stats) set.seed(14) x - rchisq(100, df = 4) hist(x) meanx-mean(x) meanx [1] 3.734162 abline(v=meanx, col=red) Cheers! On Tue, Nov 18, 2014 at 10:56 AM, Nuno Sá [email protected] wrote: Hello! My aim is to add a Median or a Mean line to an histogram plot in R. The problem is the following: The hist function from the raster package uses a maximum of 100 000 values for generating the histogram but allows me to edit the plot, so I can easily add a line within the plot using abline. The histogram function in rasterVIS uses all my dataset but does not allow editing of the plot area, so I cannot use the abline function to this line (or I do not know how to do it). If you have an alternative or a solution to work around this, I am all hears Thank you for any help in advance! Ciao! Nuno César de Sá +351 91 961 90 37 [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Does rasterVIS use base or grid graphics? The place I usually see the error that plot.new has not been called yet is when trying to use a base graphics function to augment a plot created using grid graphics. I suspect that the plot created used grid graphics and that is why abline is not working for you. There is a 'grid.abline' function in the grid package, though it does not look like it would work easily for a vertical line, so you may need to dive into the details of using grid in more detail, probably using the grid.lines function. There is also the gridBase package which helps grid and base graphics play nicely together, though it could be a bit of overkill for just adding a vertical line. On Tue, Nov 18, 2014 at 11:21 AM, Nuno S� [email protected] wrote: Hey Thomas! Yes thanks, I know it works for the normal hist function which uses 3% of all data (mind I am talking about rasters) which I am not sure that is a good enough sample, I leave that question in the air The problem is it does not work in rasterVIS (which does use all the raster). The error it gives me is the following: abline(v=m.rst,col=blue,lwd=2) Error in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...) : plot.new has not been called yet Regards!! Thank you once again! On 18 November 2014 18:04, Thomas Adams [email protected] wrote: Nuno, This works for me: require(stats) set.seed(14) x - rchisq(100, df = 4) hist(x) meanx-mean(x) meanx [1] 3.734162 abline(v=meanx, col=red) Cheers! On Tue, Nov 18, 2014 at 10:56 AM, Nuno S� [email protected] wrote: Hello! My aim is to add a Median or a Mean line to an histogram plot in R. The problem is the following: The hist function from the raster package uses a maximum of 100 000 values for generating the histogram but allows me to edit the plot, so I can easily add a line within the plot using abline. The histogram function in rasterVIS uses all my dataset but does not allow editing of the plot area, so I cannot use the abline function to this line (or I do not know how to do it). If you have an alternative or a solution to work around this, I am all hears Thank you for any help in advance! Ciao! Nuno C�sar de S� +351 91 961 90 37 [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Possibly having a reproducible example would help. I do all my spatial analysis work starting from GRASS GIS and calling R from the GRASS prompt — this works very well. On Tue, Nov 18, 2014 at 11:35 AM, Greg Snow [email protected] wrote: Nuno, Does rasterVIS use base or grid graphics? The place I usually see the error that plot.new has not been called yet is when trying to use a base graphics function to augment a plot created using grid graphics. I suspect that the plot created used grid graphics and that is why abline is not working for you. There is a 'grid.abline' function in the grid package, though it does not look like it would work easily for a vertical line, so you may need to dive into the details of using grid in more detail, probably using the grid.lines function. There is also the gridBase package which helps grid and base graphics play nicely together, though it could be a bit of overkill for just adding a vertical line. On Tue, Nov 18, 2014 at 11:21 AM, Nuno Sá [email protected] wrote: Hey Thomas! Yes thanks, I know it works for the normal hist function which uses 3% of all data (mind I am talking about rasters) which I am not sure that is a good enough sample, I leave that question in the air The problem is it does not work in rasterVIS (which does use all the raster). The error it gives me is the following: abline(v=m.rst,col=blue,lwd=2) Error in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...) : plot.new has not been called yet Regards!! Thank you once again! On 18 November 2014 18:04, Thomas Adams [email protected] wrote: Nuno, This works for me: require(stats) set.seed(14) x - rchisq(100, df = 4) hist(x) meanx-mean(x) meanx [1] 3.734162 abline(v=meanx, col=red) Cheers! On Tue, Nov 18, 2014 at 10:56 AM, Nuno Sá [email protected] wrote: Hello! My aim is to add a Median or a Mean line to an histogram plot in R. The problem is the following: The hist function from the raster package uses a maximum of 100 000 values for generating the histogram but allows me to edit the plot, so can easily add a line within the plot using abline. The histogram function in rasterVIS uses all my dataset but does not allow editing of the plot area, so I cannot use the abline function to this line (or I do not know how to do it). If you have an alternative or a solution to work around this, I am all hears Thank you for any help in advance! Ciao! Nuno César de Sá +351 91 961 90 37 [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Yes, it sorted the problem but maybe not as we were expecting, you've liberated me from the raster package 3% data limitation and now I can work with it easily as a normal histogram. Basically the trick that I got from the readRAST6 function you suggested is that it converts it into sp objects (SpatialGridDataFrame in this case i guess). From here it is easy: all I need is to convert the rasters to an sp object and from there I am back working with normal hist functions. So, using the conversion: as(inputraster,SpatialGridDataFrame) sorted it out. The rest is history. I hope they update this method into the raster package because there is no need to be limited to a 100 000 sample which is usually kind of low for RS data. Found a similar limitation on usdm package recently which chooses to needlessly convert a whole raster into a data.frame before calculating VIF between variables. For huge rasters, it would kill off the computer memory. Because of that, I guess the authors also made a limit on the sample size (which was not what was actually killing the computer memory) Dr Greg Yes, I think the rasterVIS package uses grid. maybe it can be inserted when we call the histogram function, but making it do a vertical line might be kind of complicated, I think just from looking to how it is called. Thank you guys for the suggestions!! I jump between them also, but making graphs is usually easy enough just using R. anyway, it is so great to have all these awesome languages and softwares linked together. Best of luck for all! On 18 November 2014 18:38, Thomas Adams [email protected] wrote: Nuno, Possibly having a reproducible example would help. I do all my spatial analysis work starting from GRASS GIS and calling R from the GRASS prompt — this works very well. On Tue, Nov 18, 2014 at 11:35 AM, Greg Snow [email protected] wrote: Nuno, Does rasterVIS use base or grid graphics? The place I usually see the error that plot.new has not been called yet is when trying to use a base graphics function to augment a plot created using grid graphics. I suspect that the plot created used grid graphics and that is why abline is not working for you. There is a 'grid.abline' function in the grid package, though it does not look like it would work easily for a vertical line, so you may need to dive into the details of using grid in more detail, probably using the grid.lines function. There is also the gridBase package which helps grid and base graphics play nicely together, though it could be a bit of overkill for just adding a vertical line. On Tue, Nov 18, 2014 at 11:21 AM, Nuno Sá [email protected] wrote: Hey Thomas! Yes thanks, I know it works for the normal hist function which uses 3% of all data (mind I am talking about rasters) which I am not sure that is a good enough sample, I leave that question in the air The problem is it does not work in rasterVIS (which does use all the raster). The error it gives me is the following: abline(v=m.rst,col=blue,lwd=2) Error in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...) : plot.new has not been called yet Regards!! Thank you once again! On 18 November 2014 18:04, Thomas Adams [email protected] wrote: Nuno, This works for me: require(stats) set.seed(14) x - rchisq(100, df = 4) hist(x) meanx-mean(x) meanx [1] 3.734162 abline(v=meanx, col=red) Cheers! On Tue, Nov 18, 2014 at 10:56 AM, Nuno Sá [email protected] wrote: Hello! My aim is to add a Median or a Mean line to an histogram plot in The problem is the following: The hist function from the raster package uses a maximum of 100 000 values for generating the histogram but allows me to edit the plot, can easily add a line within the plot using abline. The histogram function in rasterVIS uses all my dataset but does not allow editing of the plot area, so I cannot use the abline function to add this line (or I do not know how to do it). If you have an alternative or a solution to work around this, I am all hears Thank you for any help in advance! Ciao! Nuno César de Sá +351 91 961 90 37 [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo I miss you guys from GEOSTAT Bergen! Here's my question: I want to make a map directly from R for my manuscript, using a raster layer overlaid with spatial polygons as outlines. Typically, you use the sp.layout= argument to spplot. This takes a list of extra information to be plotted. However, spatial polygons are thought of as filled, so you should coerce then to SpatialLines first: spl - list(sp.lines, as(my_polys, SpatialLines), lwd=2, col=white) spplot(..., sp.layout=spl) or a list of lists if there are other kinds of information to add. Hope this helps, Roger Please see the attachment for what I mean: http://1drv.ms/1sYYqPg I did the upper plot in the attachement using these codes: NRI = stack(NRI.A,NRI.G) # NRI.A and NRI.G are raster layer names names(NRI) = c(Angiosperms NRI, Gymnosperms NRI) spplot(NRI) I need to put two rasters together in spplot because I want to use the same legend scale. I also want to add the outlines, so that it can look like the Photoshopped plot in the attachment. Normally, if not using spplot, I can do this: plot(NRI.G) # NRI.G is a raster layer plot(Outlines) # Outlines is a SpatialPolygons object But how can I do it for both panels in spplot? I think the panels and titles in spplot are very neat, so it is preferred if there is a straight forward way to tweak the settings. Or should I use other functions than spplot? (e.g. try coding the color scheme of the 2 plots and place them side by side?) Thank you for the help. I really appreciate those who made R able to produce publication quality plots. Cheers, Ziyu Ma PhD Student Ecoinformatics Biodiversity, Department of Bioscience, Aarhus University Ny Munkegade 114, DK-8000 Aarhus C, Denmark [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Roger Bivand Department of Economics, Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 91 00 e-mail: [email protected] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo You should be able to plot the polygon on top of your raster with a panel function. Here is an example with the meuse data: r - raster(system.file(external/test.grd, package=raster)) s - stack(r, r*2) names(s) - c('meuse', 'meuse x 2') spplot(s) data(meuse.riv) meuse.sr = SpatialPolygons(list(Polygons(list(Polygon(meuse.riv)),meuse.riv))) spplot(s, panel = function(x, y, ...) { panel.gridplot(x, y, ...) sp.polygons(meuse.sr) Is this more or less how you want it? Cheers, On 8/22/2014 9:15 AM, Ziyu Ma wrote: Dear all, I miss you guys from GEOSTAT Bergen! Here's my question: I want to make a map directly from R for my manuscript, using a raster layer overlaid with spatial polygons as outlines. Please see the attachment for what I mean: http://1drv.ms/1sYYqPg I did the upper plot in the attachement using these codes: NRI = stack(NRI.A,NRI.G) # NRI.A and NRI.G are raster layer names names(NRI) = c(Angiosperms NRI, Gymnosperms NRI) spplot(NRI) I need to put two rasters together in spplot because I want to use the same legend scale. I also want to add the outlines, so that it can look like the Photoshopped plot in the attachment. Normally, if not using spplot, I can do this: plot(NRI.G) # NRI.G is a raster layer plot(Outlines) # Outlines is a SpatialPolygons object But how can I do it for both panels in spplot? I think the panels and titles in spplot are very neat, so it is preferred if there is a straight forward way to tweak the settings. Or should I use other functions than spplot? (e.g. try coding the color scheme of the 2 plots and place them side by side?) Thank you for the help. I really appreciate those who made R able to produce publication quality plots. Cheers, Ziyu Ma PhD Student Ecoinformatics Biodiversity, Department of Bioscience, Aarhus University Ny Munkegade 114, DK-8000 Aarhus C, Denmark [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo You should be able to plot the polygon on top of your raster with a panel function. Here is an example with the meuse data: r - raster(system.file(external/test.grd, package=raster)) s - stack(r, r*2) names(s) - c('meuse', 'meuse x 2') spplot(s) data(meuse.riv) meuse.sr = SpatialPolygons(list(Polygons(list(Polygon(meuse.riv)), meuse.riv))) spplot(s, panel = function(x, y, ...) { panel.gridplot(x, y, ...) sp.polygons(meuse.sr) Is this more or less how you want it? Cheers, On 8/22/2014 9:15 AM, Ziyu Ma wrote: Dear all, I miss you guys from GEOSTAT Bergen! Here's my question: I want to make a map directly from R for my manuscript, using a raster layer overlaid with spatial polygons as outlines. Please see the attachment for what I mean: http://1drv.ms/1sYYqPg I did the upper plot in the attachement using these codes: NRI = stack(NRI.A,NRI.G) # NRI.A and NRI.G are raster layer names names(NRI) = c(Angiosperms NRI, Gymnosperms NRI) spplot(NRI) I need to put two rasters together in spplot because I want to use the legend scale. I also want to add the outlines, so that it can look like the Photoshopped plot in the attachment. Normally, if not using spplot, I can do this: plot(NRI.G) # NRI.G is a raster layer plot(Outlines) # Outlines is a SpatialPolygons object But how can I do it for both panels in spplot? I think the panels and titles in spplot are very neat, so it is preferred if there is a straight forward way to tweak the settings. Or should I use other functions than spplot? (e.g. try coding the color scheme of the 2 plots and place them side by side?) Thank you for the help. I really appreciate those who made R able to produce publication quality plots. Cheers, Ziyu Ma PhD Student Ecoinformatics Biodiversity, Department of Bioscience, Aarhus University Ny Munkegade 114, DK-8000 Aarhus C, Denmark [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo The result you got is the same as for a variogram model from vgm with a nugget effect. This is the ratio of anisotropy, which is one both for the nugget and for the correlated part of the model. See: str(vgm(10, Exp, 300)) str(vgm(10, Exp, 300, 3)) I guess you see something like the second one as autofitVariogram is not able to fit the anisotropy parameters. For that you could have a look at estimateParameters in the intamap package, which can fit a variogram model of the vgm-type with anisotropy parameters. Cheers, On 7/29/2014 6:56 PM, Hodgess, Erin wrote: Hello again! I was using autoKrige (finally successfully!) and was looking at the model results. I got anis1 of 1,1 and anis2 of 1,1. How does this relate to the anis in the vgm function, please? Thanks, [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Sorry for simple question. Is it correct to say that a CAR model is a spatial lag model and a SAR model is a error spatial model.? No, it is not. The spatial lag and spatial error models are both SAR - simultaneous autoregressive, while CAR is conditional autoregressive and is not directly comparable. See books by Ripley (1981) or Haining (2003). Roger I revised some literature where these concepts are not changed. [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Roger Bivand Department of Economics, Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: [email protected] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo unique(cbind(bigz$lon,bigz$lat)) Nothing like posting to the list to help you solve a problem. thanks, From: [email protected] [[email protected]] on behalf of Hodgess, Erin [[email protected]] Sent: Saturday, February 15, 2014 8:30 PM To: [email protected] Subject: [R-sig-Geo] question about unique in lon/lat Hello! I have space-time data frame. Here are the first few rows: name indexlon lat 1 Abilene,TX 95.285 -99.73314 32.44874 2 Abilene,TX 98.603 -99.73314 32.44874 3 Abilene,TX 100.197 -99.73314 32.44874 4 Abilene,TX 102.016 -99.73314 32.44874 5 Abilene,TX 100.000 -99.73314 32.44874 6 Abilene,TX 102.822 -99.73314 32.44874 7 Abilene,TX 107.988 -99.73314 32.44874 8 Abilene,TX 110.546 -99.73314 32.44874 9 Abilene,TX 105.897 -99.73314 32.44874 10 Abilene,TX 107.642 -99.73314 32.44874 11 Abilene,TX 108.523 -99.73314 32.44874 12 Amarillo,TX 92.486 -101.83130 35.22200 13 Amarillo,TX 95.431 -101.83130 35.22200 14 Amarillo,TX 97.485 -101.83130 35.22200 15 Amarillo,TX 99.714 -101.83130 35.22200 16 Amarillo,TX 100.000 -101.83130 35.22200 17 Amarillo,TX 103.649 -101.83130 35.22200 18 Amarillo,TX 105.998 -101.83130 35.22200 19 Amarillo,TX 109.965 -101.83130 35.22200 20 Amarillo,TX 109.325 -101.83130 35.22200 21 Amarillo,TX 112.550 -101.83130 35.22200 22 Amarillo,TX 112.541 -101.83130 35.22200 23 Austin ,TX 83.228 -97.74306 30.26715 24 Austin ,TX 83.232 -97.74306 30.26715 25 Austin ,TX 86.072 -97.74306 30.26715 26 Austin ,TX 93.592 -97.74306 30.26715 My question is: is there a way to get the unique values from lon/lat for each city, please? Thanks, [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo I tried the following example: erin1 - summary(COL.mix.eig, correlation=TRUE, Nagelkerke=TRUE) names(erin1) [1] typerho coefficientsrest.se [5] LL s2 SSE parameters [9] logLik_lm.model AIC_lm.modelmethod call [13] residuals opt tarXtary [17] y X fitted.values se.fit [21] similar ase rho.se LMtest [25] resvar zero.policy aliased listw_style [29] intervalfdHess optimHess insert [33] trs LLNullLlm timings f_calls [37] hf_callsintern_classic coeftitle Coef [41] NK Wald1 correlation correltext [45] LR1 erin1$LMtest [1,] 0.2891926 and it does indeed have the LMtest result. Or were you looking for the formula, please? Thanks, From: [email protected] [[email protected]] on behalf of Dongwoo Kang [[email protected]] Sent: Tuesday, July 09, 2013 3:47 PM To: [email protected] Subject: [R-sig-Geo] Question about LM test for residual autocorrelation in R Dear list, Hello, I am Dongwoo Kang. I am studying Spatial Econometric modeling. I've faced one question while using *spdep* package in R. I want to ask your help for my qeustion. While estimating my empirical models, I want to test whether residuals of my spatial regression models (SEM, SAR, SARAR, SDM estimated by maximum likelihood) still have spatial autocorrelation pattern. I think I have two options, 1) Moran's I test using *moran.mc* function in R, 2) Lagrange multiplier diagnostics with LMerr option using *lm.LMtests* function in R. But I also find that for SAR, SDM, *summary.sarlm* function returns LM test for residual autocorrelation by default. However, this LM test is not given for SER and SARAR. At first, I thought that Lagrange multiplier diagnostics and LM test for residual autocorrelation in *summary.sarlm* function are same tests. But in my empirical results, they give me different statistics (please see below example). - example summary.sarlm(sar2, Nagelkerke=T) Log likelihood: -3533.378 for lag model ML residual variance (sigma squared): 224.88, (sigma: 14.996) Nagelkerke pseudo-R-squared: 0.76166 Number of observations: 853 Number of parameters estimated: 29 AIC: 7124.8, (AIC for lm: 7393.3) LM test for residual autocorrelation test value: 6.8391, p-value: 0.0089184 lm.LMtests(sar2$residuals, listw=w100.listw, test=c(LMerr)) Lagrange multiplier diagnostics for spatial dependence data: residuals: sar2$residuals weights: w100.listw LMErr = 3.7108, df = 1, p-value = 0.05406 I try to find formulation of LM test for residual autocorrelation given by *summary.sarlm* function but I couldn't. Would you tell me where I can get some documents or explanations about LM test for residual autocorrelation given by *summary.sarlm*? I also want to know why LM test for residual autocorrelation is not provided in SER and SARAR models. Thank you very much for your time. Best regards, Dongwoo Kang [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo From: [email protected] [[email protected]] on behalf of Dongwoo Kang [[email protected]] Sent: Tuesday, July 09, 2013 3:47 PM To: [email protected] Subject: [R-sig-Geo] Question about LM test for residual autocorrelation in R Dear list, Hello, I am Dongwoo Kang. I am studying Spatial Econometric modeling. I've faced one question while using *spdep* package in R. I want to ask your help for my qeustion. While estimating my empirical models, I want to test whether residuals of my spatial regression models (SEM, SAR, SARAR, SDM estimated by maximum likelihood) still have spatial autocorrelation pattern. I think I have two options, 1) Moran's I test using *moran.mc* function in R, 2) Lagrange multiplier diagnostics with LMerr option using *lm.LMtests* function in R. But I also find that for SAR, SDM, *summary.sarlm* function returns LM test for residual autocorrelation by default. However, this LM test is not given for SER and SARAR. At first, I thought that Lagrange multiplier diagnostics and LM test for residual autocorrelation in *summary.sarlm* function are same tests. But in my empirical results, they give me different statistics (please see below example). - example summary.sarlm(sar2, Nagelkerke=T) Log likelihood: -3533.378 for lag model ML residual variance (sigma squared): 224.88, (sigma: 14.996) Nagelkerke pseudo-R-squared: 0.76166 Number of observations: 853 Number of parameters estimated: 29 AIC: 7124.8, (AIC for lm: 7393.3) LM test for residual autocorrelation test value: 6.8391, p-value: 0.0089184 lm.LMtests(sar2$residuals, listw=w100.listw, test=c(LMerr)) Lagrange multiplier diagnostics for spatial dependence data: residuals: sar2$residuals weights: w100.listw LMErr = 3.7108, df = 1, p-value = 0.05406 I try to find formulation of LM test for residual autocorrelation given by *summary.sarlm* function but I couldn't. Would you tell me where I can get some documents or explanations about LM test for residual autocorrelation given by *summary.sarlm*? I also want to know why LM test for residual autocorrelation is not provided in SER and SARAR models. Thank you very much for your time. Best regards, Dongwoo Kang [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Neither R nor geo, but the free program ImageJ is pretty good for this kind of image analysis. http://rsb.info.nih.gov/ij/ Cheers, David On 1 Apr 2013, at 15:13, A.P.B. Carneiro [email protected]:[email protected] wrote: Dear all, I am trying to calculate an area from a picture. I think I can use Spatial Analyst but I am not sure how to do it. I have a picture from the wing of a bird, and from there, I want to calculate the area of a white patch. The picture is with a scale bar. Is there a way to do that? There is not spatial information on the picture. Thanks in advance, On Sep 29 2012, [email protected]:[email protected] wrote: Send R-sig-Geo mailing list submissions to [email protected]:[email protected] To subscribe or unsubscribe via the World Wide Web, visit https://stat.ethz.ch/mailman/listinfo/r-sig-geo or, via email, send a message with subject or body 'help' to [email protected] You can reach the person managing the list at [email protected] When replying, please edit your Subject line so it is more specific than Re: Contents of R-sig-Geo digest... Today's Topics: 1. Re: plotting holes of esri shape files (Edzer Pebesma) 2. Re: basic question about plotKML (Struve,Juliane) 3. Re: basic question about plotKML (Edzer Pebesma) 4. Re: basic question about plotKML (Struve,Juliane) 5. Re: Create a raster brick from a netcdf using all levels (Robert J. Hijmans) 6. Re: NaN / NA values in SpatialPixelDataFrame conversion (Robert J. Hijmans) 7. Re: Error when trying to extract values for points in a SpatialPointsDataFrame (Robert J. Hijmans) Message: 1 Date: Fri, 28 Sep 2012 13:18:18 +0200 From: Edzer Pebesma [email protected] To: [email protected] [email protected] Subject: Re: [R-sig-Geo] plotting holes of esri shape files Message-ID: [email protected] Content-Type: text/plain; charset=ISO-8859-1 Stefan, with the help of http://www.stat.auckland.ac.nz/~paul/R/Paths/MurrellPaths.pdf I corrected this in spplot. To run it, you'd need to compile from source on svn, or source this in your R session: https://r-forge.r-project.org/scm/viewvc.php/*checkout*/pkg/sp/R/spplot.R?root=rspatial or wait for the next release on CRAN (due soon, I hope). spplot now also supports holes in polygons when the polygons are the main object to plot, contain holes, and something in the background should shine through. Building on the example in https://stat.ethz.ch/pipermail/r-sig-geo/2012-September/016116.html p1 = spplot(mask, col.regions='lightgrey', sp.layout = list(sp.grid, as(dem_ascii, SpatialGridDataFrame), col = bpy.colors(), first = TRUE), main = spplot of mask) p2 = spplot(dem_ascii, col.regions=bpy.colors(), main = spplot of dem, sp.layout = list(sp.polygons, mask, fill = 'lightgray', first = FALSE)) png('spplot.png', 600, 600) print(p1, split = c(1,1,1,2), more = TRUE) print(p2, split = c(1,2,1,2)) dev.off() gives http://ifgi.uni-muenster.de/~epebe_01/spplot.png as output. Mind the legends. With best regards, On 09/20/2012 10:32 AM, sluedtke wrote: Dear friends, am trying to overlay a polygon data set as a SpatialPolygonDateFrame to an existing plot. The data set has a hole in the middle, and I do not want to have any filling in there. If a check the shape in a GIS, everything looks fine, but using the following commands to do the plot, the hole is filled. That is the code I am using, the data is attached. Thanks in advance! Stefan library(sp) library(rgdal) library(raster) library(latticeExtra) #read the DEM from the ASCII format dem_ascii=raster(readGDAL(./zerafshan_mod_dem.asc)) # read additional data for printing mask=readOGR(./, mask_zerafshan) #plot this stuff plot_dem=spplot(dem_ascii, col.regions=terrain.colors(400), colorkey=list(width=1, space=right))+ layer(sp.polygons(mask,fill=lightgray)) print(plot_dem) mask_zerafshan.dbf (110 bytes) http://r-sig-geo.2731867.n2.nabble.com/attachment/7580992/0/mask_zerafshan.dbf mask_zerafshan.prj (143 bytes) http://r-sig-geo.2731867.n2.nabble.com/attachment/7580992/1/mask_zerafshan.prj mask_zerafshan.shp (173K) http://r-sig-geo.2731867.n2.nabble.com/attachment/7580992/2/mask_zerafshan.shp mask_zerafshan.shx (150 bytes) http://r-sig-geo.2731867.n2.nabble.com/attachment/7580992/3/mask_zerafshan.shxRe: [R-sig-Geo] Question about geoTIFF
2013-03-19 Thread Pascal Oettli If you use the raster package to read the geoTIFF, you can use getValues (whole data) or extract (spatial selection). Pascal On 19/03/13 17:02, Daniel Flø wrote: I have downloaded forestry inventory file in geoTIFF format successfully, and the goal is to extract the data of the cells in the geoTIFF file for use in R calculations. So far, I use the raster package to handle the TIFF file (display map etc), but I cannot figure out how to extract the data and store it to a .csv file. Any one that can help? Thank you. Daniel Flø PhD Fellow, Forest Health The Norwegian Forest and Landscape Institute Pb 115, NO-1431 Ås Office: (+47) 64 94 90 28 Mobile: (+47) 91 14 52 74 Skype:floe.daniel www.skogoglandskap.nohttp://www.skogoglandskap.no/ [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo I don't see filename defined, or any packages required. Is it the png files that are a problem or the 'kmlsetup' output? (whatever that is) See the animation package for extremely simplifed image frame creation btw On Wednesday, March 6, 2013, Hodgess, Erin wrote: Hello again! I'm running into a problem with the image function within a loop. When I put it to the screen, it works fine, but when I put it to a PNG file, the picture doesn't change. Here is the code: for(i in 1:11) { j1 - (i-1)*n2 + 1 j2 - i*n2 predx - data.frame(pred=kr3$var1.pred[j1:j2],lon=kr3@sp @coords[,1],lat=kr3@sp@coords[,2]) coordinates(predx) - ~lon+lat proj4string(predx) - proj4string(kr3) predx.wgs84 - spTransform(predx,CRS(+proj=longlat +datum=WGS84)) llSPix - as(llSGDF, SpatialPixelsDataFrame) (st - system.time(llSPix$pred - idw(pred ~ 1, loc=predx.wgs84,newdata = llSPix, nmax = 16)$var1.pred)) cat(stuff,i,\n) f1 - paste(pred,i,a.png,sep=) #Works fine image(llSPix,pred,col=bpy.colors(),main=f1) assign(file1,f1) cat(filen,file1,\n) png(filename = file1, width = GRD.wgs84$width, height = GRD.wgs84$height,bg=transparent) par(mar = c(0, 0, 0, 0), xaxs = i, yaxs = i) image(llSPix, pred, col = bpy.colors()) dev.off() kmlsetup(obj = GRD.wgs84, kmlfile =filename, imagefile = file1, name =as.character(2000+i)) Has anyone run into this before, please? Thanks, [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] javascript:; https://stat.ethz.ch/mailman/listinfo/r-sig-geo Michael Sumner Hobart, Australia e-mail: [email protected] [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo When using png in a loop, especially one where the loop index is an integer, I usually make use of the %02d in the file name, e.g., png(pred_%02d_a.png,...) Clint BowmanINTERNET: [email protected] Air Quality Modeler INTERNET: [email protected] Department of Ecology VOICE: (360) 407-6815 PO Box 47600FAX:(360) 407-7534 Olympia, WA 98504-7600 USPS: PO Box 47600, Olympia, WA 98504-7600 Parcels:300 Desmond Drive, Lacey, WA 98503-1274 On Wed, 6 Mar 2013, Michael Sumner wrote: I don't see filename defined, or any packages required. Is it the png files that are a problem or the 'kmlsetup' output? (whatever that is) See the animation package for extremely simplifed image frame creation btw On Wednesday, March 6, 2013, Hodgess, Erin wrote: Hello again! I'm running into a problem with the image function within a loop. When I put it to the screen, it works fine, but when I put it to a PNG file, the picture doesn't change. Here is the code: for(i in 1:11) { j1 - (i-1)*n2 + 1 j2 - i*n2 predx - data.frame(pred=kr3$var1.pred[j1:j2],lon=kr3@sp @coords[,1],lat=kr3@sp@coords[,2]) coordinates(predx) - ~lon+lat proj4string(predx) - proj4string(kr3) predx.wgs84 - spTransform(predx,CRS(+proj=longlat +datum=WGS84)) llSPix - as(llSGDF, SpatialPixelsDataFrame) (st - system.time(llSPix$pred - idw(pred ~ 1, loc=predx.wgs84,newdata = llSPix, nmax = 16)$var1.pred)) cat(stuff,i,\n) f1 - paste(pred,i,a.png,sep=) #Works fine image(llSPix,pred,col=bpy.colors(),main=f1) assign(file1,f1) cat(filen,file1,\n) png(filename = file1, width = GRD.wgs84$width, height = GRD.wgs84$height,bg=transparent) par(mar = c(0, 0, 0, 0), xaxs = i, yaxs = i) image(llSPix, pred, col = bpy.colors()) dev.off() kmlsetup(obj = GRD.wgs84, kmlfile =filename, imagefile = file1, name =as.character(2000+i)) Has anyone run into this before, please? Thanks, [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] javascript:; https://stat.ethz.ch/mailman/listinfo/r-sig-geo From: Michael Sumner [[email protected]] Sent: Tuesday, March 05, 2013 2:41 PM To: Hodgess, Erin Cc: [email protected] Subject: Re: [R-sig-Geo] question about image function I don't see filename defined, or any packages required. Is it the png files that are a problem or the 'kmlsetup' output? (whatever that is) See the animation package for extremely simplifed image frame creation btw On Wednesday, March 6, 2013, Hodgess, Erin wrote: Hello again! I'm running into a problem with the image function within a loop. When I put it to the screen, it works fine, but when I put it to a PNG file, the picture doesn't change. Here is the code: for(i in 1:11) { j1 - (i-1)*n2 + 1 j2 - i*n2 predx - data.frame(pred=kr3$var1.pred[j1:j2],lon=kr3@sp@coords[,1],lat=kr3@sp@coords[,2]) coordinates(predx) - ~lon+lat proj4string(predx) - proj4string(kr3) predx.wgs84 - spTransform(predx,CRS(+proj=longlat +datum=WGS84)) llSPix - as(llSGDF, SpatialPixelsDataFrame) (st - system.time(llSPix$pred - idw(pred ~ 1, loc=predx.wgs84,newdata = llSPix, nmax = 16)$var1.pred)) cat(stuff,i,\n) f1 - paste(pred,i,a.png,sep=) #Works fine image(llSPix,pred,col=bpy.colors(),main=f1) assign(file1,f1) cat(filen,file1,\n) png(filename = file1, width = GRD.wgs84$width, height = GRD.wgs84$height,bg=transparent) par(mar = c(0, 0, 0, 0), xaxs = i, yaxs = i) image(llSPix, pred, col = bpy.colors()) dev.off() kmlsetup(obj = GRD.wgs84, kmlfile =filename, imagefile = file1, name =as.character(2000+i)) Has anyone run into this before, please? Thanks, [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Michael Sumner Hobart, Australia e-mail: [email protected]:[email protected] [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo From: Clint Bowman [[email protected]] Sent: Tuesday, March 05, 2013 2:46 PM To: Michael Sumner Cc: Hodgess, Erin; [email protected] Subject: Re: [R-sig-Geo] question about image function When using png in a loop, especially one where the loop index is an integer, I usually make use of the %02d in the file name, e.g., png(pred_%02d_a.png,...) Clint BowmanINTERNET: [email protected] Air Quality Modeler INTERNET: [email protected] Department of Ecology VOICE: (360) 407-6815 PO Box 47600FAX:(360) 407-7534 Olympia, WA 98504-7600 USPS: PO Box 47600, Olympia, WA 98504-7600 Parcels:300 Desmond Drive, Lacey, WA 98503-1274 On Wed, 6 Mar 2013, Michael Sumner wrote: I don't see filename defined, or any packages required. Is it the png files that are a problem or the 'kmlsetup' output? (whatever that is) See the animation package for extremely simplifed image frame creation btw On Wednesday, March 6, 2013, Hodgess, Erin wrote: Hello again! I'm running into a problem with the image function within a loop. When I put it to the screen, it works fine, but when I put it to a PNG file, the picture doesn't change. Here is the code: for(i in 1:11) { j1 - (i-1)*n2 + 1 j2 - i*n2 predx - data.frame(pred=kr3$var1.pred[j1:j2],lon=kr3@sp @coords[,1],lat=kr3@sp@coords[,2]) coordinates(predx) - ~lon+lat proj4string(predx) - proj4string(kr3) predx.wgs84 - spTransform(predx,CRS(+proj=longlat +datum=WGS84)) llSPix - as(llSGDF, SpatialPixelsDataFrame) (st - system.time(llSPix$pred - idw(pred ~ 1, loc=predx.wgs84,newdata = llSPix, nmax = 16)$var1.pred)) cat(stuff,i,\n) f1 - paste(pred,i,a.png,sep=) #Works fine image(llSPix,pred,col=bpy.colors(),main=f1) assign(file1,f1) cat(filen,file1,\n) png(filename = file1, width = GRD.wgs84$width, height = GRD.wgs84$height,bg=transparent) par(mar = c(0, 0, 0, 0), xaxs = i, yaxs = i) image(llSPix, pred, col = bpy.colors()) dev.off() kmlsetup(obj = GRD.wgs84, kmlfile =filename, imagefile = file1, name =as.character(2000+i)) Has anyone run into this before, please? Thanks, [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] javascript:; https://stat.ethz.ch/mailman/listinfo/r-sig-geo Michael Sumner Hobart, Australia e-mail: [email protected] [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo https://stat.ethz.ch/pipermail/r-sig-geo/2013-February/017577.html (where the same error message appeared, but in Italian!) On 02/23/2013 04:57 AM, Hodgess, Erin wrote: Hello everyone: I have an STFDF object and I would like to produce a variogramST. I'm using the version of spacetime from 2/19. However, I'm getting a weird error message. Here is the output: str(city4.st) Formal class 'STFDF' [package spacetime] with 4 slots ..@ data :'data.frame': 230 obs. of 1 variable: .. ..$ index: num [1:230] 100 100 100 100 100 100 100 100 100 100 ... ..@ sp :Formal class 'SpatialPoints' [package sp] with 3 slots .. .. ..@ coords : num [1:23, 1:2] -99.7 -101.8 -97.7 -94.1 -97.5 ... .. .. .. ..- attr(*, dimnames)=List of 2 .. .. .. .. ..$ : NULL .. .. .. .. ..$ : chr [1:2] lon lat .. .. ..@ bbox : num [1:2, 1:2] -106.5 25.9 -94 35.2 .. .. .. ..- attr(*, dimnames)=List of 2 .. .. .. .. ..$ : chr [1:2] lon lat .. .. .. .. ..$ : chr [1:2] min max .. .. ..@ proj4string:Formal class 'CRS' [package sp] with 1 slots .. .. .. .. ..@ projargs: chr +proj=lcc +lat_1=27.5 +lat_2=35 +lat_0=18 +lon_0=-100 +x_0=150 +y_0=500 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +t| __truncated__ ..@ time :An ‘xts’ object on 2001-01-01/2010-01-01 containing: Data: int [1:10, 1] 1 2 3 4 5 6 7 8 9 10 - attr(*, dimnames)=List of 2 ..$ : NULL ..$ : chr timeIndex Indexed by objects of class: [Date] TZ: UTC xts Attributes: ..@ endTime: POSIXct[1:10], format: 2002-01-01 2003-01-01 2004-01-01 2005-01-01 ... variogramST(index~1,city4.st) Error in if (nrow(x) nrow(y)) matrix(t(apply(x, 1, spDiN1, y = y, ll = longlat)), : argument is of length zero Edzer Pebesma Institute for Geoinformatics (ifgi), University of Münster Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de http://www.52north.org/geostatistics [email protected] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Erin, see https://stat.ethz.ch/**pipermail/r-sig-geo/2013-**February/017577.htmlhttps://stat.ethz.ch/pipermail/r-sig-geo/2013-February/017577.html (where the same error message appeared, but in Italian!) On 02/23/2013 04:57 AM, Hodgess, Erin wrote: Hello everyone: I have an STFDF object and I would like to produce a variogramST. I'm using the version of spacetime from 2/19. However, I'm getting a weird error message. Here is the output: str(city4.st) Formal class 'STFDF' [package spacetime] with 4 slots ..@ data :'data.frame': 230 obs. of 1 variable: .. ..$ index: num [1:230] 100 100 100 100 100 100 100 100 100 100 ... ..@ sp :Formal class 'SpatialPoints' [package sp] with 3 slots .. .. ..@ coords : num [1:23, 1:2] -99.7 -101.8 -97.7 -94.1 -97.5 .. .. .. ..- attr(*, dimnames)=List of 2 .. .. .. .. ..$ : NULL .. .. .. .. ..$ : chr [1:2] lon lat .. .. ..@ bbox : num [1:2, 1:2] -106.5 25.9 -94 35.2 .. .. .. ..- attr(*, dimnames)=List of 2 .. .. .. .. ..$ : chr [1:2] lon lat .. .. .. .. ..$ : chr [1:2] min max .. .. ..@ proj4string:Formal class 'CRS' [package sp] with 1 slots .. .. .. .. ..@ projargs: chr +proj=lcc +lat_1=27.5 +lat_2=35 +lat_0=18 +lon_0=-100 +x_0=150 +y_0=500 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +t| __truncated__ ..@ time :An xts object on 2001-01-01/2010-01-01 containing: Data: int [1:10, 1] 1 2 3 4 5 6 7 8 9 10 - attr(*, dimnames)=List of 2 ..$ : NULL ..$ : chr timeIndex Indexed by objects of class: [Date] TZ: UTC xts Attributes: ..@ endTime: POSIXct[1:10], format: 2002-01-01 2003-01-01 2004-01-01 2005-01-01 ... variogramST(index~1,city4.st) Error in if (nrow(x) nrow(y)) matrix(t(apply(x, 1, spDiN1, y = y, ll = longlat)), : argument is of length zero Any suggestions would be much appreciated. Thanks, Sincerely, [[alternative HTML version deleted]] __**_ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/**listinfo/r-sig-geohttps://stat.ethz.ch/mailman/listinfo/r-sig-geo Edzer Pebesma Institute for Geoinformatics (ifgi), University of Münster Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de http://www.52north.org/**geostatisticshttp://www.52north.org/geostatistics [email protected] __**_ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/**listinfo/r-sig-geohttps://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Formal class 'STFDF' [package spacetime] with 4 slots ..@ data :'data.frame': 8250 obs. of 1 variable: .. ..$ var1.pred: num [1:8250] 95 95.1 95.3 95.4 95.6 ... ..@ sp :Formal class 'SpatialPoints' [package sp] with 3 slots .. .. ..@ coords : num [1:825, 1:2] -108 -108 -107 -106 -106 ... .. .. .. ..- attr(*, dimnames)=List of 2 .. .. .. .. ..$ : NULL .. .. .. .. ..$ : chr [1:2] lon lat .. .. ..@ bbox : num [1:2, 1:2] -108 24 -92 36 .. .. .. ..- attr(*, dimnames)=List of 2 .. .. .. .. ..$ : chr [1:2] lon lat .. .. .. .. ..$ : chr [1:2] min max .. .. ..@ proj4string:Formal class 'CRS' [package sp] with 1 slots .. .. .. .. ..@ projargs: chr +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 ..@ time :An ‘xts’ object on 2001-01-01/2010-01-01 containing: Data: int [1:10, 1] 1 2 3 4 5 6 7 8 9 10 - attr(*, dimnames)=List of 2 ..$ : NULL ..$ : chr timeIndex Indexed by objects of class: [Date] TZ: UTC xts Attributes: ..@ endTime: POSIXct[1:10], format: 2002-01-01 2003-01-01 ... I would like to draw a wireframe, but I keep getting the following: str(pred.grd) Formal class 'STF' [package spacetime] with 3 slots ..@ sp :Formal class 'SpatialPoints' [package sp] with 3 slots .. .. ..@ coords : num [1:825, 1:2] -108 -108 -107 -106 -106 ... .. .. .. ..- attr(*, dimnames)=List of 2 .. .. .. .. ..$ : NULL .. .. .. .. ..$ : chr [1:2] lon lat .. .. ..@ bbox : num [1:2, 1:2] -108 24 -92 36 .. .. .. ..- attr(*, dimnames)=List of 2 .. .. .. .. ..$ : chr [1:2] lon lat .. .. .. .. ..$ : chr [1:2] min max .. .. ..@ proj4string:Formal class 'CRS' [package sp] with 1 slots .. .. .. .. ..@ projargs: chr +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 ..@ time :An ‘xts’ object on 2001-01-01/2010-01-01 containing: Data: int [1:10, 1] 1 2 3 4 5 6 7 8 9 10 - attr(*, dimnames)=List of 2 ..$ : NULL ..$ : chr timeIndex Indexed by objects of class: [Date] TZ: UTC xts Attributes: ..@ endTime: POSIXct[1:10], format: 2002-01-01 2003-01-01 ... wireframe(pred.grd@coords[,1],pred.grd@coords[,2],pred.grd@time,data=kr1$var1.pred) Error in wireframe(pred.grd@coords[, 1], pred.grd@coords[, 2], pred.grd@time, : no slot of name coords for this object of class STF the error message is pretty clear; showClass(STF) would have told you where to look, as does your output above: assuming kr1 and pred.grd refer to the same geometry, you should have used pred.grd@sp@coords and index(pred.grd@time). Simpler would be to use as(kr1, data.frame) and specify your wireframe from there, using the formula var1.pred~lon+lat . I don't understand what you expect lattice::wireframe to do when you pass it four vectors. Edzer Pebesma Institute for Geoinformatics (ifgi), University of Münster Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de http://www.52north.org/geostatistics [email protected] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geoRe: [R-sig-Geo] question on the difference between spdep function spautolm() and glm() with autocovariate
2013-02-12 Thread Li, Han Hey Roger, thank you very much for your explanation. I found out from R help that errorsarlm() and spautolm(..., family=SAR) do the same analysis. So I was afraid this (two functions using the same model and doing the same analysis) might happen somewhere else. I never took any regression analysis or spatial statistic class. In general I have been using your book Applied Spatial Data Analysis with R as the guide. I have this following question: You mentioned that the model for glm()+auto is y = rho W y + Xb + e Based on what I know, this is also the model for function lagsarlm() So do glm()+auto and lagsarlm() do the same/similar analysis? The second question is: If I use several different models, for example here, I used CAR, SAR, glm()+auto, and lagsarlm(), is it wise to use AIC to decide which model fits my data better? Based on your previous explanation, I understand now glm()+auto will cause problems to interpret the coefficients. But how about CAR and SAR, or lagsarlm(), assuming it is different from glm()+auto? Also after trying all these models (SAR, CAR, glm()+auto, and lagsarlm()), if some variables are significant in all of these models, can I interpret that these variables are more important than other variables which are only significant in one or two models? For example, variable 1 is significant in all 4 models for predicting the presence of bats; variable 2 is only significant in SAR. Can I say that variable 1 is more important than variable 2? I would like to point out that my research is about bats in urban environment. There is very little biological information to help me interpret the results. Thanks again. Han Li Ph.D. Candidate Department of Biology Baylor University Waco, TX 76798-7388 Phone: (254) 710-2151 Fax: (254) 710-2969 [email protected]:[email protected] From: Roger Bivand [[email protected]] Sent: Tuesday, February 12, 2013 1:05 AM To: Li, Han Cc: [email protected] Subject: Re: [R-sig-Geo] question on the difference between spdep function spautolm() and glm() with autocovariate On Tue, 12 Feb 2013, Li, Han wrote: Dear list, I am currently working on spatial autoregression modeling for my dissertation research. I want to use regression models to identify socioeconomic/landscape variables (15 total, var1$bat_survey - var15$bay_survey) that can affect the presence/absence of bats (p/a$bat_survey). Since spatial autocorrelation exists in my P/A data, I tried different spatial models. My question is: If I use the same way (same neighboring criteria, same weight style) to define neighbors, and build model #1 spatial simultaneous autoregression model (SAR) by function spautolm(), and model #2 glm() with autocovariate generated by function autocov_dist(), should I expect the same result, or not? No, because obviously they are different models: spautolm: y = Xb + u, u = lambda W u + e, e ~ N(0, sigma2 I) glm+auto: y = rho W y + Xb + e Interpreting the latter is subject to great difficulty (see ?impacts) because the DGP is: (I - rho W) y = Xb + e, y = (I - rho W)^{-1} (Xb + e) so the b coefficients cannot be interpreted directly. In addition, the glm estimate of rho is biassed because it is not constrained to its feasible range (so that (I - rho W) can be inverted). Using geo-additive models, ML in the spautolm case, or others, is easier to handle because the fitted coefficients don't interact. It is only in some settings that the observations interact with each other directly, more often the autocorrelation is in the residuals. Hope this clarifies, Roger I understood that if I use glm() with autocovariate it will include one more variable (the autocovariate) in the result. I also learned that glm() is more a predicting model and spautolm() is more an explanatory model. But I am not sure whether the significant variables selected by these two models will be same. ##r code example## model_1 - spautolm (p/a ~ var1 + var2 + ... + var15, data = bat_survey, listw = neighbor_regime1, family=SAR) autocov_model_2 - autocov_dist (p/a$bat_survey, xy = coords, style = W, type = one) model_2 - glm (p/a ~ var1 + var2 + ... + var15 + autocov_model_2, family = binomial, data = bat_survey) Thanks in advance. Your insight will be deeply appreciated. Han Li Department of Biology Baylor University Waco, TX 76798-7388 Phone: (254) 710-2151 Fax: (254) 710-2969 [email protected]:[email protected] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Roger Bivand Department of Economics, NHH Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: [email protected] R-sig-Geo mailing list [email protected] https://stat.ethz.chRe: [R-sig-Geo] question on the difference between spdep function spautolm() and glm() with autocovariate
2013-02-11 Thread Roger Bivand On Tue, 12 Feb 2013, Li, Han wrote: Dear list, I am currently working on spatial autoregression modeling for my dissertation research. I want to use regression models to identify socioeconomic/landscape variables (15 total, var1$bat_survey - var15$bay_survey) that can affect the presence/absence of bats (p/a$bat_survey). Since spatial autocorrelation exists in my P/A data, I tried different spatial models. My question is: If I use the same way (same neighboring criteria, same weight style) to define neighbors, and build model #1 spatial simultaneous autoregression model (SAR) by function spautolm(), and model #2 glm() with autocovariate generated by function autocov_dist(), should I expect the same result, or not? No, because obviously they are different models: spautolm: y = Xb + u, u = lambda W u + e, e ~ N(0, sigma2 I) glm+auto: y = rho W y + Xb + e Interpreting the latter is subject to great difficulty (see ?impacts) because the DGP is: (I - rho W) y = Xb + e, y = (I - rho W)^{-1} (Xb + e) so the b coefficients cannot be interpreted directly. In addition, the glm estimate of rho is biassed because it is not constrained to its feasible range (so that (I - rho W) can be inverted). Using geo-additive models, ML in the spautolm case, or others, is easier to handle because the fitted coefficients don't interact. It is only in some settings that the observations interact with each other directly, more often the autocorrelation is in the residuals. Hope this clarifies, Roger I understood that if I use glm() with autocovariate it will include one more variable (the autocovariate) in the result. I also learned that glm() is more a predicting model and spautolm() is more an explanatory model. But I am not sure whether the significant variables selected by these two models will be same. ##r code example## model_1 - spautolm (p/a ~ var1 + var2 + ... + var15, data = bat_survey, listw = neighbor_regime1, family=SAR) autocov_model_2 - autocov_dist (p/a$bat_survey, xy = coords, style = W, type = one) model_2 - glm (p/a ~ var1 + var2 + ... + var15 + autocov_model_2, family = binomial, data = bat_survey) Thanks in advance. Your insight will be deeply appreciated. Han Li Department of Biology Baylor University Waco, TX 76798-7388 Phone: (254) 710-2151 Fax: (254) 710-2969 [email protected]:[email protected] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Roger Bivand Department of Economics, NHH Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: [email protected] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo From: [email protected] [[email protected]] on behalf of Hodgess, Erin Sent: Saturday, February 09, 2013 10:10 PM To: [email protected] Subject: [R-sig-Geo] question about STIDF objects Dear R-sig Geo People: I have a question about the space-time package, please. I'm working through the worked examples from the JRSS paper. My question is: can I make the wind.data object an STIDF instead of an STFDF, please? I would like it in that form so I can use plotKML to produce those cool Google Earth KML files. Thanks, [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Just checking - are you trying to install an up-to-date package onto a 2.14, or are you compiling from the archives? Cheers, Roman On Sat, Jan 12, 2013 at 3:15 PM, Hugo Loyola [email protected]: I am trying to install the rgdal package on R that runs under Ubuntu 12.04. Everything goes smooth until this happen: installing to /home/rosa/R/i686-pc-linux-gnu-library/2.14/rgdal/libs ** inst ** preparing package for lazy loading ** help *** installing help indices ** building package indices ... *** tangling vignette sources ... âOGR_shape_encoding.Rnwâ using âUTF-8â ** testing if installed package can be loaded *** caught illegal operation *** address 0xb6001b34, cause 'illegal operand' Traceback: 1: dyn.load(file, DLLpath = DLLpath, ...) 2: library.dynam(lib, package, package.lib) 3: loadNamespace(package, c(which.lib.loc, lib.loc)) 4: doTryCatch(return(expr), name, parentenv, handler) 5: tryCatchOne(expr, names, parentenv, handlers[[1L]]) 6: tryCatchList(expr, classes, parentenv, handlers) Illegal instruction (core dumped) ERROR: loading failed * removing â/home/rosa/R/i686-pc-linux-gnu-library/2.14/rgdalâ The downloaded packages are in â/tmp/RtmpwQ513H/downloaded_packagesâ Warning message: In install.packages(rgdal) : installation of package ârgdalâ had non-zero exit status What is possibly the cause of this fail? Thank you, in advance [[alternative HTML version deleted]] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo there are options that come fairly close: What you could do is add columns to a SpatialPixelsDataFrame: a$colb = b$col a$colc = c$col but of course you need to check yourself that this make sense (i.e. the pixels of a, b and c are identical). Another way would be working with SpatialGridDataFrame's: fullgrid(a) = TRUE fullgrid(b) = TRUE fullgrid(c) = TRUE abc = cbind(a,b,c) which combines layers into one object (similar to what raster::stack does), and does check for equality of all grid topologies. Best regards, On 12/04/2012 11:34 AM, Michela Giusti wrote: Hello everybody, I have many SpatialPixelsDataFrame in my R project and I am wondering if there is a command like the command stack for raster to do the same with this object class. Can you help me? Thank you. Michela Edzer Pebesma Institute for Geoinformatics (ifgi), University of Münster Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de http://www.52north.org/geostatistics [email protected] R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Good, please post your question on the list after checking all the references, for example for lm.morantest, lm.morantest.sad and lm.LMtests; lm.moantest.exact is covered in Bivand, Roger; Müller, Werner G.; Reder, Markus. Power calculations for global and local Moran's l. Computational Statistics Data Analysis 2009 ;Volum 53.(8) s. 2859-2872. The Tiefelsdorf 2002 reference is very relevant, as are works cited there. In general, normality is not a real issue, unless the shape of the residuals is very badly behaved, but you can always run some simulations. Roger On Fri, 16 Nov 2012, Marcio Pupin Mello wrote: Dear Prof. Roger, thanks for your reply. I've just sent a new post to the list and it worked now. I haven't read the papers where tests were introduced. Would them be those ones you cited in the help for summary.sarlm {spdep}? But I read in LeSage and Pace (2009) that a CI for the estimated parameters could be found through simulation processes. Apparently, results of the summary.sarlm showed that the test for the Intercept is Z based (Gaussian assumption, though), but I am not sure about the LR and Wald tests for Rho as well as LM test for residual autocorrelation. Bests, Marcio www.dsr.inpe.br/~mello LeSage, J. and Pace, R.K. Introduction to Spatial Econometrics. CRC Press: Boca Raton, USA, 2009. 354 p. On 12/22/11 1:21 PM, Roger Bivand wrote: On Thu, 22 Dec 2011, Dimitrios Efthymiou wrote: Hi Roger, thanks for your help. Now I get a result, but many variables that are significant at the error and lag models, are not significant here. In addition, few variables have statistically significant direct impact, others indirect and others total. Simulated z-values: DirectIndirect Total log(sqm) 134.33514612 6.95584281 27.94811234 parking 9.01409655 0.51593572 1.58416761 chimney 4.20407856 1.15585250 1.67326910 auto_heat5.37231195 0.66981923 1.30766715 klima4.75291342 2.42475203 2.99642118 dist_cbd 1.86144465 -2.04747304 -1.08141630 metro -0.29541702 2.75180894 3.47650020 This question does not have to do with the package's functions, but I would appreciate any help because I'm not totally sure how to explain these results. Please refer to the advice given by J. Paul Elhorst on the openspace list: http://groups.google.com/group/openspace-list/browse_thread/thread/9381788c8f5c5352 (click on: Show quoted text) In addition, if you look at the summary output for your error model with Hausman=TRUE, you'll see whether an error model is acceptable if the Common Factor LR test indicates error rather than spatial Durbin. The Hausman test is described in LeSage Pace 2009. Hope this helps, Roger I'm using the 'spdep' package of R, in a server with Ubuntu. I'm trying to estimate a lagsarlm mixed model using 8.000 observations house prices. until now, I can estimate the errorraslm, lagsarlm (lagged) and sacsarlm without problem, but when I try to estimate the mixed lagsarlm, NAs appear at the Std. Error, z-value and p columns of many variables. This means that you have issues in the close correlations between the variables and lagged variables. If you use a trs= argument of MC traces, the finite difference estimate of the Hessian will be augmented as shown in LeSage Pace (2009), and as the help page for lagsarlm says, the problem may go away. You need the series of traces of powers of W anyway for impacts(), so no extra effort is required. If the coefficient itself is dropped, a variable and its lag are aliased, but this isn't your case here. Hope this helps, Roger Does anyone know what causes this? I'm using k-nearest neighbours = 12, the S-coding variance - stabilising, and the LU method (I concluded to these numbers and methods after experimentation in the other models). Characteristics of weights list object: Neighbour list object: Number of regions: 8066 Number of nonzero links: 185518 Percentage nonzero weights: 0.2851475 Average number of links: 23 Non-symmetric neighbours list Weights style: S Weights constants summary: n nn S0 S1 S2 S 8066 65060356 8066 630.5066 32982.99 other attached packages: [1] ape_2.8 ade4_1.4-17 car_2.0-11 [4] survival_2.36-10nnet_7.3-1 gstat_1.0-10 [7] spacetime_0.5-7 xts_0.8-2 zoo_1.7-6 [10] rasterVis_0.10-7hexbin_1.26.0Re: [R-sig-Geo] Question about space-time analysis routines
2012-11-05 Thread Virgilio Gómez-Rubio I am curious if the list can point me to some other libraries in R, besides INLA, that are capable of using the spacetime data directly in a regression INLA is a good software to fit spatiotemporal models as it is very flexible. framework, or perhaps in a spatio-temporal disease clustering framework. I know about the DCluster library and use it some, but would be very grateful for any other ideas. You may want to check DClusterm : https://r-forge.r-project.org/projects/dclusterm/ It is underdevelopment but it implements space-time clustering and spatial clustering based on GLMs. Best, Virgilio R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo I am curious if the list can point me to some other libraries in R, besides INLA, that are capable of using the spacetime data directly in a regression INLA is a good software to fit spatiotemporal models as it is very flexible. framework, or perhaps in a spatio-temporal disease clustering framework. I know about the DCluster library and use it some, but would be very grateful for any other ideas. You may want to check DClusterm : https://r-forge.r-project.org/projects/dclusterm/ It is underdevelopment but it implements space-time clustering and spatial clustering based on GLMs. Best, Virgilio R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo I am curious if the list can point me to some other libraries in R, besides INLA, that are capable of using the spacetime data directly in a regression INLA is a good software to fit spatiotemporal models as it is very flexible. framework, or perhaps in a spatio-temporal disease clustering framework. I know about the DCluster library and use it some, but would be very grateful for any other ideas. You may want to check DClusterm : https://r-forge.r-project.org/projects/dclusterm/ It is underdevelopment but it implements space-time clustering and spatial clustering based on GLMs. Best, Virgilio R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo Corey Sparks, PhD Assistant professor Department of Demography The University of Texas at San Antonio 501 West Cesar E Chavez Blvd San Antonio TX 78207 Corey.sparks 'at' utsa.edu 210 458 3166 Latitude: 29.423614 / Longitude: -98.504282 R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo1 2 >
推荐文章