Subscribe to R-sig-geo feed
This is an archive for R-sig-geo, not a forum. It is not possible to post through Nabble - you may not start a new thread nor follow up an existing thread. If you wish to post, but are not subscribed to the list through the list homepage, subscribe first (through the list homepage, not via Nabble) and post from your subscribed email address. Until 2015-06-20, subscribers could post through Nabble, but policy has been changed as too many non-subscribers misunderstood the interface.
Updated: 2 hours 1 min ago

Re: Maintain SRID with st_write to postgis db

Wed, 03/15/2017 - 15:43
Hi Chris,

Thanks for that insight - I hadn't caught that part of the vignette (it had
admittedly been a bit since I last looked at it), and I appreciate your
interpretation. What you say makes sense to me, reading through that, and I
can naturally add code to set the SRID in PostGIS pretty easily.

Best,
Mike

On Wed, Mar 15, 2017 at 5:21 PM, chris english <
[hidden email]> wrote:

> Michael,
>
> this from the sf vignettes on read write to databases:
>
> The dsn and layer arguments to st_read and st_write denote a data source
> name and optionally a layer name. Their exact interpretation as well as the
> options they support vary per driver, the GDAL driver documentation is best
> consulted for this. For instance, a PostGIS table in database postgis might
> be read by
>
> meuse <- st_read("PG:dbname=postgis", "meuse")
>
> where the PG: string indicates this concerns the PostGIS driver, followed
> by database name, and possibly port and user credentials.
>
> st_read typically reads the coordinate reference system as proj4string,
> but not the EPSG (SRID). GDAL cannot retrieve SRID (EPSG code) from
> proj4string strings, and, when needed, it has to be set by the user.
>
> I just re-read this this morning for my own understanding, and the
> statements regarding st_read would appear to apply as explicitly to
> st_write as both or mediated by GDAL, that processes proj4string(s) but not
> epsg. So it looks like manually setting or resetting epsg is part of the
> workflow. Further down in the crs section is discussion about how within a
> layer both proj4string and epsg must be the same, or NA. This may localize
> where the 900914 is being applied, i.e. in PostGIS to settle the table
> parameters.
>
> HTH,
>
> Chris
>
>
>
>
>
> On Wed, Mar 15, 2017 at 3:50 PM, Michael Treglia <[hidden email]>
> wrote:
>
>> Hi All,
>>
>> Been working to import and manipulate a CSV file with point data
>> (lat/long), and then export to a PostGIS db.
>>
>> Overall, successful, but one thing I'd like to fix - when I write out the
>> layer to postgis, the SRID is not maintained. The final EPSG/SRID should
>> be
>> 2263, but when I check in PostGIS, it comes up as 900914.
>>
>> Below is my code and sessionInfo, and the data are from here:
>> https://data.cityofnewyork.us/Public-Safety/NYPD-Complaint-D
>> ata-Current-YTD/5uac-w243
>> (downloaded as plain old CSV)
>>
>> Anything I might be missing? Thanks in advance for giving a quick look!
>> Mike
>>
>>
>> ##Start Code
>>
>> #load packages
>> library(sf)
>> library(RPostgreSQL)
>>
>> #read data
>> crime_current <- read.csv("NYPD_Complaint_Data_Current_YTD.csv",
>> stringsAsFactors = FALSE)
>>
>> #format time columns for easier reading in postgres (I think), as keeping
>> as date format caused problems in export
>> crime_current$CMPLNT_FR_TIME <-
>> as.character(as.POSIXct(paste(crime_current$CMPLNT_FR_DT,
>> crime_current$CMPLNT_FR_TM), format="%m/%d/%Y\ %H:%M", tz=""))
>> crime_current$CMPLNT_TO_TIME <-
>> as.character(as.POSIXct(paste(crime_current$CMPLNT_TO_DT,
>> crime_current$CMPLNT_TO_TM), format="%m/%d/%Y\ %H:%M", tz=""))
>> crime_current$RPT_DT <- as.character(as.POSIXct(crime_current$RPT_DT,
>> format="%m/%d/%Y", tz=""))
>>
>> #convert to sf object
>> crime_current.sf <- st_as_sf(crime_current, coords = c("Longitude",
>> "Latitude"), crs = 4326)
>> #reproject to EPSG 2263
>> crime_current.sf <- st_transform(crime_current.sf, crs=2263)
>>
>> #write to postgres
>> st_write(crime_current.sf, "PG:dbname=mydb user=user host=xx.xx.xx.xx",
>> 'health_safety.crime_current')
>> ###End Code
>>
>>
>>
>> > sessionInfo()
>> R version 3.3.1 (2016-06-21)
>> Platform: x86_64-pc-linux-gnu (64-bit)
>> Running under: Ubuntu 14.04.5 LTS
>>
>> locale:
>>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
>> LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
>>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
>>  LC_PAPER=en_US.UTF-8       LC_NAME=C
>>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
>> LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>>
>> attached base packages:
>> [1] stats     graphics  grDevices utils     datasets  methods   base
>>
>> other attached packages:
>> [1] sp_1.2-3          RPostgreSQL_0.4-1 DBI_0.6           sf_0.3-4
>>
>> loaded via a namespace (and not attached):
>> [1] tools_3.3.1     units_0.4-2     Rcpp_0.12.9     udunits2_0.13
>> grid_3.3.1      lattice_0.20-33
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
        [[alternative HTML version deleted]]

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

Re: Maintain SRID with st_write to postgis db

Wed, 03/15/2017 - 15:21
Michael,

this from the sf vignettes on read write to databases:

The dsn and layer arguments to st_read and st_write denote a data source
name and optionally a layer name. Their exact interpretation as well as the
options they support vary per driver, the GDAL driver documentation is best
consulted for this. For instance, a PostGIS table in database postgis might
be read by

meuse <- st_read("PG:dbname=postgis", "meuse")

where the PG: string indicates this concerns the PostGIS driver, followed
by database name, and possibly port and user credentials.

st_read typically reads the coordinate reference system as proj4string, but
not the EPSG (SRID). GDAL cannot retrieve SRID (EPSG code) from proj4string
strings, and, when needed, it has to be set by the user.

