Re: [R-sig-Geo] Error reading hdf files with rgdal

[Kelly McCaffrey]

Here are the startup messages and the sessionInfo() output:

> library(raster)
Loading required package: sp
> library(rgdal)
rgdal: version: 1.2-5, (SVN revision 648)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
 Path to GDAL shared files:
C:/Users/kmccaffrey2011/Documents/R/win-library/3.3/rgdal/gdal
 Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
 Path to PROJ.4 shared files:
C:/Users/kmccaffrey2011/Documents/R/win-library/3.3/rgdal/proj
 Linking to sp version: 1.2-4
> library(gdalUtils)
> sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] stats graphics  grDevices utils
[5] datasets  methods   base

other attached packages:
[1] gdalUtils_2.0.1.7 rgdal_1.2-5
[3] raster_2.5-8  sp_1.2-4

loaded via a namespace (and not attached):
 [1] tools_3.3.2   Rcpp_0.12.8
 [3] R.methodsS3_1.7.1 codetools_0.2-15
 [5] grid_3.3.2iterators_1.0.8
 [7] foreach_1.4.3 R.utils_2.5.0
 [9] R.oo_1.21.0   lattice_0.20-34


Kelly

On Fri, Jan 13, 2017 at 3:06 PM, Roger Bivand  wrote:

> Please provide the output of the startup messages from the loaded
> packages, and of sessionInfo().
>
> Roger Bivand
> Norwegian School of Economics
> Bergen, Norway
>
> Fra: Kelly McCaffrey
> Sendt: fredag 13. januar, 20.51
> Emne: [R-sig-Geo] Error reading hdf files with rgdal
> Til: R-sig-geo mailing list
>
> Hello all, I'm working with some .hdf files containing MODIS satellite
> data, but I have been unable to read them into R using rgdal and gdalUtils
> since the last rgdal update. Has anyone else experienced these issues and
> does anyone have a suggested fix? The code I'm attempting to use and error
> message are as follows: library(raster) library(rgdal) library(gdalUtils)
> file
>



-- 

*Kelly McCaffrey*
Graduate Student
Department of Biological Sciences
Florida Institute of Technology
150 W. University Blvd.
Melbourne, FL, 32901

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] Error reading hdf files with rgdal

[Roger Bivand]

Please provide the output of the startup messages from the loaded packages, and 
of sessionInfo().

Roger Bivand
Norwegian School of Economics
Bergen, Norway

Fra: Kelly McCaffrey
Sendt: fredag 13. januar, 20.51
Emne: [R-sig-Geo] Error reading hdf files with rgdal
Til: R-sig-geo mailing list

Hello all, I'm working with some .hdf files containing MODIS satellite data, 
but I have been unable to read them into R using rgdal and gdalUtils since the 
last rgdal update. Has anyone else experienced these issues and does anyone 
have a suggested fix? The code I'm attempting to use and error message are as 
follows: library(raster) library(rgdal) library(gdalUtils) file

[[alternative HTML version deleted]]

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

[R-sig-Geo] Error reading hdf files with rgdal

[Kelly McCaffrey]

Hello all,
I'm working with some .hdf files containing MODIS satellite data, but I
have been unable to read them into R using rgdal and gdalUtils since the
last rgdal update. Has anyone else experienced these issues and does anyone
have a suggested fix? The code I'm attempting to use and error message are
as follows:

library(raster)
library(rgdal)
library(gdalUtils)
file<-"MODWCW_2017010_DAILY_AQUA_CHLORA_EC_1KM.hdf"
sds<-get_subdatasets(file)
Error in split1[[1]] : subscript out of bounds

I've tried this code with .hdf files from different sources, and they
return the same error. I've attempted to download and use an older version
of rgdal as a workaround, but I am having trouble obtaining its software
dependencies and would prefer a fix compatible with the current version of
the package.
Thanks!
Kelly

-- 

*Kelly McCaffrey*
Graduate Student
Department of Biological Sciences
Florida Institute of Technology
150 W. University Blvd.
Melbourne, FL, 32901

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] MODIS package released to CRAN

[Forrest Stevens]

Great stuff to finally see it make it!  Well done to all involved!

On Fri, Jan 13, 2017 at 8:24 AM Matteo Mattiuzzi 
wrote:

> Dear list,
>
>
> We would like to announce the release of the MODIS package to CRAN. This
> was finally possible thanks to the constant efforts of Florian Detsch, whom
> I would like to say thank you!
>
>
> Some months ago we decided to move the development from R-Forge to GitHub
> and we are thankful for any contributions in the form of ideas,
> bug-reports, pull requests or whatever helps.
>
>
> GitHub repository: https://github.com/MatMatt/MODIS
>
>
>
> We hope you will enjoy!
>
>
> Kind regards,
>
> Matteo
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] Elevation data

[Miluji Sb]

Dear Loic,

Thank you for your reply!

Re. aggregation, I agree - I will use either mean or median.

gcp_grid are points as below

structure(list(longitude = c(-179L, -178L, -177L, -177L, -177L,
-176L), latitude = c(-15L, -15L, -14L, 51L, 52L, -22L)), .Names =
c("longitude",
"latitude"), row.names = c("1", "2", "3", "4", "5", "6"), class =
"data.frame")

Thank you for your suggestion on the additional options, I will look them
up. Would you say I am on the right path after the changes? Thanks again!

Sincerely,

Milu

On Fri, Jan 13, 2017 at 3:35 PM, Loïc Dutrieux  wrote:

>
>
> On 13/01/2017 10:59, Miluji Sb wrote:
> > Thank you for your reply. This is what I did:
> >
> > ###
> > library(data.table)
> > library(raster)
> > library(rgdal)
> > library(foreign)
> >
> > elevation_world <- getData('worldclim', var='alt', res=2.5)
> >
> > # Aggregate Elevation to 1 degree
> > elevation_world_1deg <- aggregate(elevation_world, fact = 24, fun = sum)
>
> Aggregating with fun=sum is a bit strange for elevation. mean or median
> would certainly be a better choice.
>
> >
> > # Extract by lon lat (1° x 1° - gcp_grid)
> > elevation <- cbind(gcp_grid, alt = extract(elevation_world_1deg,
> gcp_grid))
>
> gcp_grid are points or polygons? If they are polygons, there's no need
> for the aggregation step above. Also have a look at df = TRUE, which
> returns a dataframe, and sp=TRUE, which returns a Spatial*DataFrame with
> extracted values cbinded to the attributes of the original
> Spatial*DataFrame. (both are arguments of raster::extract)
>
> Cheers,
> Loïc
>
> >
> > elevation <- as.data.frame(elevation)
> > ###
> >
> > Is this correct? Thanks again!
> >
> > Sincerely,
> >
> > Milu
> >
> > On Fri, Jan 13, 2017 at 3:11 AM, Bacou, Melanie  wrote:
> >
> >> R raster::getData("SRTM", ...) will return elevation rasters at 90m
> >> resolution.
> >> See:
> >> https://www.rdocumentation.org/packages/raster/versions/2.5-
> >> 8/topics/getData
> >> http://www.cgiar-csi.org/data/srtm-90m-digital-elevation-database-v4-1
> >>
> >> --Mel.
> >>
> >>
> >> On 1/11/2017 6:26 AM, Miluji Sb wrote:
> >>
> >>> Dear Michael,
> >>>
> >>> Thank you for your reply and the suggestions.
> >>>
> >>> Ideally, I would like a raster from which I can extract elevation at
> 1° x
> >>> 1° resolution. I do not have much experience with working with DEM but
> >>> have
> >>> work with data such as GPW.
> >>>
> >>> I will definitely look at the datasets. Could you kindly suggest one
> that
> >>> I
> >>> could convert to raster and extract? Hope that's not a silly question.
> >>>
> >>> Sincerely,
> >>>
> >>> Milu
> >>>
> >>> On Wed, Jan 11, 2017 at 1:51 AM, Michael Sumner 
> >>> wrote:
> >>>
> >>> Passed through? Maybe you want ?raster::extract
> 
>  There are a few versions of global elevation on CRAN, necessarily at
> low
>  resolution but no overall summary afaik (someone should do this :).
> 
>  This is one: https://cloud.r-project.org/
> web/packages/GEOmap/index.html
> 
>  If you have the stomach for development versions of packages see
> elevatr:
>  https://github.com/jhollist/elevatr
> 
>  I tend to have the high-resolution files at hand because we use them
>  constantly, the main ones are Gebco14/Gebco08 and Etopo1/Etopo2 (from
>  Smith-Sandwell).
> 
>  There's a reasonable overview here, you probably should find a
> specific
>  data set that is at the resolution you are after already, and you can
>  cite
>  its derivation for your work:
> 
>  http://vterrain.org/Elevation/global.html
> 
>  Cheers, Mike.
> 
>  On Wed, 11 Jan 2017 at 10:20 Miluji Sb  wrote:
> 
>  Dear all.
> >
> > Is there a way to download global elevation data at the 1° x 1°
> > resolution
> > in R using a given set of coordinates?
> >
> > I know about the getData() function but can many coordinates be
> passed
> > through this? Thanks!
> >
> > Sincerely,
> >
> > Milu
> >
> >  [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> >
> > University of Tasmania Electronic Communications Policy (December,
> > 2014).
> > This email is confidential, and is for the intended recipient only.
> > Access, disclosure, copying, distribution, or reliance on any of it
> by
> > anyone outside the intended recipient organisation is prohibited and
> > may be
> > a criminal offence. Please delete if obtained in error and email
> > confirmation to the sender. The views expressed in this email are not
> > necessarily the views of the University of Tasmania, unless clearly
> > intended otherwise.
> > .
> >
> > --
>  Dr. Michael 

Re: [R-sig-Geo] Elevation data

[Loïc Dutrieux]

On 13/01/2017 10:59, Miluji Sb wrote:
> Thank you for your reply. This is what I did:
> 
> ###
> library(data.table)
> library(raster)
> library(rgdal)
> library(foreign)
> 
> elevation_world <- getData('worldclim', var='alt', res=2.5)
> 
> # Aggregate Elevation to 1 degree
> elevation_world_1deg <- aggregate(elevation_world, fact = 24, fun = sum)

Aggregating with fun=sum is a bit strange for elevation. mean or median
would certainly be a better choice.

> 
> # Extract by lon lat (1° x 1° - gcp_grid)
> elevation <- cbind(gcp_grid, alt = extract(elevation_world_1deg, gcp_grid))

gcp_grid are points or polygons? If they are polygons, there's no need
for the aggregation step above. Also have a look at df = TRUE, which
returns a dataframe, and sp=TRUE, which returns a Spatial*DataFrame with
extracted values cbinded to the attributes of the original
Spatial*DataFrame. (both are arguments of raster::extract)

Cheers,
Loïc

> 
> elevation <- as.data.frame(elevation)
> ###
> 
> Is this correct? Thanks again!
> 
> Sincerely,
> 
> Milu
> 
> On Fri, Jan 13, 2017 at 3:11 AM, Bacou, Melanie  wrote:
> 
>> R raster::getData("SRTM", ...) will return elevation rasters at 90m
>> resolution.
>> See:
>> https://www.rdocumentation.org/packages/raster/versions/2.5-
>> 8/topics/getData
>> http://www.cgiar-csi.org/data/srtm-90m-digital-elevation-database-v4-1
>>
>> --Mel.
>>
>>
>> On 1/11/2017 6:26 AM, Miluji Sb wrote:
>>
>>> Dear Michael,
>>>
>>> Thank you for your reply and the suggestions.
>>>
>>> Ideally, I would like a raster from which I can extract elevation at 1° x
>>> 1° resolution. I do not have much experience with working with DEM but
>>> have
>>> work with data such as GPW.
>>>
>>> I will definitely look at the datasets. Could you kindly suggest one that
>>> I
>>> could convert to raster and extract? Hope that's not a silly question.
>>>
>>> Sincerely,
>>>
>>> Milu
>>>
>>> On Wed, Jan 11, 2017 at 1:51 AM, Michael Sumner 
>>> wrote:
>>>
>>> Passed through? Maybe you want ?raster::extract

 There are a few versions of global elevation on CRAN, necessarily at low
 resolution but no overall summary afaik (someone should do this :).

 This is one: https://cloud.r-project.org/web/packages/GEOmap/index.html

 If you have the stomach for development versions of packages see elevatr:
 https://github.com/jhollist/elevatr

 I tend to have the high-resolution files at hand because we use them
 constantly, the main ones are Gebco14/Gebco08 and Etopo1/Etopo2 (from
 Smith-Sandwell).

 There's a reasonable overview here, you probably should find a specific
 data set that is at the resolution you are after already, and you can
 cite
 its derivation for your work:

 http://vterrain.org/Elevation/global.html

 Cheers, Mike.

 On Wed, 11 Jan 2017 at 10:20 Miluji Sb  wrote:

 Dear all.
>
> Is there a way to download global elevation data at the 1° x 1°
> resolution
> in R using a given set of coordinates?
>
> I know about the getData() function but can many coordinates be passed
> through this? Thanks!
>
> Sincerely,
>
> Milu
>
>  [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
> University of Tasmania Electronic Communications Policy (December,
> 2014).
> This email is confidential, and is for the intended recipient only.
> Access, disclosure, copying, distribution, or reliance on any of it by
> anyone outside the intended recipient organisation is prohibited and
> may be
> a criminal offence. Please delete if obtained in error and email
> confirmation to the sender. The views expressed in this email are not
> necessarily the views of the University of Tasmania, unless clearly
> intended otherwise.
> .
>
> --
 Dr. Michael Sumner
 Software and Database Engineer
 Australian Antarctic Division
 203 Channel Highway
 Kingston Tasmania 7050 Australia


 [[alternative HTML version deleted]]
>>>
>>> ___
>>> R-sig-Geo mailing list
>>> R-sig-Geo@r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>>
> 
>   [[alternative HTML version deleted]]
> 
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

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

[R-sig-Geo] MODIS package released to CRAN

[Matteo Mattiuzzi]

Dear list,


We would like to announce the release of the MODIS package to CRAN. This
was finally possible thanks to the constant efforts of Florian Detsch, whom
I would like to say thank you!


Some months ago we decided to move the development from R-Forge to GitHub
and we are thankful for any contributions in the form of ideas,
bug-reports, pull requests or whatever helps.


GitHub repository: https://github.com/MatMatt/MODIS



We hope you will enjoy!


Kind regards,

Matteo

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] shortestPath in a loop

[marta azores]

Thanks for your quick reply to my question,

This ten points are ordered by time (date + time), and I know ,the points
are all connected.  I didn't include the column with the time in this
example, to make it simple.

After read your messages, I was looking the minimum spanning tree for the
vertices.
The advantage with the function "shortestPath" is that you have as result a
spatial line output, which I need it.
However, the minimum spanning tree give as output an ordiplot.

Is it possible to transform the ordiplot into spatialLines?

Thanks for your answers

Marta



2017-01-12 20:03 GMT-01:00 Eric Carr :

> Generically independent of R, A graph and the connectivity between points
> needs defined before a shortest path algorithm can be applied.  If it
> assumes all points are connected, than the shortest path will be a straight
> line.  What you are looking for is some sort of  minimum spanning tree for
> the vertices.
> Eric
>
> On Thu, Jan 12, 2017 at 11:17 AM, marta azores 
> wrote:
>
>> Dear forum members,
>>
>> I would like to know how join several points with the aim to track a ship.
>>
>> After reading the documentation of some packages, I decided to use the
>> function shortestPath, but I only got the line between the first and the
>> last location of my points list. I need the complete survey, including also
>> the middle points. I try a loop to build the survey of the boats using
>> their locations, but It didn't work to me.
>>
>>
>> Any idea?
>>
>> Thanks in advance,
>>
>> Marta
>>
>> #script# it's also attached in a R.file: question loop2.R
>> ##
>> #
>> #raster# it's attached
>> azoTS1<- raster("C:/Users/Documents/azoTS1.tif")#wgs84
>> #
>> #10 points# it's attached
>> boat <- read.table("C:/Users/Documents/10pontos.csv", header=TRUE,
>> sep=",", na.strings="NA", dec=".", strip.white=TRUE)#
>> head(boat)
>>
>> #raster to transitionlayer
>> trCostS4<- transition(1/azoTS1, mean, directions=4)
>>
>> # points to spatialpointsdataframe
>> x=boat$Long1
>> y=boat$Lat1
>> coords = cbind(x, y)
>> plot(coords)
>> sp = SpatialPoints(coords, proj4string=CRS("+proj=longlat +ellps=WGS84
>> +datum=WGS84"), bbox = NULL)
>> sp
>> spdf=SpatialPointsDataFrame(sp,boat)
>> spdf
>> nrow(spdf)
>> plot(sp,axes=TRUE)
>> plot(spdf,add=TRUE, axes=TRUE)
>>
>> #shortestpath
>>
>> ## 1) this script only join the first point of the list and the last one,
>> and the points in the middle are not used.
>> CostpathSPdf <- shortestPath(trCostS4, spdf[1,], spdf[10,],
>> output="SpatialLines")
>> plot(CostpathSPdf,add=TRUE,axes=TRUE,col=2)#R_plot1.png (it's attached)
>>
>> ## 2) this script didn't work to me
>>
>> #first way from website: http://stackoverflow.com/quest
>> ions/8127066/loop-or-sapply-function-for-multiple-least-
>> cost-analysis-in-r?answertab=active#tab-top
>> for(i in 1:nrow(spdf)) {
>>   # Computation
>>   Costpath <- shortestPath(trCostS4, spdf[i,], spdf[10,],
>> output="SpatialLines")
>>   plot(Costpath)
>>
>> }
>>
>> #Error in validObject(.Object) :
>> #invalid class “SpatialLines” object: bbox should never contain infinite
>> values
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] Elevation data

[Miluji Sb]

Thank you for your reply. This is what I did:

###
library(data.table)
library(raster)
library(rgdal)
library(foreign)

elevation_world <- getData('worldclim', var='alt', res=2.5)

# Aggregate Elevation to 1 degree
elevation_world_1deg <- aggregate(elevation_world, fact = 24, fun = sum)

# Extract by lon lat (1° x 1° - gcp_grid)
elevation <- cbind(gcp_grid, alt = extract(elevation_world_1deg, gcp_grid))

elevation <- as.data.frame(elevation)
###

Is this correct? Thanks again!

Sincerely,

Milu

On Fri, Jan 13, 2017 at 3:11 AM, Bacou, Melanie  wrote:

> R raster::getData("SRTM", ...) will return elevation rasters at 90m
> resolution.
> See:
> https://www.rdocumentation.org/packages/raster/versions/2.5-
> 8/topics/getData
> http://www.cgiar-csi.org/data/srtm-90m-digital-elevation-database-v4-1
>
> --Mel.
>
>
> On 1/11/2017 6:26 AM, Miluji Sb wrote:
>
>> Dear Michael,
>>
>> Thank you for your reply and the suggestions.
>>
>> Ideally, I would like a raster from which I can extract elevation at 1° x
>> 1° resolution. I do not have much experience with working with DEM but
>> have
>> work with data such as GPW.
>>
>> I will definitely look at the datasets. Could you kindly suggest one that
>> I
>> could convert to raster and extract? Hope that's not a silly question.
>>
>> Sincerely,
>>
>> Milu
>>
>> On Wed, Jan 11, 2017 at 1:51 AM, Michael Sumner 
>> wrote:
>>
>> Passed through? Maybe you want ?raster::extract
>>>
>>> There are a few versions of global elevation on CRAN, necessarily at low
>>> resolution but no overall summary afaik (someone should do this :).
>>>
>>> This is one: https://cloud.r-project.org/web/packages/GEOmap/index.html
>>>
>>> If you have the stomach for development versions of packages see elevatr:
>>> https://github.com/jhollist/elevatr
>>>
>>> I tend to have the high-resolution files at hand because we use them
>>> constantly, the main ones are Gebco14/Gebco08 and Etopo1/Etopo2 (from
>>> Smith-Sandwell).
>>>
>>> There's a reasonable overview here, you probably should find a specific
>>> data set that is at the resolution you are after already, and you can
>>> cite
>>> its derivation for your work:
>>>
>>> http://vterrain.org/Elevation/global.html
>>>
>>> Cheers, Mike.
>>>
>>> On Wed, 11 Jan 2017 at 10:20 Miluji Sb  wrote:
>>>
>>> Dear all.

 Is there a way to download global elevation data at the 1° x 1°
 resolution
 in R using a given set of coordinates?

 I know about the getData() function but can many coordinates be passed
 through this? Thanks!

 Sincerely,

 Milu

  [[alternative HTML version deleted]]

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


 University of Tasmania Electronic Communications Policy (December,
 2014).
 This email is confidential, and is for the intended recipient only.
 Access, disclosure, copying, distribution, or reliance on any of it by
 anyone outside the intended recipient organisation is prohibited and
 may be
 a criminal offence. Please delete if obtained in error and email
 confirmation to the sender. The views expressed in this email are not
 necessarily the views of the University of Tasmania, unless clearly
 intended otherwise.
 .

 --
>>> Dr. Michael Sumner
>>> Software and Database Engineer
>>> Australian Antarctic Division
>>> 203 Channel Highway
>>> Kingston Tasmania 7050 Australia
>>>
>>>
>>> [[alternative HTML version deleted]]
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] sp::spDists*

[Edzer Pebesma]

A new version of gstat, fixing the same bug, is now also on CRAN.

On 09/01/17 11:39, Roger Bivand wrote:
> New versions of spdep and spgwr fixing the same bug are now on CRAN.
> 
> Roger
> 
> On Thu, 8 Dec 2016, Edzer Pebesma wrote:
> 
>> While building support for geosphere distance functions in sf and
>> comparing results from geosphere::distMeeuws with those of sp::spDists,
>> I found a typo in the C code underlying the great circle distance
>> computation functions in sp (sp::spDists, sp::spDistsN1) [*]:
>>
>> https://github.com/edzer/sp/commit/d8374ff7efc6735cba9a054748c602bed0672f23
>>
>>
>> Before:
>>> library(sf)
>> Linking to GEOS 3.5.0, GDAL 2.1.0, proj.4 4.9.2
>>> library(sp)
>>> library(units)
>>>
>>> x = st_sfc(
>> + st_point(c(0,0)),
>> + st_point(c(1,0)),
>> + st_point(c(2,0)),
>> + st_point(c(3,0)),
>> + crs = 4326
>> + )
>>>
>>> y = st_sfc(
>> + st_point(c(0,10)),
>> + st_point(c(1,0)),
>> + st_point(c(2,0)),
>> + st_point(c(3,0)),
>> + st_point(c(4,0)),
>> + crs = 4326
>> + )
>>>
>>> st_distance(x, y)
>> Units: m
>>[,1] [,2] [,3] [,4] [,5]
>> [1,] 1105855 111319.5 222639.0 333958.5 445278.0
>> [2,] 387  0.0 111319.5 222639.0 333958.5
>> [3,] 1127822 111319.5  0.0 111319.5 222639.0
>> [4,] 1154693 222639.0 111319.5  0.0 111319.5
>>>
>>> d.sf = st_distance(x, y, geosphere::distMeeus)
>>> d.sp = spDists(as(x, "Spatial"), as(y, "Spatial"))
>>> units(d.sp) = make_unit("km")
>>> d.sf - d.sp
>> Units: m
>> [,1] [,2] [,3] [,4] [,5]
>> [1,] 1851.990 0.00e+00 0.00e+00 0.00e+000
>> [2,] 1842.938 0.00e+00 0.00e+00 2.910383e-110
>> [3,] 1816.564 0.00e+00 0.00e+00 1.455192e-110
>> [4,] 1775.032 2.910383e-11 1.455192e-11 0.00e+000
>>>
>>
>>
>>
>> After:
>>
>>> library(sf)
>> Linking to GEOS 3.5.0, GDAL 2.1.0, proj.4 4.9.2
>>> library(sp)
>>> library(units)
>>>
>>> x = st_sfc(
>> + st_point(c(0,0)),
>> + st_point(c(1,0)),
>> + st_point(c(2,0)),
>> + st_point(c(3,0)),
>> + crs = 4326
>> + )
>>>
>>> y = st_sfc(
>> + st_point(c(0,10)),
>> + st_point(c(1,0)),
>> + st_point(c(2,0)),
>> + st_point(c(3,0)),
>> + st_point(c(4,0)),
>> + crs = 4326
>> + )
>>>
>>> st_distance(x, y)
>> Units: m
>>[,1] [,2] [,3] [,4] [,5]
>> [1,] 1105855 111319.5 222639.0 333958.5 445278.0
>> [2,] 387  0.0 111319.5 222639.0 333958.5
>> [3,] 1127822 111319.5  0.0 111319.5 222639.0
>> [4,] 1154693 222639.0 111319.5  0.0 111319.5
>>>
>>> d.sf = st_distance(x, y, geosphere::distMeeus)
>>> d.sp = spDists(as(x, "Spatial"), as(y, "Spatial"))
>>> units(d.sp) = make_unit("km")
>>> d.sf - d.sp
>> Units: m
>>  [,1] [,2] [,3] [,4] [,5]
>> [1,] -2.328306e-10 0.00e+00 0.00e+00 0.00e+000
>> [2,]  2.328306e-10 0.00e+00 0.00e+00 2.910383e-110
>> [3,]  0.00e+00 0.00e+00 0.00e+00 1.455192e-110
>> [4,]  0.00e+00 2.910383e-11 1.455192e-11 0.00e+000
>>>
>>
>> [*] assuming the source code of
>> http://www.abecedarical.com/javascript/script_greatcircle.html is
>> correct; I don't have access to the Meeuws book mentioned there under
>> (1).
>>
>> The same C function is used by gstat; I corrected that as well.
>>
> 

-- 
Edzer Pebesma
Institute for Geoinformatics  (ifgi),  University of Münster
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers & Geosciences:   http://elsevier.com/locate/cageo/



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