I just re-read this this morning for my own understanding, and the
statements regarding st_read would appear to apply as explicitly to
st_write as both or mediated by GDAL, that processes proj4string(s) but not
epsg. So it looks like manually setting or resetting epsg is part of the
workflow. Further down in the crs section is discussion about how within a
layer both proj4string and epsg must be the same, or NA. This may localize
where the 900914 is being applied, i.e. in PostGIS to settle the table
parameters.

HTH,

Chris





On Wed, Mar 15, 2017 at 3:50 PM, Michael Treglia <[hidden email]> wrote:

> Hi All,
>
> Been working to import and manipulate a CSV file with point data
> (lat/long), and then export to a PostGIS db.
>
> Overall, successful, but one thing I'd like to fix - when I write out the
> layer to postgis, the SRID is not maintained. The final EPSG/SRID should be
> 2263, but when I check in PostGIS, it comes up as 900914.
>
> Below is my code and sessionInfo, and the data are from here:
> https://data.cityofnewyork.us/Public-Safety/NYPD-Complaint-
> Data-Current-YTD/5uac-w243
> (downloaded as plain old CSV)
>
> Anything I might be missing? Thanks in advance for giving a quick look!
> Mike
>
>
> ##Start Code
>
> #load packages
> library(sf)
> library(RPostgreSQL)
>
> #read data
> crime_current <- read.csv("NYPD_Complaint_Data_Current_YTD.csv",
> stringsAsFactors = FALSE)
>
> #format time columns for easier reading in postgres (I think), as keeping
> as date format caused problems in export
> crime_current$CMPLNT_FR_TIME <-
> as.character(as.POSIXct(paste(crime_current$CMPLNT_FR_DT,
> crime_current$CMPLNT_FR_TM), format="%m/%d/%Y\ %H:%M", tz=""))
> crime_current$CMPLNT_TO_TIME <-
> as.character(as.POSIXct(paste(crime_current$CMPLNT_TO_DT,
> crime_current$CMPLNT_TO_TM), format="%m/%d/%Y\ %H:%M", tz=""))
> crime_current$RPT_DT <- as.character(as.POSIXct(crime_current$RPT_DT,
> format="%m/%d/%Y", tz=""))
>
> #convert to sf object
> crime_current.sf <- st_as_sf(crime_current, coords = c("Longitude",
> "Latitude"), crs = 4326)
> #reproject to EPSG 2263
> crime_current.sf <- st_transform(crime_current.sf, crs=2263)
>
> #write to postgres
> st_write(crime_current.sf, "PG:dbname=mydb user=user host=xx.xx.xx.xx",
> 'health_safety.crime_current')
> ###End Code
>
>
>
> > sessionInfo()
> R version 3.3.1 (2016-06-21)
> Platform: x86_64-pc-linux-gnu (64-bit)
> Running under: Ubuntu 14.04.5 LTS
>
> locale:
>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
> LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
>  LC_PAPER=en_US.UTF-8       LC_NAME=C
>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
> LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] sp_1.2-3          RPostgreSQL_0.4-1 DBI_0.6           sf_0.3-4
>
> loaded via a namespace (and not attached):
> [1] tools_3.3.1     units_0.4-2     Rcpp_0.12.9     udunits2_0.13
> grid_3.3.1      lattice_0.20-33
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
        [[alternative HTML version deleted]]

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

Maintain SRID with st_write to postgis db

Wed, 03/15/2017 - 13:50
Hi All,

Been working to import and manipulate a CSV file with point data
(lat/long), and then export to a PostGIS db.

Overall, successful, but one thing I'd like to fix - when I write out the
layer to postgis, the SRID is not maintained. The final EPSG/SRID should be
2263, but when I check in PostGIS, it comes up as 900914.

Below is my code and sessionInfo, and the data are from here:
https://data.cityofnewyork.us/Public-Safety/NYPD-Complaint-Data-Current-YTD/5uac-w243
(downloaded as plain old CSV)

Anything I might be missing? Thanks in advance for giving a quick look!
Mike


##Start Code

#load packages
library(sf)
library(RPostgreSQL)

#read data
crime_current <- read.csv("NYPD_Complaint_Data_Current_YTD.csv",
stringsAsFactors = FALSE)

#format time columns for easier reading in postgres (I think), as keeping
as date format caused problems in export
crime_current$CMPLNT_FR_TIME <-
as.character(as.POSIXct(paste(crime_current$CMPLNT_FR_DT,
crime_current$CMPLNT_FR_TM), format="%m/%d/%Y\ %H:%M", tz=""))
crime_current$CMPLNT_TO_TIME <-
as.character(as.POSIXct(paste(crime_current$CMPLNT_TO_DT,
crime_current$CMPLNT_TO_TM), format="%m/%d/%Y\ %H:%M", tz=""))
crime_current$RPT_DT <- as.character(as.POSIXct(crime_current$RPT_DT,
format="%m/%d/%Y", tz=""))

#convert to sf object
crime_current.sf <- st_as_sf(crime_current, coords = c("Longitude",
"Latitude"), crs = 4326)
#reproject to EPSG 2263
crime_current.sf <- st_transform(crime_current.sf, crs=2263)

#write to postgres
st_write(crime_current.sf, "PG:dbname=mydb user=user host=xx.xx.xx.xx",
'health_safety.crime_current')
###End Code



> sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.5 LTS

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 LC_PAPER=en_US.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C
LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

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

other attached packages:
[1] sp_1.2-3          RPostgreSQL_0.4-1 DBI_0.6           sf_0.3-4

loaded via a namespace (and not attached):
[1] tools_3.3.1     units_0.4-2     Rcpp_0.12.9     udunits2_0.13
grid_3.3.1      lattice_0.20-33

        [[alternative HTML version deleted]]

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

Sentienl 2 gdal translate

Wed, 03/15/2017 - 12:49
Hi,

I have Sentinel 2 images (jp2) that I want to convert to GTiff using gdal_translate in R (from gdalutils package). Here is the code I’m using

Image_Path<-“/path/to/wd“

S2_JP2_List <- list.files(Image_Path, full.names = TRUE, pattern = ".jp2$")

for (i in 1:length(S2_JP2_List)) {
gdal_translate(src_dataset = S2_JP2_List[[i]], dst_dataset = paste("Band", i , "tif"), ot = "UInt16", of = “GTiff")
    }

When running this code, the process starts without any error but its never ending neither producing any output

I've tried to do it in gdal and the same code works, except that it is necessary to skip  the default driver since it cannot manage large jp2 files.

gdal_translate -of GTiff /path/to/input/S2_b1.jp2  /path/to/output/S2_b1converted.tif --config GDAL_SKIP JP2ECW

Any idea how could I run this in R?

Thanks
M



        [[alternative HTML version deleted]]

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

Getis Ord G test in 3 dimensions

Wed, 03/15/2017 - 11:43
> Dear all

> I am working on a 3D home range estimator, we want to calculated home ranges and core areas on primates that were followed sleeping site to sleeping site. We collected GPS locations and animals height at equal intervals along the day.
> We already calculates home ranges and core area in 2D using Delaunay triangulation and then applying a hot spot analysis in ArcGis( Getis Ord G*) using the perimeter of the triangles. Then, to calculate it on 3D we first created a 3D convex hull formed by tetrahedrons using animals locations (coordinates in meters)  and height.

> we are wondering if it is possible to use the function Local G or Global G test (package Spdep) for 3D locatios (X,Y,Z) as we previously  did in 2D (Hotspot with rendering, ArgGis)

Thanks in advance
>
>
> Juanma
                       
 
Juan Manuel José Domínguez
Postdoctoral Researcher        
                                                 
Conservation Ecology Program
School of Bioresources & Technology
King Mongkut´s University of Technology Thonburi

49 Soi Tienthalay 25
Bangkhuntien-Chaithalay Road
Thakham, Bangkhuntien
Bangkok 10150
Thailand
http://cons-ecol-kmutt.weebly.com/
        [[alternative HTML version deleted]]

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

Re: loop memory problem

Wed, 03/15/2017 - 11:37
Hi again,

I read in the guide of foreach that include the packages was a good idea.
Now the parallel loop is runing but at the end I got only one spatial line
(the first) I don't understand why is not adding the other spatiallines.

Any idea?

###3)dopar%
library(doSNOW)
registerDoSEQ()
registerDoParallel(cores=10)
getDoParWorkers()
line2<- gdistance::shortestPath(trCostS16, pos@coords[13,],pos@coords[12,],
output="SpatialLines")#

x<-foreach(i=2:13,.packages=c("gdistance","doParallel","foreach","base","sp"))
%dopar% {
  nextSegment2=gdistance::shortestPath(trCostS16,
pos@coords[i,],pos@coords[i+1,],
output="SpatialLines")
  line2 =nextSegment2+line2
  print(i)}

stopCluster(cl)
registerDoSEQ()

Thanks in advance,
Marta
####files
https://drive.google.com/open?id=0BwqSBe1Yq-FBUWVBOUdvaThEU1k

        [[alternative HTML version deleted]]

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

Re: Same mask for two rasters, but obtaining different extents in R

Wed, 03/15/2017 - 11:21
I re-checked the original extent and here it is: This is data obtained from
WorldClim.

> bio5
class       : RasterLayer
dimensions  : 3600, 8640, 31104000  (nrow, ncol, ncell)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -180, 180, -60, *90.00001 * (xmin, xmax, ymin, ymax) #This
varies at the 6th decimal place.
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
data source : C:\Users\rameshv\Downloads\Climate
Stability\Data_LGM_Present\Present\1_RawData\bio_5
names       : bio_5
values      : -86, 489  (min, max)
attributes  :
        ID COUNT
 from: -86     1
 to  : 489    38


> lg5
class       : RasterLayer
dimensions  : 3600, 8640, 31104000  (nrow, ncol, ncell)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
data source : C:\Users\rameshv\Downloads\Climate
Stability\Data_LGM_Present\LGM\1_RawData\cclgmbi5.tif
names       : cclgmbi5
values      : -204, 519  (min, max)

What do you suggest?

Also, any suggestions on aligning them before I crop them? I wouldn't want
to projectRaster using one of the rasters, as that creates extra number of
rows (when I coerce the raster to a Spatial Pixels Data Frame).



On Wed, Mar 15, 2017 at 1:13 PM, Sarah Goslee <[hidden email]>
wrote:

> The most likely explanation seems to me that the rasters were not
> aligned before cropping, so they are not aligned after cropping.
>
> Have you checked the original extent and resolution?
>
> Sarah
>
> On Tue, Mar 14, 2017 at 5:31 PM, Vijay Ramesh via R-sig-Geo
> <[hidden email]> wrote:
> > I loaded a shapefile in R, and cropped and masked it with the same mask,
> > but I get different extents. Any suggestions?
> >
> > Code below:
> >
> >     library(raster)
> >     library(rgdal)
> >     library(GISTools)
> >     library(sp)
> >     library(maptools)
> >
> >     ##Loading the first mask
> >     Mask <- readOGR("C:\\Users\\rameshv\\Downloads\\Climate
> > Stability\\SPV_Mask\\mask.shp")
> >     proj4string(Mask) <- CRS("+init=epsg:4326")
> >
> >     #Loading the Rasters For Present and LGM and clipping them to the
> right
> > extent
> >     bio5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate
> > Stability\\Data_LGM_Present\\Present\\1_RawData\\bio_5")
> >     proj4string(bio5) <- CRS("+init=epsg:4326")
> >     lg5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate
> > Stability\\Data_LGM_Present\\LGM\\1_RawData\\cclgmbi5.tif")
> >     proj4string(lg5) <- CRS("+init=epsg:4326")
> >
> >     ##Cropping by using the Crop and Mask functions
> >     cr<- crop(bio5, extent(Mask))
> >     bio5 <- mask(cr, Mask)
> >
> >     cr2<- crop(lg5, extent(Mask))
> >     lg5<- mask(cr2, Mask)
> >
> >     > bio5
> >      class       : RasterLayer
> >      dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
> >      resolution  : 0.04166667, 0.04166667  (x, y)
> >      extent      : 72.45835, 82.70835, 8.083337, 22.20834  (xmin, xmax,
> > ymin, ymax)
> >      coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
> > +ellps=WGS84 +towgs84=0,0,0
> >      data source : in memory
> >      names       : bio_5
> >      values      : 207, 432  (min, max)
> >      attributes  :
> >         ID COUNT
> >      from: -86     1
> >      to  : 489    38
> >
> >      > lg5
> >      class       : RasterLayer
> >      dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
> >      resolution  : 0.04166667, 0.04166667  (x, y)
> >      extent      : 72.45833, 82.70833, 8.083333, 22.20833  (xmin, xmax,
> > ymin, ymax)
> >      coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
> > +ellps=WGS84 +towgs84=0,0,0
> >      data source : in memory
> >      names       : cclgmbi5
> >      values      : 187, 392  (min, max)
> >
> >
>
> --
> Sarah Goslee
> http://www.functionaldiversity.org
>


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

        [[alternative HTML version deleted]]

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

Re: Same mask for two rasters, but obtaining different extents in R

Wed, 03/15/2017 - 11:13
The most likely explanation seems to me that the rasters were not
aligned before cropping, so they are not aligned after cropping.

Have you checked the original extent and resolution?

Sarah

On Tue, Mar 14, 2017 at 5:31 PM, Vijay Ramesh via R-sig-Geo
<[hidden email]> wrote:
> I loaded a shapefile in R, and cropped and masked it with the same mask,
> but I get different extents. Any suggestions?
>
> Code below:
>
>     library(raster)
>     library(rgdal)
>     library(GISTools)
>     library(sp)
>     library(maptools)
>
>     ##Loading the first mask
>     Mask <- readOGR("C:\\Users\\rameshv\\Downloads\\Climate
> Stability\\SPV_Mask\\mask.shp")
>     proj4string(Mask) <- CRS("+init=epsg:4326")
>
>     #Loading the Rasters For Present and LGM and clipping them to the right
> extent
>     bio5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate
> Stability\\Data_LGM_Present\\Present\\1_RawData\\bio_5")
>     proj4string(bio5) <- CRS("+init=epsg:4326")
>     lg5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate
> Stability\\Data_LGM_Present\\LGM\\1_RawData\\cclgmbi5.tif")
>     proj4string(lg5) <- CRS("+init=epsg:4326")
>
>     ##Cropping by using the Crop and Mask functions
>     cr<- crop(bio5, extent(Mask))
>     bio5 <- mask(cr, Mask)
>
>     cr2<- crop(lg5, extent(Mask))
>     lg5<- mask(cr2, Mask)
>
>     > bio5
>      class       : RasterLayer
>      dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
>      resolution  : 0.04166667, 0.04166667  (x, y)
>      extent      : 72.45835, 82.70835, 8.083337, 22.20834  (xmin, xmax,
> ymin, ymax)
>      coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
> +ellps=WGS84 +towgs84=0,0,0
>      data source : in memory
>      names       : bio_5
>      values      : 207, 432  (min, max)
>      attributes  :
>         ID COUNT
>      from: -86     1
>      to  : 489    38
>
>      > lg5
>      class       : RasterLayer
>      dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
>      resolution  : 0.04166667, 0.04166667  (x, y)
>      extent      : 72.45833, 82.70833, 8.083333, 22.20833  (xmin, xmax,
> ymin, ymax)
>      coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
> +ellps=WGS84 +towgs84=0,0,0
>      data source : in memory
>      names       : cclgmbi5
>      values      : 187, 392  (min, max)
>
>
--
Sarah Goslee
http://www.functionaldiversity.org

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

Re: Same mask for two rasters, but obtaining different extents in R

Wed, 03/15/2017 - 07:58
Hi,

Might you show us what bio5 and lg5 look like before you do the masking?

Cheers,
Ben


> On Mar 14, 2017, at 5:31 PM, Vijay Ramesh via R-sig-Geo <[hidden email]> wrote:
>
> I loaded a shapefile in R, and cropped and masked it with the same mask,
> but I get different extents. Any suggestions?
>
> Code below:
>
>    library(raster)
>    library(rgdal)
>    library(GISTools)
>    library(sp)
>    library(maptools)
>
>    ##Loading the first mask
>    Mask <- readOGR("C:\\Users\\rameshv\\Downloads\\Climate
> Stability\\SPV_Mask\\mask.shp")
>    proj4string(Mask) <- CRS("+init=epsg:4326")
>
>    #Loading the Rasters For Present and LGM and clipping them to the right
> extent
>    bio5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate
> Stability\\Data_LGM_Present\\Present\\1_RawData\\bio_5")
>    proj4string(bio5) <- CRS("+init=epsg:4326")
>    lg5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate
> Stability\\Data_LGM_Present\\LGM\\1_RawData\\cclgmbi5.tif")
>    proj4string(lg5) <- CRS("+init=epsg:4326")
>
>    ##Cropping by using the Crop and Mask functions
>    cr<- crop(bio5, extent(Mask))
>    bio5 <- mask(cr, Mask)
>
>    cr2<- crop(lg5, extent(Mask))
>    lg5<- mask(cr2, Mask)
>
>> bio5
>     class       : RasterLayer
>     dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
>     resolution  : 0.04166667, 0.04166667  (x, y)
>     extent      : 72.45835, 82.70835, 8.083337, 22.20834  (xmin, xmax,
> ymin, ymax)
>     coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
> +ellps=WGS84 +towgs84=0,0,0
>     data source : in memory
>     names       : bio_5
>     values      : 207, 432  (min, max)
>     attributes  :
>        ID COUNT
>     from: -86     1
>     to  : 489    38
>
>> lg5
>     class       : RasterLayer
>     dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
>     resolution  : 0.04166667, 0.04166667  (x, y)
>     extent      : 72.45833, 82.70833, 8.083333, 22.20833  (xmin, xmax,
> ymin, ymax)
>     coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
> +ellps=WGS84 +towgs84=0,0,0
>     data source : in memory
>     names       : cclgmbi5
>     values      : 187, 392  (min, max)
>
>
> --
> Data Manager,
> Barbara Han Lab,
> Cary Institute of Ecosystem Studies,
> 2801 Sharon Turnpike, Millbrook, NY 12545
> Lab Site : http://www.hanlab.science/
> Personal Site : http://evolecol.weebly.com/
> Phone : (845)-677-7600 Ext: 241
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org

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

Same mask for two rasters, but obtaining different extents in R

Tue, 03/14/2017 - 15:31
I loaded a shapefile in R, and cropped and masked it with the same mask,
but I get different extents. Any suggestions?

Code below:

    library(raster)
    library(rgdal)
    library(GISTools)
    library(sp)
    library(maptools)

    ##Loading the first mask
    Mask <- readOGR("C:\\Users\\rameshv\\Downloads\\Climate
Stability\\SPV_Mask\\mask.shp")
    proj4string(Mask) <- CRS("+init=epsg:4326")

    #Loading the Rasters For Present and LGM and clipping them to the right
extent
    bio5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate
Stability\\Data_LGM_Present\\Present\\1_RawData\\bio_5")
    proj4string(bio5) <- CRS("+init=epsg:4326")
    lg5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate
Stability\\Data_LGM_Present\\LGM\\1_RawData\\cclgmbi5.tif")
    proj4string(lg5) <- CRS("+init=epsg:4326")

    ##Cropping by using the Crop and Mask functions
    cr<- crop(bio5, extent(Mask))
    bio5 <- mask(cr, Mask)

    cr2<- crop(lg5, extent(Mask))
    lg5<- mask(cr2, Mask)

    > bio5
     class       : RasterLayer
     dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
     resolution  : 0.04166667, 0.04166667  (x, y)
     extent      : 72.45835, 82.70835, 8.083337, 22.20834  (xmin, xmax,
ymin, ymax)
     coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
+ellps=WGS84 +towgs84=0,0,0
     data source : in memory
     names       : bio_5
     values      : 207, 432  (min, max)
     attributes  :
        ID COUNT
     from: -86     1
     to  : 489    38

     > lg5
     class       : RasterLayer
     dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
     resolution  : 0.04166667, 0.04166667  (x, y)
     extent      : 72.45833, 82.70833, 8.083333, 22.20833  (xmin, xmax,
ymin, ymax)
     coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
+ellps=WGS84 +towgs84=0,0,0
     data source : in memory
     names       : cclgmbi5
     values      : 187, 392  (min, max)


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

        [[alternative HTML version deleted]]

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

Re: loop memory problem

Tue, 03/14/2017 - 11:02
> #answer Marcel in r-sig-geo forum:
>> #1) maybe you can first make your calculations and afterwards plot the
>> result.
>> #It might be faster this way.
>>
>> #2) maybe you can use the parallel for each loop in order to use the
>> whole performance of your cpu
>>
>> #3) loops in R are really slow in general.. you could also think about
>> some fancy stuff like some R compiler or Rcpp package
>>
>> #4) in order to analyze your memory consumption, you should have a look
>> on your system resources (depending on your OS: task manager (win) or htop
>> (linux))
>>
>> #############
>> #1)MARCEL: maybe you can first make your calculations and afterwards plot
>> the result.
>> #It might be faster this way.
>> ##MARTA: I did that, but the difference between the calculations inside
>> the loop of afterwards is less than a second.
>> #
>>
>       #############
> #2) MARCEL:maybe you can use the parallel or each loop in order to use the
>> whole performance of your cpu
>> #MARTA:I  developed 5 different loops, the old one (2.A) it works
>> properly but is too slow. I followed your suggestion with the parallel
>> loops to increase the speed. I wrote a loop (2.B) with %do%, which is works
>> but the time is the #same as the old loop. The other loops with %doPar%,
>> didn't work at all. Only the simple (2.D) runs, but is not what I want. The
>> other two ( 2.C and 2.E) didn't run, they have errors. The 2.C "task 1
>> failed - "non-numeric argument #to binary operator"  and the 2.E "task 1
>> failed - "no method for coercing this S4 class to a vector".
>> #
>> #Any idea of how repair this loop?
>>
>> #########
>> #library
>> ########
>> library(parallel);library(foreach);library(doParallel);
>> library(gdistance);library(raster)
>> library(rgdal);library(rgeos)
>> library(sp)
>> #######
>> #data
>> ###########
>> path_data<-"E:/Q11/"
>> # cost raster sea and island
>> costa6Azo <- raster("E:/Q11/costa6Azo_projected.tif")
>> transitioncosta6Azo <- transition(costa6Azo, min, directions=16)#porque
>> min????
>> trCostS16 <- gdistance::geoCorrection(transitioncosta6Azo, type="c")
>> #points
>> boat <- read.table(paste0(path_data,"boat2905.csv"), header=TRUE,
>> sep=";", na.strings="NA", dec=".", strip.white=TRUE)
>> pos<-as.data.frame(cbind(boat$Lat1,boat$Long1,boat$Ref))
>> str(pos)
>> sp::coordinates(pos) <- ~V2+V1
>> sp::proj4string(pos) <-CRS("+proj=longlat +datum=WGS84 +no_defs")
>>
>> pos<-sp::spTransform(pos, CRS("+proj=utm +zone=26 +ellps=intl
>> +towgs84=-104,167,-38,0,0,0,0 +units=m +no_defs"))
>> # The aim of the loop: calculating the track of the whole sailing
>> path#############################################################
>> #
>> line<- gdistance::shortestPath(trCostS16, pos@coords[1,],pos@coords[2,],
>> output="SpatialLines")#
>> glength(1-2=4887.737m);(2-3=12590.01);(12-11=11360.39m);(12-13=9453.001m)
>> lines(line,col=5)
>> #2. A) old
>> loop##########################################################################################################################
>> #it works, but it's too slow
>>
>
> ## here we start with the for-loop
>> for (i in (seq(2,length(pos) - 1))) {
>>   # calculation of the rest of the segements
>>   nextSegment<- gdistance::shortestPath(trCostS16, pos@coords
>> [i,],pos@coords[i+1,], output="SpatialLines")
>>   # simple addition combines the single spatialline segements
>>   line <- nextSegment + line
>>   # we plot each new segment
>>   lines(nextSegment)
>> }
>> # note that we have now ten combined line features in this SpatialLines
>> object
>> line
>> gLength(line)#110747.2
>>
>> #2.B### new parallel
>> loop##################################################################################################3
>> #it works, but it's too slow##%do%
>> #
>> line<- gdistance::shortestPath(trCostS16, pos@coords[1,],pos@coords[2,],
>> output="SpatialLines")#
>> glength(1-2=4887.737m);(2-3=12590.01);(12-11=11360.39m);(12-13=9453.001m)
>>
>> x <- foreach(i=2:13) %do%
>>   {nextSegment<-gdistance::shortestPath(trCostS16, pos@coords
>> [i,],pos@coords[i-1,], output="SpatialLines")
>>     line <- nextSegment + line
>>  }
>> x      ##  it works!!!! 1+ 12 spatialLines!!
>> ##
>
> #     #2.C#parallel loop ##%dopar%

> #without success
>>
> registerDoParallel()
>
> getDoParWorkers()
>> line<- gdistance::shortestPath(trCostS16, pos@coords[1,],pos@coords[2,],
>> output="SpatialLines")#
>> glength(1-2=4887.737m);(2-3=12590.01);(12-11=11360.39m);(12-13=9453.001m)
>> #function
>> #
>> funMTM<-function(){
>> nextSegment<-gdistance::shortestPath(trCostS16, pos@coords[i,],pos@coords[i-1,],
>> output="SpatialLines")
>> line <- nextSegment + line
>> }
>> getDoParWorkers()
>> ptime <- system.time({
>>   result <- foreach(i=2:13) %dopar% funMTM()
>>   })
>> ptime#Error in funMTM() :
>> #task 1 failed - "non-numeric argument to binary operator"
>> #
>>
> #
> #2.D#loop %dopar% simple
>>
>        # it works, but I need the spatiallines output, not a list.

> #If I run the %dopar% only with the functions shortestPath without
>> increase the lines into the SpatialLines. I get a list with 13 features(
>> indidivual spatialLines). However, I need an SpatialLines object, with 13
>> spatiallines inside.
>> registerDoParallel()
>> registerDoSEQ()
>> registerDoParallel(cores=10)
>> getDoParWorkers()
>> system.time(foreach(i=2:13) %dopar% gdistance::shortestPath(trCostS16,
>> pos@coords[i,],pos@coords[i-1,], output="SpatialLines"))
>> #user  system elapsed
>> #0.74    0.69   54.81
>>
> #    #2.E#loop %dopar% with the function defined in the loop
      #without success:

> system.time(foreach(i=2:13) %dopar% {
>>   nextSegment=gdistance::shortestPath(trCostS16, pos@coords
>> [i,],pos@coords[i-1,], output="SpatialLines")
>>             line =nextSegment + line})
>> #Error in { :
>> #task 1 failed - "no method for coercing this S4 class to a vector"
>> stopCluster(cl)
>>
>> ##############
>> #3)MARCEL: loops in R are really slow in general.. you could also think
>> about some fancy stuff like some R compiler or Rcpp package
>> #MARTA##### Functions #####
>> # byte code compilation
>> library(compiler)
>> myfunc<-gdistance::shortestPath(trCostS16, pos@coords[1,],pos@coords[2,],
>> output="SpatialLines")
>> myFuncCmp <- cmpfun(myfunc)
>> system.time({
>>   output <- SpatialLines(LinesList = , 1, FUN=myFuncCmp)
>> })
>> #############
>>
>>       #############
> #4) MARCEL: in order to analyze your memory consumption, you should have a
>> look on your system resources (depending on your OS: task manager (win) or
>> htop (linux))
>> #MARTA:I run the loop with the task manager's window open, and never over
>> pass the 30% of the CPU's memory .
>> ####files
>> https://drive.google.com/open?id=0BwqSBe1Yq-FBUWVBOUdvaThEU1k
>>
>>
> I haven't my solution yet but I'm closer now, your suggestions were really helpful.
Marta

        [[alternative HTML version deleted]]

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

Re: Error in bbox(newdata) : object not a >= 2-column array

Tue, 03/14/2017 - 03:00
Looks to me as if you made a typo in your script, possibly with a quote,
but I'd need to see the full script, and preferably have access the
necessary data in order to run it, to be completely sure.

On 14/03/17 09:54, [hidden email] wrote:
>> Dear David,
> Dear R users,
> Now R calculates var1.pred as a constatnt:
> kriged=krige(Nmin~1,locationsDD,grid,model=fit)
>
>> kriged["var1.pred"]
>
>
> 925  (339003.4, 2988160)  11.22117
> 926  (339013.4, 2988160)  11.22117
> 927  (339023.4, 2988160)  11.22117
> 928  (339033.4, 2988160)  11.22117
> 929  (339043.4, 2988160)  11.22117
> 930  (339053.4, 2988160)  11.22117
> 931  (339063.4, 2988160)  11.22117
> 932  (339073.4, 2988160)  11.22117
> 933  (339083.4, 2988160)  11.22117
> ........................................
> 7482 (339293.4, 2989180)  11.22117
> 7483 (339303.4, 2989180)  11.22117
> 7484 (339313.4, 2989180)  11.22117
> 7485 (339323.4, 2989180)  11.22117
> 7486 (339333.4, 2989180)  11.22117
> 7487 (339343.4, 2989180)  11.22117
> 7488 (339353.4, 2989180)  11.22117
>>
> spplot(kriged["var1.pred"])
>  Error: unexpected string constant in:
> "
> spplot(kriged[""
>
> And it cann't give ssplot
> Please if anyone can help me!
>
>
> I was getting this error as well.  The cause was either not having the
>> coordinates and coordinate reference system (CRS) set or having the wrong
>> CRS.  I did both at the same time so I am not sure which it was.
>>
>> Is your 'grid' a SpatialPixels object?
>>
>> Did you run a line like 'coordinates(grid) <- ~x + y'  ?
>>
>> If so, did you then set the CRS?
>>
>>
>>
>> On Sat, Mar 11, 2017 at 7:38 AM, <[hidden email]> wrote:
>>
>>> Dear R-Users,
>>> Please anyone to help me with this error:
>>>
>>> kriged=krige(Nmin~1,locationsDD,grid,model=vgm(1,"Exp",150,1))
>>> Error in bbox(newdata) : object not a >= 2-column array
>>>
>>> Thank you so much in advance!
>>> Kind Regards
>>> Nina Philipova
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>>
>>
>> --
>> David Warner
>> Research Fisheries Biologist
>> U.S.G.S. Great Lakes Science Center
>> 1451 Green Road
>> Ann Arbor, MI 48105
>> 734-214-9392
>>
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo--
Edzer Pebesma
Institute for Geoinformatics  (ifgi),  University of Münster
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers & Geosciences:   http://elsevier.com/locate/cageo/


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

Re: Error in bbox(newdata) : object not a >= 2-column array

Tue, 03/14/2017 - 02:54
> Dear David,
Dear R users,
Now R calculates var1.pred as a constatnt:
kriged=krige(Nmin~1,locationsDD,grid,model=fit)

> kriged["var1.pred"]


925  (339003.4, 2988160)  11.22117
926  (339013.4, 2988160)  11.22117
927  (339023.4, 2988160)  11.22117
928  (339033.4, 2988160)  11.22117
929  (339043.4, 2988160)  11.22117
930  (339053.4, 2988160)  11.22117
931  (339063.4, 2988160)  11.22117
932  (339073.4, 2988160)  11.22117
933  (339083.4, 2988160)  11.22117
........................................
7482 (339293.4, 2989180)  11.22117
7483 (339303.4, 2989180)  11.22117
7484 (339313.4, 2989180)  11.22117
7485 (339323.4, 2989180)  11.22117
7486 (339333.4, 2989180)  11.22117
7487 (339343.4, 2989180)  11.22117
7488 (339353.4, 2989180)  11.22117
>
spplot(kriged["var1.pred"])
 Error: unexpected string constant in:
"
spplot(kriged[""

And it cann't give ssplot
Please if anyone can help me!


I was getting this error as well.  The cause was either not having the
> coordinates and coordinate reference system (CRS) set or having the wrong
> CRS.  I did both at the same time so I am not sure which it was.
>
> Is your 'grid' a SpatialPixels object?
>
> Did you run a line like 'coordinates(grid) <- ~x + y'  ?
>
> If so, did you then set the CRS?
>
>
>
> On Sat, Mar 11, 2017 at 7:38 AM, <[hidden email]> wrote:
>
>> Dear R-Users,
>> Please anyone to help me with this error:
>>
>> kriged=krige(Nmin~1,locationsDD,grid,model=vgm(1,"Exp",150,1))
>> Error in bbox(newdata) : object not a >= 2-column array
>>
>> Thank you so much in advance!
>> Kind Regards
>> Nina Philipova
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
>
> --
> David Warner
> Research Fisheries Biologist
> U.S.G.S. Great Lakes Science Center
> 1451 Green Road
> Ann Arbor, MI 48105
> 734-214-9392
>
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Variogram & Kriging.txt (3K) Download Attachment

Re: Avoid axis labels, ticks and marks in ggRGB() plot

Mon, 03/13/2017 - 07:25
Yes! Thanks!
I knew there had to be a better way. For my set of plots I can set
t <- theme(
  axis.text = element_blank(),
  axis.ticks = element_blank(),
  axis.title = element_blank())

and then

ggRGB(rlogo, r=1, g=2, b=3) + t +ggtitle("patatin")

Agus



On Mon, Mar 13, 2017 at 2:05 PM, Benjamin Leutner
<[hidden email]> wrote:
> The problem with your first approach is that there is no "data layer" which
> could span the coordinate axes.
> Since ggRGB by default uses an annotation_raster() in order to keep your
> "fill" dimension free for other data,  the coordinate limits are empty,
> hence the plotted raster is lost in the dimensionless void and you get see
> nothing ;-)
>
> I would recommend to start your plot with ggRGB like this:
>
> ggRGB(rlogo, r=1, g=2, b=3) +
>   theme(
>     axis.text = element_blank(),
>     axis.ticks = element_blank(),
>     axis.title = element_blank())
>
>
> If you do want to use it with ggLayer = TRUE, there should be at least one
> "geom" with actual data in it (or you'd have to set ggRGB(...,
> geom_raster=TRUE), but I don't see any benefit in that and it would give you
> a huge legend with one entry per rgb color...)
>
> best,
> Benjamin
>
>
> On 03/13/2017 12:10 PM, Agustin Lobo wrote:
>>
>> In order to save space, I prefer to avoid coordinates in the plot
>> generated by ggRGB. Is there any way? I've tried:
>>
>> ggplot() +
>> ggRGB(rlogo, r=1, g=2, b=3,ggLayer=TRUE) +
>>    theme(
>>      axis.text = element_blank(),
>>      axis.ticks = element_blank(),
>>      axis.title = element_blank())
>>
>> but get an empty plot.
>> I can circumvent the problem as follows, but it is kind of ugly and
>> complicated,
>> there must be another way:
>>
>> a <- as(extent(rlogo), 'SpatialPolygons')
>> y <- fortify(a,region=id)
>> ggplot(data=y,aes(y=lat, x=long)) +
>>    ggRGB(rlogo, r=1, g=2, b=3,ggLayer=TRUE) +
>>    geom_polygon(colour="yellow", fill=NA) +
>>    guides(fill=FALSE) +
>>    coord_fixed(ratio=1) +
>>    theme(
>>      axis.text = element_blank(),
>>      axis.ticks = element_blank(),
>>      axis.title = element_blank())
>>
>> Agus
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
> Benjamin Leutner M.Sc.
> PhD candidate
>
> Department of Remote Sensing
> University of Wuerzburg
> Campus Hubland Nord 86
> 97074 Wuerzburg, Germany
>
> Tel: +49-(0)931-31 89594
> Fax: +49-(0)931-31 89594-0
> Email: [hidden email]
> Web: http://www.fernerkundung.uni-wuerzburg.de
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Avoid axis labels, ticks and marks in ggRGB() plot

Mon, 03/13/2017 - 07:05
The problem with your first approach is that there is no "data layer"
which could span the coordinate axes.
Since ggRGB by default uses an annotation_raster() in order to keep your
"fill" dimension free for other data,  the coordinate limits are empty,
hence the plotted raster is lost in the dimensionless void and you get
see nothing ;-)

I would recommend to start your plot with ggRGB like this:

ggRGB(rlogo, r=1, g=2, b=3) +
   theme(
     axis.text = element_blank(),
     axis.ticks = element_blank(),
     axis.title = element_blank())


If you do want to use it with ggLayer = TRUE, there should be at least
one "geom" with actual data in it (or you'd have to set ggRGB(...,
geom_raster=TRUE), but I don't see any benefit in that and it would give
you a huge legend with one entry per rgb color...)

best,
Benjamin

On 03/13/2017 12:10 PM, Agustin Lobo wrote:
> In order to save space, I prefer to avoid coordinates in the plot
> generated by ggRGB. Is there any way? I've tried:
>
> ggplot() +
> ggRGB(rlogo, r=1, g=2, b=3,ggLayer=TRUE) +
>    theme(
>      axis.text = element_blank(),
>      axis.ticks = element_blank(),
>      axis.title = element_blank())
>
> but get an empty plot.
> I can circumvent the problem as follows, but it is kind of ugly and complicated,
> there must be another way:
>
> a <- as(extent(rlogo), 'SpatialPolygons')
> y <- fortify(a,region=id)
> ggplot(data=y,aes(y=lat, x=long)) +
>    ggRGB(rlogo, r=1, g=2, b=3,ggLayer=TRUE) +
>    geom_polygon(colour="yellow", fill=NA) +
>    guides(fill=FALSE) +
>    coord_fixed(ratio=1) +
>    theme(
>      axis.text = element_blank(),
>      axis.ticks = element_blank(),
>      axis.title = element_blank())
>
> Agus
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Benjamin Leutner M.Sc.
PhD candidate

Department of Remote Sensing
University of Wuerzburg
Campus Hubland Nord 86
97074 Wuerzburg, Germany

Tel: +49-(0)931-31 89594
Fax: +49-(0)931-31 89594-0
Email: [hidden email]
Web: http://www.fernerkundung.uni-wuerzburg.de

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

Avoid axis labels, ticks and marks in ggRGB() plot

Mon, 03/13/2017 - 05:10
In order to save space, I prefer to avoid coordinates in the plot
generated by ggRGB. Is there any way? I've tried:

ggplot() +
ggRGB(rlogo, r=1, g=2, b=3,ggLayer=TRUE) +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank())

but get an empty plot.
I can circumvent the problem as follows, but it is kind of ugly and complicated,
there must be another way:

a <- as(extent(rlogo), 'SpatialPolygons')
y <- fortify(a,region=id)
ggplot(data=y,aes(y=lat, x=long)) +
  ggRGB(rlogo, r=1, g=2, b=3,ggLayer=TRUE) +
  geom_polygon(colour="yellow", fill=NA) +
  guides(fill=FALSE) +
  coord_fixed(ratio=1) +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank())

Agus

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

Does ggmap work with sf objects through geom_sf?

Sun, 03/12/2017 - 07:29
Dear list members,

Does ggmap work with sf objects through geom_sf?

I got the following error:

nc <- st_read(system.file("shape/nc.shp", package="sf"))

nclocation <- c(-80, 34)

ncmap <- get_map(location = nclocation, zoom = 6)

ggmap(ncmap) + geom_sf(data = nc)

Error in eval(expr, envir, enclos) : object 'lon' not found


--
*Manuel Spínola, Ph.D.*
Instituto Internacional en Conservación y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
[hidden email] <[hidden email]>
[hidden email]
Teléfono: (506) 8706 - 4662
Personal website: Lobito de río <https://sites.google.com/site/lobitoderio/>
Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/>

        [[alternative HTML version deleted]]

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

Re: Error in bbox(newdata) : object not a >= 2-column array

Sat, 03/11/2017 - 09:13
I was getting this error as well.  The cause was either not having the
coordinates and coordinate reference system (CRS) set or having the wrong
CRS.  I did both at the same time so I am not sure which it was.

Is your 'grid' a SpatialPixels object?

Did you run a line like 'coordinates(grid) <- ~x + y'  ?

If so, did you then set the CRS?



On Sat, Mar 11, 2017 at 7:38 AM, <[hidden email]> wrote:

> Dear R-Users,
> Please anyone to help me with this error:
>
> kriged=krige(Nmin~1,locationsDD,grid,model=vgm(1,"Exp",150,1))
> Error in bbox(newdata) : object not a >= 2-column array
>
> Thank you so much in advance!
> Kind Regards
> Nina Philipova
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>


--
David Warner
Research Fisheries Biologist
U.S.G.S. Great Lakes Science Center
1451 Green Road
Ann Arbor, MI 48105
734-214-9392

        [[alternative HTML version deleted]]

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

Re: Error in bbox(newdata) : object not a >= 2-column array

Sat, 03/11/2017 - 08:03
On Sat, Mar 11, 2017 at 12:38 PM,  <[hidden email]> wrote:
> Dear R-Users,
> Please anyone to help me with this error:
>
> kriged=krige(Nmin~1,locationsDD,grid,model=vgm(1,"Exp",150,1))
> Error in bbox(newdata) : object not a >= 2-column array
>
> Thank you so much in advance!

 Hard to debug without your data. The error message is telling us that
something is wrong with your input, but we can't see your data....

 If you can't share your data, at least share what you get from
summary(grid) and summary(locationsDD).

 You should also confirm what package `krige` is coming from - I
assume its gstat.

If you can get back to the mailing list with this information, we
might be able to help, otherwise we're guessing!

Barry

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

Error in bbox(newdata) : object not a >= 2-column array

Sat, 03/11/2017 - 06:38
Dear R-Users,
Please anyone to help me with this error:

kriged=krige(Nmin~1,locationsDD,grid,model=vgm(1,"Exp",150,1))
Error in bbox(newdata) : object not a >= 2-column array

Thank you so much in advance!
Kind Regards
Nina Philipova
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
variogram2.txt (2K) Download Attachment

Pages