Re: [R-sig-Geo] Merging Shape files

[Zia Ahmed]

You welcome.

If you are interested in geospatial analysis in python, you check this:

https://zia207.github.io/geospatial-python.io/

Best
Zia

On Wed, Mar 18, 2020, 10:45 AM Dexter Locke  wrote:

> Zia: this is a fabulous resource!!
>
> Thanks for sharing,
> Dexter
>
>
> On Wed, Mar 18, 2020 at 10:37 AM Zia Ahmed  wrote:
>
>> Please  check this:
>>
>>
>> https://zia207.github.io/geospatial-r-github.io/geoprocessing-vector-data.html
>> Best
>> Zia
>>
>> On Tue, Mar 17, 2020, 10:39 AM Dexter Locke 
>> wrote:
>>
>>> Within st_join() see the options like "left = TRUE".
>>>
>>> You may also want to look at st_intersects()
>>>
>>> Here are other helpful reference websites:
>>> https://r-spatial.github.io/sf/reference/geos_binary_pred.html
>>>
>>> https://cran.r-project.org/web/packages/sf/vignettes/sf3.html
>>>
>>> -Dexter
>>>
>>>
>>> On Tue, Mar 17, 2020 at 10:29 AM El Mechry El Koudouss <
>>> eelkoudo...@fordham.edu> wrote:
>>>
>>> > Dear readers,
>>> > I downloaded a shapefile with data on crime incidents in Washington DC
>>> in
>>> > 2016 (available here:
>>> > https://opendata.dc.gov/datasets/crime-incidents-in-2016) and another
>>> > shapefile with data on parks and recreation areas, also in Washington
>>> DC
>>> > (available here:
>>> > https://opendata.dc.gov/datasets/parks-and-recreation-areas).
>>> > I was able to read both shapefiles using st_read() from the sf
>>> package. My
>>> > question is, how can I merge the two data sets? I tried st_join() but
>>> it
>>> > simply resulted in empty columns from the parks dataset getting added
>>> to
>>> > the crime data. Any guidance would be much appreciated.
>>> > library(sf)
>>> > parks <- st_read("Parks_and_Recreation_Areas.shp")
>>> > crime <- st_read("Crime_Incidents_in_2016.shp")
>>> > df <- st_join(crime, parks)
>>> >
>>> > --
>>> > El Mechry, El Koudouss (Meshry)
>>> > Graduate Research Assistant
>>> > Center for International Policy Studies
>>> > Fordham University
>>> > Department of Economics
>>> > Website: www.meshry.com
>>> >
>>> > [[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
>>>
>>

[[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] Merging Shape files

[Zia Ahmed]

Please  check this:

https://zia207.github.io/geospatial-r-github.io/geoprocessing-vector-data.html
Best
Zia

On Tue, Mar 17, 2020, 10:39 AM Dexter Locke  wrote:

> Within st_join() see the options like "left = TRUE".
>
> You may also want to look at st_intersects()
>
> Here are other helpful reference websites:
> https://r-spatial.github.io/sf/reference/geos_binary_pred.html
>
> https://cran.r-project.org/web/packages/sf/vignettes/sf3.html
>
> -Dexter
>
>
> On Tue, Mar 17, 2020 at 10:29 AM El Mechry El Koudouss <
> eelkoudo...@fordham.edu> wrote:
>
> > Dear readers,
> > I downloaded a shapefile with data on crime incidents in Washington DC in
> > 2016 (available here:
> > https://opendata.dc.gov/datasets/crime-incidents-in-2016) and another
> > shapefile with data on parks and recreation areas, also in Washington DC
> > (available here:
> > https://opendata.dc.gov/datasets/parks-and-recreation-areas).
> > I was able to read both shapefiles using st_read() from the sf package.
> My
> > question is, how can I merge the two data sets? I tried st_join() but it
> > simply resulted in empty columns from the parks dataset getting added to
> > the crime data. Any guidance would be much appreciated.
> > library(sf)
> > parks <- st_read("Parks_and_Recreation_Areas.shp")
> > crime <- st_read("Crime_Incidents_in_2016.shp")
> > df <- st_join(crime, parks)
> >
> > --
> > El Mechry, El Koudouss (Meshry)
> > Graduate Research Assistant
> > Center for International Policy Studies
> > Fordham University
> > Department of Economics
> > Website: www.meshry.com
> >
> > [[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
>

[[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] Extract VIIRS data in R

[Zia Ahmed]

This is might help you:

https://zia207.github.io/geospatial-r-github.io/netCDF-data-processing.html

Zia

On Wed, Mar 18, 2020, 9:34 AM Frederico Faleiro  wrote:

> Hi Miluji,
>
> the R in general do not handle very well with netCDF. However I
> recommend you use the netcdf4 package before try others.
> This link is a good start:
>
> http://geog.uoregon.edu/bartlein/courses/geog490/week04-netCDF.html#get-coordinate-including-time-variables
> .
> The coordinates of this file is not in longitude and latitude, it use
> "number_of_lines" for X Axis and "number_of_pixels" for Y Axis. This way
> the difference between the cells is always one and raster package
> interprete the difference as one degree and return an warning about it.
> I open the file in the software Panoply (design to work with netCDF)  and
> everything is OK apparently.
> Check bellow some exploration of the data in R.
>
> library(raster)
> library(ncdf4)
> library(rgdal)
>
> # netcdf4 package
> netcdf_file <- "VNP02DNB_NRT.A2020069.1048.001.nc"
> nc <- nc_open(netcdf_file)
> print(nc)
> # check your variable is structured as:
> observation_data/DNB_observations[number_of_pixels,number_of_lines]
> # variables available
> names(nc[["var"]])
> # properties from your variable
> ncatt_get(nc, "observation_data/DNB_observations")
> # get the variables values and check the structure
> obs_matrix <- ncvar_get(nc, "observation_data/DNB_observations")
> str(obs_matrix)
>
> # raster package
> nc_r <- brick(netcdf_file, lvar=0, values=TRUE,
> varname="observation_data/DNB_observations")
> # [1] "vobjtovarid4:  WARNING  I was asked to get a varid for
> dimension named number_of_pixels BUT this dimension HAS NO DIMVAR! Code
> will probably fail at this point"
> # [1] "vobjtovarid4:  WARNING  I was asked to get a varid for
> dimension named number_of_lines BUT this dimension HAS NO DIMVAR! Code will
> probably fail at this point"
>
> Best regards,
>
> Frederico Faleiro
> Postdoctoral Researcher in the INCT-EECBio (https://www.eecbio.ufg.br/)
> Federal University of Goiás | Brazil
> RG: https://www.researchgate.net/profile/Frederico_Faleiro
>
> On Tue, Mar 17, 2020 at 11:39 AM Miluji Sb  wrote:
>
> > Dear all,
> >
> > Hope everyone is keeping safe.
> >
> > I am trying to extract/read VIIRS nighttime lights data but the output
> > seems rather strange.
> >
> > I have uploaded the file here
> >  so
> as
> > not exceed the size limit of the email.
> >
> > ## Code ##
> > library(ncdf4)
> > library(rgdal)
> >
> > netcdf_file <- c("VNP02DNB_NRT.A2020069.1048.001.nc")
> > nl <- brick(netcdf_file, lvar=0, values=TRUE,
> > varname="observation_data/DNB_observations")
> >
> > ## Detail ##
> > class  : RasterLayer
> > dimensions : 3232, 4064, 13134848  (nrow, ncol, ncell)
> > resolution : 1, 1  (x, y)
> > extent : 0.5, 4064.5, 0.5, 3232.5  (xmin, xmax, ymin, ymax)
> > crs: NA
> > source : C:/Users/shour/Desktop/VNP02DNB_NRT.A2020069.1048.001.nc
> > names  : DNB.observations.at.pixel.locations
> > zvar   : observation_data/DNB_observations
> >
> > The spatial resolution is supposed to be 750m but this shows 1°. What am
> I
> > doing wrong? Any help will be greatly appreciated. Thank you.
> >
> > Sincerely,
> >
> > Millu
> >
> > [[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
>

[[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] R-Tutorial: Geospatial Data Science in R (Resent)

[Zia Ahmed]

Hi All,

I have published an online tutorial related to spatial-data analysis with R
(https://zia207.github.io/geospatial-r-github.io/). Nothing is new here, I
have just organized several R-code and data that I have used in my several
publications. Most of the codes were written with the help of postings in
several online blogs: such as R-sig-Geo
<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>, Stack Overflow
<https://stackoverflow.com/>, and R bloggers
<https://www.r-bloggers.com/update-can-we-predict-flu-outcome-with-machine-learning-in-r/>
and on-line tutorials such as Spatial Data Science
<https://rspatial.org/index.html> and Geostatistics & Open-source
statistical computing
<http://www.css.cornell.edu/faculty/dgr2/teach/degeostats.html>.  I think
It would help someone we have no prior knowledge of GIS, remote sensing or
any other area of geoinformatics, but have some experience in R-coding. The
data used in this tutorial also available for download.



I appreciate any feedback for improving this tutorial.



I appreciate any feedback for improving this tutorial.



Best



Zia Ahmed

University at Buffalo



Tutorial consist following topics:



*1. **Spatial Data Processing*
<https://zia207.github.io/geospatial-r-github.io/about.html>

   - Reading and Writing Spatial Data
   
<https://zia207.github.io/geospatial-r-github.io/read-write-spatial-data.html>
  - Vector data
  - Raster data
   - Map Projection and Coordinate Reference Systems
   
<https://zia207.github.io/geospatial-r-github.io/map-projection-coordinate-reference-systems.html>
  - Geographic coordinate system (GCS)
  - Projected coordinate system
  - Coordinate Reference System in R
   - Geoprocessing of Vector data
   
<https://zia207.github.io/geospatial-r-github.io/geoprocessing-vector-data.html>
  - Clipping
  - Union
  - Dissolve
  - Intersect
  - Erase
  - Convex Hull
  - Buffer
   - Working with Spatial Point Data
   
<https://zia207.github.io/geospatial-r-github.io/working-with-spatial-point-data.html>
  - Create a Spatial Point Data Frame
  - Extract Environmental Covariates to SPDF
  - Create a Prediction Grid
  - Exploratory Data Analysis
  - Plot Data on Web Map
   - Working with Spatial Polygon Data
   
<https://zia207.github.io/geospatial-r-github.io/working-with-spatial-polygon.html>
  - Data Processing
  - Visualization
  - Animation of Time Series Data
   - Working with Raster Data
   
<https://zia207.github.io/geospatial-r-github.io/working-with-raster-data.html>
  - Basic Raster Operation
  - Clipping
  - Reclassification
  - Focal Statistics
  - Raster Algebra
  - Aggregation
  - Resample
  - Mosaic
  - Convert Raster to Point Data
  - Convert Point Data to Raster
  - Raster Stack and Raster Brick
  - Digital Terrain Modeling
 - Slope
 - Aspect
 - Hillshade
 - Terrain Ruggedness Index
 - Topographic Position Index
 - Roughness
 - Curvature
 - Flow Direction
  - netCDF Data Processing
   <https://zia207.github.io/geospatial-r-github.io/netCDF-data-processing.html>



*2. **Spatial Statistics*
<https://zia207.github.io/geospatial-r-github.io/spatial-statistics.html>

   - Spatial Autocorrelation
   
<https://zia207.github.io/geospatial-r-github.io/spatial-autocorrelation.html>
  - Moran’s I
  - Geary’s C
  - Getis’s Gi
   - Point Pattern Analysis
   <https://zia207.github.io/geospatial-r-github.io/point-pattern-analysis.html>
   - Geographically Weighted Mmodels
   
<https://zia207.github.io/geospatial-r-github.io/geographically-weighted-models.html>
  - Geographically Weighted Summary Statistics
  
<https://zia207.github.io/geospatial-r-github.io/geographically-weighted-summary-statistics.html>
  - Geographically Weighted Principal Components Analysis
  
<https://zia207.github.io/geospatial-r-github.io/geographically-weighted-principal-components-analysis.html>
  - Geographically Weighted Regression
  
<https://zia207.github.io/geospatial-r-github.io/geographically-weighted-regression.html>
 - Geographically Weighted OLS Regression
 
<https://zia207.github.io/geospatial-r-github.io/geographically-weighted-ols-regression.html>
 - Geographically Weighted Poisson Regression
 
<https://zia207.github.io/geospatial-r-github.io/geographically-weighted-poisson-regression.html>
 - Global and local (Geographically Weighted) Random Forest
 
<https://zia207.github.io/geospatial-r-github.io/geographically-wighted-random-forest.html>



*3. **Spatial Interpolation*
<https://zia207.github.io/geospatial-r-github.io/spatial-interpolation.html>



· Spatial Interpolation
<https://zia207.github.io/ge

[R-sig-Geo] R-tutorial: Geospatial Data Science with R

[Zia Ahmed]

Dear All,


Best


*Zia Ahmed, PhD*

Research Associate Professor (Data & Visualization)

RENEW (Research and Education in eNergy, Environment and Water) Institute
<http://www.buffalo.edu/renew.html>

University at Buffalo <http://www.buffalo.edu/>



Tutorial consist following topics:



*1. Spatial Data Processing
<https://zia207.github.io/geospatial-r-github.io/about.html>*

   - Reading and Writing Spatial Data
   
<https://zia207.github.io/geospatial-r-github.io/read-write-spatial-data.html>
  - Vector data
  - Raster data
   - Map Projection and Coordinate Reference Systems
   
<https://zia207.github.io/geospatial-r-github.io/map-projection-coordinate-reference-systems.html>
  - Geographic coordinate system (GCS)
  - Projected coordinate system
  - Coordinate Reference System in R
   - Geoprocessing of Vector data
   
<https://zia207.github.io/geospatial-r-github.io/geoprocessing-vector-data.html>
  - Clipping
  - Union
  - Dissolve
  - Intersect
  - Erase
  - Convex Hull
  - Buffer
   - Working with Spatial Point Data
   
<https://zia207.github.io/geospatial-r-github.io/working-with-spatial-point-data.html>
  - Create a Spatial Point Data Frame
  - Extract Environmental Covariates to SPDF
  - Create a Prediction Grid
  - Exploratory Data Analysis
  - Plot Data on Web Map
   - Working with Spatial Polygon Data
   
<https://zia207.github.io/geospatial-r-github.io/working-with-spatial-polygon.html>
  - Data Processing
  - Visualization
  - Animation of Time Series Data
   - Working with Raster Data
   
<https://zia207.github.io/geospatial-r-github.io/working-with-raster-data.html>
  - Basic Raster Operation
  - Clipping
  - Reclassification
  - Focal Statistics
  - Raster Algebra
  - Aggregation
  - Resample
  - Mosaic
  - Convert Raster to Point Data
  - Convert Point Data to Raster
  - Raster Stack and Raster Brick
  - Digital Terrain Modeling
 - Slope
 - Aspect
 - Hillshade
 - Terrain Ruggedness Index
 - Topographic Position Index
 - Roughness
 - Curvature
 - Flow Direction
  - netCDF Data Processing
   <https://zia207.github.io/geospatial-r-github.io/netCDF-data-processing.html>



*2. Spatial Statistics
<https://zia207.github.io/geospatial-r-github.io/spatial-statistics.html>*

   - Spatial Autocorrelation
   
<https://zia207.github.io/geospatial-r-github.io/spatial-autocorrelation.html>
  - Moran’s I
  - Geary’s C
  - Getis’s Gi
   - Point Pattern Analysis
   <https://zia207.github.io/geospatial-r-github.io/point-pattern-analysis.html>
   - Geographically Weighted Mmodels
   
<https://zia207.github.io/geospatial-r-github.io/geographically-weighted-models.html>
  - Geographically Weighted Summary Statistics
  
<https://zia207.github.io/geospatial-r-github.io/geographically-weighted-summary-statistics.html>
  - Geographically Weighted Principal Components Analysis
  
<https://zia207.github.io/geospatial-r-github.io/geographically-weighted-principal-components-analysis.html>
  - Geographically Weighted Regression
  
<https://zia207.github.io/geospatial-r-github.io/geographically-weighted-regression.html>
 - Geographically Weighted OLS Regression
 
<https://zia207.github.io/geospatial-r-github.io/geographically-weighted-ols-regression.html>
 - Geographically Weighted Poisson Regression
 
<https://zia207.github.io/geospatial-r-github.io/geographically-weighted-poisson-regression.html>
 - Global and local (Geographically Weighted) Random Forest
 
<https://zia207.github.io/geospatial-r-github.io/geographically-wighted-random-forest.html>



*3. Spatial Interpolation
<https://zia207.github.io/geospatial-r-github.io/spatial-interpolation.html>*



· Spatial Interpolation
<https://zia207.github.io/geospatial-r-github.io/spatial-interpolation.html>

oDeterministic Methods for Spatial Interpolation
<https://zia207.github.io/geospatial-r-github.io/deterministic-methods-for-spatial-interpolation.html>

§  Polynomial Trend Surface

§  Proximity Analysis-Thiessen Polygons

§  Nearest Neighbor Interpolation

§  Inverse Distance Weighted

§  Thin Plate Spline

oGeostatistical Methods for Spatial Interpolation
<https://zia207.github.io/geospatial-r-github.io/geostatistical-methods-for-spatial-interpolation.html>

§  Semivariogram Modeling
<https://zia207.github.io/geospatial-r-github.io/semivariogram-modeling.html>

§  Kriging <https://zia207.github.io/geospatial-r-github.io/kriging.html>

§  Ordinary Kriging
<https://zia207.github.io/geospatial-r-github.io/ordinary-kriging.html>

§  Universal Kriging
<https://zia207.github.io/geospatial-r-github.io/univer

Re: [R-sig-Geo] Plot a map in r (csv and shapefile)

[Zia Ahmed]

Please check this below link.  It might help you what you need to do.
Zia

https://zia207.github.io/geospatial-data-science.github.io/

On Wed, Jun 19, 2019, 6:40 AM Lara Silva  wrote:

> Hello,
>
> I am new in R and I need to plot the occurrence of a species.
> I have the presences of species (CSV) and the shapefile with the limits of
> study area.
>
> Which code should I use?
>
> Regards,
>
> Lara
>
> <
> https://www.avast.com/sig-email?utm_medium=email_source=link_campaign=sig-email_content=webmail
> >
> Sem
> vírus. www.avast.com
> <
> https://www.avast.com/sig-email?utm_medium=email_source=link_campaign=sig-email_content=webmail
> >
> <#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>
>
> [[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] Problems with SpatialPixelsDataFrame objects

[Zia Ahmed]

Dear List,
I am trying to calculate

*GLCM textures of all bands of a raster stack using GLCM package.  Using
following code it runs perfectly. The code created  four stacks for each
bands,  each stack contains for 8 raster objects,   but I am facing
difficulty  to get results for further analysis.  I appreciate if some one
help me out how to extract  four raster stacks from the results (glcm.all).*

*Thanks*

*Zia*



> library(glcm)> library(raster)> library(doParallel)> > names(L5TSR_1986)[1] 
> "b1" "b2" "b3" "b4"> > start.time <- Sys.time()> > foreach(rasname = 
> iter(names(L5TSR_1986)), .packages = "raster") %dopar% {+   glcm.all <-  
> glcm(L5TSR_1986[[rasname]])+   +   }[[1]]
class   : RasterStack
dimensions  : 167, 213, 35571, 8  (nrow, ncol, ncell, nlayers)
resolution  : 30, 30  (x, y)
extent  : 826245, 832635, 1107825, 1112835  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=16 +ellps=WGS84 +units=m +no_defs
names   :   glcm_mean, glcm_variance, glcm_homogeneity,
glcm_contrast, glcm_dissimilarity, glcm_entropy, glcm_second_moment,
glcm_correlation
min values  :  0.05902778,3.10906455,   0.02218517,
0., 0.,   0., 0.,
   -Inf
max values  :   0.6927083,   497.8806062,1.000,
118.556,  9.667,2.1972246,  1.000,
 Inf


[[2]]
class   : RasterStack
dimensions  : 167, 213, 35571, 8  (nrow, ncol, ncell, nlayers)
resolution  : 30, 30  (x, y)
extent  : 826245, 832635, 1107825, 1112835  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=16 +ellps=WGS84 +units=m +no_defs
names   :   glcm_mean, glcm_variance, glcm_homogeneity,
glcm_contrast, glcm_dissimilarity, glcm_entropy, glcm_second_moment,
glcm_correlation
min values  :  0.04687500,1.84165521,   0.03395896,
0., 0.,   0., 0.,
   -Inf
max values  :   0.736,   556.6514305,1.000,
94.778,  8.778,2.1972246,  1.000,
Inf


[[3]]
class   : RasterStack
dimensions  : 167, 213, 35571, 8  (nrow, ncol, ncell, nlayers)
resolution  : 30, 30  (x, y)
extent  : 826245, 832635, 1107825, 1112835  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=16 +ellps=WGS84 +units=m +no_defs
names   :   glcm_mean, glcm_variance, glcm_homogeneity,
glcm_contrast, glcm_dissimilarity, glcm_entropy, glcm_second_moment,
glcm_correlation
min values  :  0.03993056,0.91840278,   0.02775318,
0., 0.,   0., 0.,
   -Inf
max values  :   0.7291667,   572.3779206,1.000,
147.556, 10.889,2.1972246,  1.000,
 Inf


[[4]]
class   : RasterStack
dimensions  : 167, 213, 35571, 8  (nrow, ncol, ncell, nlayers)
resolution  : 30, 30  (x, y)
extent  : 826245, 832635, 1107825, 1112835  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=16 +ellps=WGS84 +units=m +no_defs
names   :   glcm_mean, glcm_variance, glcm_homogeneity,
glcm_contrast, glcm_dissimilarity, glcm_entropy, glcm_second_moment,
glcm_correlation
min values  :  0.07812500,3.79071422,   0.02196611,
0., 0.,   0., 0.,
   -Inf
max values  :   0.8975694,   803.7921278,1.000,
160.111, 10.667,2.1972246,  1.000,
 Inf

> > end.time <- Sys.time()> time.taken <- end.time - start.time> time.takenTime 
> > difference of 7.402473 secs


>


On Sun, Nov 20, 2016 at 4:50 AM, Pedro Perez  wrote:

> Hi everybody,
>
> The following two scripts will generate a "SpatialPixelDataFrame" object:
>
> # FIRST
> library(rgdal)
> elev.grid <- readGDAL("whatever.asc")
> elev.grid <- as(elev.grid, "SpatialPixelsDataFrame")
>
> # SECOND
> library(raster)
> library(SDMTools)
> library(adehabitat)
> elev.grid <- raster("whatever.asc")
> elev.grid.asc <- asc.from.raster(elev.grid)
> elev.grid.SPDF <- asc2spixdf(elev.grid.asc)
>
>
> HOWEVER, the first one excedes the capability of my computing
> resources when applying it to big (15000 x 16000) grids, and the
> second one generates an object which I can't use for some further
> analyses. For example, when I use it for krige purposes
>
> x <- krige(V3~var, points, elev.grid)
>
> I get the following:
>
> Error in model.frame.default(terms(formula), as(data, "data.frame"),
> na.action = na.fail) :
>   invalid type (closure) for variable 'var'
>
> I will be really thankful if somebody is kind enough to tell me how to
> fix it, whether providing me a trick to handle the memory/capability
> issue of the first case, or fixing the error generated by the second
> case.
>
> THANKS A LOT IN ADVANCE!!!
>
> Paolo
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>



-- 
Zia Uddin Ahmed, PhD
CIMMYT

[[alternative 

Re: [R-sig-Geo] Modify spplot frame borders, padding and font

[Zia Ahmed]

Hi Iain,

par.settings is a not attribute of spplot, it is an attribute of lattice graph. 
It works fine with ssplot. See the following examples.
Thanks, zia

library(lattice)

library(gstat)
data(meuse.grid)

gridded(meuse.grid) = ~x+y
names(meuse.grid)
spplot(meuse.grid[4],main = list(Soil types,cex=1.5),
par.settings=list(axis.line=list(col=NA)))

On 5/31/2012 3:17 AM, Edzer Pebesma wrote:
 library(lattice)
 trellis.par.set(axis.line=list(col=NA))
 library(sp)
 loadMeuse()
 spplot(meuse.grid[1])

 takes away the border, but not only around the map area, also around the
 legend colour block.

 On 05/30/2012 10:18 PM, Dillingham, Iain wrote:
 Thank you for the prompt reply. Forgive me, but how did you know that 
 par.settings was an attribute to spplot? I'm struggling with interpreting 
 the documentation. For example, if I wish to remove the outline and 
 background color from each frame title, where would I look for the 
 appropriate parameters and values?

 AFAICS, par.settings does nothing in the code fragment below.

 The spplot documentation says:

 Description:

   Lattice (trellis) plot methods for spatial data with attributes

 (and much later:)

   for useful values see the appropriate documentation of xyplot and
   levelplot.

 so reading into how lattice plots work, e.g. modifying parameters with
 trellis.par.set, or reading the documentation of xyplot and levelplot in
 package lattice might help. It assumes that the words lattice and
 trellis ring a bell.


 Thank you once again,

 Iain
 
 From: Zia Ahmed [z...@cornell.edu]
 Sent: 30 May 2012 20:08
 To: Dillingham, Iain; sig-geo
 Subject: Re: [R-sig-Geo] Modify spplot frame borders, padding and font

 You can add following line to get a map without frame.
 par.settings=list(axis.line=list(col=NA))
 For tittle,
 main=list((a) tittle,cex=.6)
 For size and space between frames (for nine frames):
 windows(width=5.5, height=8)
 tiff( file=Fig.tif,width=5.5, height=8,units = in, pointsize = 12, 
 res=1600,
 restoreConsole = T,compression =  lzw,bg=transparent)

 print(a,position = c(0.05,0.5,1,1),split=c(1,1,4,2), 
 more=TRUE,panel.width=list(1.5, inches),panel.height=list(1.5, inches))
 print(b,position = 
 c(0.15,0.5,1,1),split=c(2,1,4,2),more=TRUE,panel.width=list(1.5, 
 inches),panel.height=list(1.5, inches))
 print(c,position = c(0.45,0.5,1,1), split=c(3,1,4,2), 
 more=TRUE,panel.width=list(1.5, inches),panel.height=list(1.5, inches))
 print(d,position = c(0.05,0.5,1,0.65),split=c(1,1,4,2), 
 more=TRUE,panel.width=list(1.5, inches),panel.height=list(1.5, inches))
 print(e,position = 
 c(0.15,0.5,1,0.65),split=c(2,1,4,2),more=TRUE,panel.width=list(1.5, 
 inches),panel.height=list(1.5, inches))
 print(f,position = c(0.45,0.5,1,0.65), split=c(3,1,4,2), 
 more=TRUE,panel.width=list(1.5, inches),panel.height=list(1.5, inches))
 print(g,position = c(0.05,0.5,1,0.30),split=c(1,1,4,2), 
 more=TRUE,panel.width=list(1.5, inches),panel.height=list(1.5, inches))
 print(h,position = 
 c(0.15,0.5,1,0.30),split=c(2,1,4,2),more=TRUE,panel.width=list(1.5, 
 inches),panel.height=list(1.5, inches))
 print(i,position = c(0.45,0.5,1,0.30), split=c(3,1,4,2), 
 more=TRUE,panel.width=list(1.5, inches),panel.height=list(1.5, inches))
 print(j,position = c(0.05,0,1,0.15),split=c(1,1,4,2), 
 more=TRUE,panel.width=list(1.5, inches),panel.height=list(1.5, inches))
 print(k,position = 
 c(0.15,0,1,0.15),split=c(2,1,4,2),more=TRUE,panel.width=list(1.5, 
 inches),panel.height=list(1.5, inches))
 print(l,position = c(0.45,0,1,0.15), split=c(3,1,4,2), 
 more=TRUE,panel.width=list(1.5, inches),panel.height=list(1.5, inches))
 dev.off()

 On 5/30/2012 2:55 PM, Dillingham, Iain wrote:

 Hello everyone,

 I have used spplot to draw six small-multiple choropleth maps of London. 
 Although I am happy with the result, I would like to recolour (and 
 potentially remove) the outline from each frame; to add some padding between 
 each map and its containing frame; and finally to change the font used to 
 display the heading in each frame.

 I have searched the documentation. However, I am new to R and I have been 
 unable to uncover the mechanism that would allow me to modify these 
 attributes. Any assistance would be gratefully received!

 Many thanks,

 Iain

 --
 Iain Dillingham
 School of Informatics
 City University London
 EC1V 0HB
 http://www.soi.city.ac.uk/~abhm372/
 ___
 R-sig-Geo mailing list
 R-sig-Geo@r-project.orgmailto:R-sig-Geo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-geo



 --
 -
 Zia Ahmed, PhD
 Research Associate
 Department of Crop and Soil Sciences
 1002 Bradfield Hall, Cornell University
 Ithaca, NY 14853-4203
 t. 607.255.9387
 f. 607.255.3207
 email z...@cornell.edumailto:z...@cornell.edu

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

[R-sig-Geo] Asking Help again: Re: Help: Create several ESRI Shape files in a Loop

[Zia Ahmed]

I like to write a r-function to create 100 shape files in a loop.I 
have three data files. Each file has 100 data fields (like sim1, sim2, 
sim3sim100) with a unique ID fields (Thana_ID). First I want to
create 100 data frame from these files and then merge or join each of 
this 100 files with a polygon files individually to create 100 ESRI 
shape files.   I am trying to do this with  two steps. First,  I have 
created  100 csv files. But,  I when I am trying to merge these  files 
with a shape file (Step 2),  it show error massage.  It create only 
three shape files (sim1, sim2, sim3). I have not able figured out  the 
source of error in code. Any suggestion?  Thanks

Error in writeOGR(thana, dsn = /GWPR/Shape_Files, paste(sim, i, 
.shp,  :
   Non-unique field names

#---
library(spdep)
library(gstat)
library(maptools)
library(maps)
library(rgdal)
library(foreign)
##
setwd(H:/GWPR/csv_files)

id - read.csv(id_file.csv, as.is = TRUE)
gas - read.csv('thana_sis_logGAS_intake_246.csv', as.is = TRUE)
total - read.csv(thana_sis_logTotal_intake_246.csv, as.is = TRUE)
was - read.csv(thana_sis_logWAS_intake_246.csv, as.is = TRUE)

thana - readOGR(., bd_thana_246) #
OGR data source with driver: ESRI Shapefile
Source: ., layer: bd_thana_246
with 246 features and 1 fields
Feature type: wkbPolygon with 2 dimensions
  str(thana)
Formal class 'SpatialPolygonsDataFrame' [package sp] with 5 slots
   ..@ data   :'data.frame': 246 obs. of  1 variable:
   .. ..$ THANA_ID: num [1:246] 4 20 31 37 171 46 282 305 463 67 ...
   ..@ polygons   :List of 246


  # Step: 1: CSV files

for (i in 1:100){
key - paste('sim', i, sep = )
write.csv(data.frame(THANA_ID = id$THANA_ID
, PatNum = id$PatNum
, LogWAS = was[[key]]
, LogGAS = gas[[key]]
, LogTotal = total[[key]]
, stringsAsFactors = FALSE
  )
  , paste(key, '.csv', sep = '')
   )}

  # Step: 2: Merge with

for (i in 1:100){
 x - read.csv(paste('sim', i, '.csv', sep = ''), as.is = TRUE)
 thana@data - merge(thana@data, x, by.x=THANA_ID,by.y=THANA_ID, 
all.x = TRUE,
sort = FALSE)
 writeOGR(thana, dsn = /GWPR/Shape_Files, paste('sim',  i,.shp, 
sep =
''), driver=ESRI Shapefile)
}

-- 
-
Zia Ahmed, PhD
Research Associate
Department of Crop and Soil Sciences
1002 Bradfield Hall, Cornell University
Ithaca, NY 14853-4203
t. 607.255.9387
f. 607.255.3207
email z...@cornell.edu



[[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] Help: Create several ESRI Shape files in a Loop

[Zia Ahmed]

I like to write a r-function to create 100 shape files in a loop.I 
have three data files. Each file has 100 data fields (like sim1, sim2, 
sim3sim100) with a unique ID fields (Thana_ID). First I want to 
create 100 data frame from these files and then merge or join each of 
this 100 files with a polygon files individually to create 100 ESRI 
shape files. I am trying to this  following way. But I need do this in a 
loop.  Any suggestion will be appreciated.

Thanks
Zia

#--

library(sp)
library(spdep)
library(gstat)
library(maptools)
library(maps)
library(rgdal)

##
setwd(H:/GWPR)

gas-read.csv(thana_sis_logGAS_intake_246.csv,header=TRUE)
was-read.csv(thana_sis_logWAS_intake_246.csv,header=TRUE)
total-read.csv(thana_sis_logTotal_intake_246.csv,header=TRUE)
names(thana)
names(was)
names(gas)
names(total

# Create 100 data-frames by selecting one column from each file.

sim1-cbind(THANA_ID=thana$THANA_ID, logWAS=was$sim1, logGAs=gas$sim1, 
logTotal=total$sim1)

sim2-cbind(THANA_ID=thana$THANA_ID, logWAS=was$sim2, logGAs=gas$sim2, 
logTotal=total$sim2)
.
.
.
sim100-cbind(THANA_ID=thana$THANA_ID, logWAS=was$sim100, 
logGAs=gas$sim100, logTotal=total$sim100)



#Join each  data framewith a shape file and then write #ESRIShape files.
# Total number of shapefiles willbe100.

sim1.thana-readShapePoly(bd_thana_246.shp)

sim1.thana@data - merge(thana@data,sim1,by.x=THANA_ID, 
by.y=THANA_ID, all.x=T, sort=F)

writeOGR(sim1.thana,dsn=/GWPR/Shape_Files, layer=sim1.thana, 
driver=ESRI Shapefile)

sim2.thana-readShapePoly(bd_thana_246.shp)

sim2.thana@data - merge(thana@data,sim2,by.x=THANA_ID, 
by.y=THANA_ID, all.x=T, sort=F)

writeOGR(sim2.thana,dsn=/GWPR/Shape_Files, layer=sim2.thana, 
driver=ESRI Shapefile)

.

.

.

sim100.thana-readShapePoly(bd_thana_246.shp)

sim100.thana@data - merge(thana@data,sim100,by.x=THANA_ID, 
by.y=THANA_ID, all.x=T, sort=F)

writeOGR(sim100.thana,dsn=/GWPR/Shape_Files, layer=sim100.thana, 
driver=ESRI Shapefile)



-- 
-
Zia Ahmed, PhD
Research Associate
Department of Crop and Soil Sciences
1002 Bradfield Hall, Cornell University
Ithaca, NY 14853-4203
t. 607.255.9387
f. 607.255.3207
email z...@cornell.edu



[[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] Help: Create several ESRI Shape files in a Loop

[Zia Ahmed]

I have tried your code, but it create only one shape file (sim100.shp). 
Thanks
Zia

for (i in 1:100)
sim-cbind(THANA_ID=thana$THANA_ID, logWAS=was[,i], 
logGAs=gas[,i],logTotal=total[,i])
sim.shp - thana
sim.shp@data - merge(thana@data,sim,by.x=THANA_ID,by.y=THANA_ID, 
all.x=T, sort=F)
writeOGR(sim.shp,dsn=/GWPR/Shape_Files, 
layer=paste(sim_,i,.shp,sep=),driver=ESRI Shapefile)



On 3/12/2012 10:20 AM, Raphael Saldanha wrote:
 Ops! Take off the $ on this line
 sim-cbind(THANA_ID=thana$THANA_ID, logWAS=was$[,i], 
 logGAs=gas$[,i],logTotal=total$[,i])

 This mus be:
 sim-cbind(THANA_ID=thana$THANA_ID, logWAS=was[,i], 
 logGAs=gas[,i],logTotal=total[,i])


 On Mon, Mar 12, 2012 at 11:13 AM, Raphael Saldanha 
 saldanha.plan...@gmail.com mailto:saldanha.plan...@gmail.com wrote:

 Hi Zia,

 There are two methods were you can work: the for statement and the
 tapply function.

 First, take a good look on these: ?for and ?tapply

 You can try something like this:

 for (i in 1:100)
 sim-cbind(THANA_ID=thana$THANA_ID, logWAS=was$[,i],
 logGAs=gas$[,i],logTotal=total$[,i])
 sim.shp - thana
 sim.shp@data -
 merge(thana@data,sim,by.x=THANA_ID,by.y=THANA_ID, all.x=T, sort=F)
 writeOGR(sim.shp,dsn=/GWPR/Shape_Files,
 layer=paste(sim_,i,.shp,sep=),driver=ESRI Shapefile)

 A better way is use create a list of variables, so i will assume
 the name of the variable.


 On Mon, Mar 12, 2012 at 10:53 AM, Zia Ahmed z...@cornell.edu
 mailto:z...@cornell.edu wrote:

 I like to write a r-function to create 100 shape files in a
 loop.I
 have three data files. Each file has 100 data fields (like
 sim1, sim2,
 sim3sim100) with a unique ID fields (Thana_ID). First I
 want to
 create 100 data frame from these files and then merge or join
 each of
 this 100 files with a polygon files individually to create 100
 ESRI
 shape files. I am trying to this  following way. But I need do
 this in a
 loop.  Any suggestion will be appreciated.

 Thanks
 Zia

 #--

 library(sp)
 library(spdep)
 library(gstat)
 library(maptools)
 library(maps)
 library(rgdal)

 ##
 setwd(H:/GWPR)

 gas-read.csv(thana_sis_logGAS_intake_246.csv,header=TRUE)
 was-read.csv(thana_sis_logWAS_intake_246.csv,header=TRUE)
 total-read.csv(thana_sis_logTotal_intake_246.csv,header=TRUE)
 names(thana)
 names(was)
 names(gas)
 names(total

 # Create 100 data-frames by selecting one column from each file.

 sim1-cbind(THANA_ID=thana$THANA_ID, logWAS=was$sim1,
 logGAs=gas$sim1,
 logTotal=total$sim1)

 sim2-cbind(THANA_ID=thana$THANA_ID, logWAS=was$sim2,
 logGAs=gas$sim2,
 logTotal=total$sim2)
 .
 .
 .
 sim100-cbind(THANA_ID=thana$THANA_ID, logWAS=was$sim100,
 logGAs=gas$sim100, logTotal=total$sim100)



 #Join each  data framewith a shape file and then write
 #ESRIShape files.
 # Total number of shapefiles willbe100.

 sim1.thana-readShapePoly(bd_thana_246.shp)

 sim1.thana@data - merge(thana@data,sim1,by.x=THANA_ID,
 by.y=THANA_ID, all.x=T, sort=F)

 writeOGR(sim1.thana,dsn=/GWPR/Shape_Files, layer=sim1.thana,
 driver=ESRI Shapefile)

 sim2.thana-readShapePoly(bd_thana_246.shp)

 sim2.thana@data - merge(thana@data,sim2,by.x=THANA_ID,
 by.y=THANA_ID, all.x=T, sort=F)

 writeOGR(sim2.thana,dsn=/GWPR/Shape_Files, layer=sim2.thana,
 driver=ESRI Shapefile)

 .

 .

 .

 sim100.thana-readShapePoly(bd_thana_246.shp)

 sim100.thana@data - merge(thana@data,sim100,by.x=THANA_ID,
 by.y=THANA_ID, all.x=T, sort=F)

 writeOGR(sim100.thana,dsn=/GWPR/Shape_Files,
 layer=sim100.thana,
 driver=ESRI Shapefile)



 --
 -
 Zia Ahmed, PhD
 Research Associate
 Department of Crop and Soil Sciences
 1002 Bradfield Hall, Cornell University
 Ithaca, NY 14853-4203
 t. 607.255.9387
 f. 607.255.3207
 email z...@cornell.edu mailto:z...@cornell.edu



[[alternative HTML version deleted]]

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




 -- 
 Atenciosamente,

 Raphael Saldanha
 saldanha.plan...@gmail.com mailto:saldanha.plan...@gmail.com




 -- 
 Atenciosamente,

 Raphael Saldanha
 saldanha.plan...@gmail.com mailto:saldanha.plan...@gmail.com

-- 
-
Zia Ahmed, PhD

[R-sig-Geo] Population weighted Centroids?

[Zia Ahmed]

Is there anyway to calculate population weighted centroids of a spatial 
polygon   in R? Any idea will be appreciated. Thanks

Zia


--
-
Zia Ahmed, PhD
Research Associate
Department of Crop and Soil Sciences
1002 Bradfield Hall, Cornell University
Ithaca, NY 14853-4203
t. 607.255.9387
f. 607.255.3207
email z...@cornell.edu

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

[R-sig-Geo] GGWR (family = Poission): Prediction of unknown location

[Zia Ahmed]

Hello,
I am trying to predict expected count for a geographical location  (file 
name:  data_181 - test data)  from a another data set (data_65 - 
work_data)  using GGWR poission model.  Following R code I used.  But No 
predicted counts were  found  in output SDF  file.  I am new in GGWR.   
Help will be appreciated. I attach my data sets here. Thanks

Zia

library(spgwr)
data(data_65)   # work data set
data(data_181) #  test data set

bw - ggwr.sel(PopRisk_100TH ~ log10(WAS)*log10(GAS)+ 
offset(log(population_1000)), coords=cbind(data_65$x, data_65$y),

  data=data_65, family=poisson())
coordinates(data.181)-~x+y
ggwr.181 - ggwr(PopRisk_100TH ~ log10(WAS)*log10(GAS) + 
offset(log(population_1000)), coords=cbind(data.65$x, data.65$y),
  data=data.65,family=poisson(),bandwidth=bw,type = 
c(response),fit.points = data.181)


 ggwr.181$SDF


	sum.w 	X.Intercept. 	log10.wasMean. 	log10.gasMean. 
log10.wasMean..log10.gasMean. 	dispersion 	response_resids 	x 	y
1 	62.50764 	-2.36702 	1.45678 	-2.42768 	5.031622 	1 	NA 	2779558 
531287.3
2 	63.20595 	-2.36062 	1.451657 	-2.39074 	4.995947 	1 	NA 	2767150 
575319.5
3 	62.87379 	-2.36683 	1.456869 	-2.42598 	5.026145 	1 	NA 	2798559 
563289.6
4 	63.13653 	-2.36322 	1.453891 	-2.4054 	5.007724 	1 	NA 	2784218 
576455.7

5   62.55473-2.369911.459364-2.44354
5.0425961   NA  2807095 546223
6 	62.85015 	-2.36264 	1.453169 	-2.40282 	5.009307 	1 	NA 	2762510 
547772.6
7 	62.67004 	-2.36523 	1.455308 	-2.4175 	5.022374 	1 	NA 	2773146 
538941.8
8 	62.40817 	-2.36266 	1.452969 	-2.40357 	5.013347 	1 	NA 	2745574 
520650.3
9 	62.02654 	-2.36357 	1.453601 	-2.40908 	5.020167 	1 	NA 	2739576 
502023.7
10 	63.01464 	-2.35797 	1.449239 	-2.37628 	4.9865 	1 	NA 	2739005 
557044.6
11 	62.8022 	-2.35943 	1.450378 	-2.38489 	4.995307 	1 	NA 	2739341 
542981.4
12 	62.95506 	-2.36003 	1.45098 	-2.38804 	4.996592 	1 	NA 	2749492 
553074.4



--
-
Zia Ahmed, PhD
Research Associate
Department of Crop and Soil Sciences
1002 Bradfield Hall, Cornell University
Ithaca, NY 14853-4203
t. 607.255.9387
f. 607.255.3207
email z...@cornell.edu




data_65.csv
Description: MS-Excel spreadsheet


data_181.csv
Description: MS-Excel spreadsheet
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Problem with create several maps and save a PNG files in a loop from a sppolygon object

[Zia Ahmed]

I am trying to create several maps and save as PNG files in a loop from 
a sppolygon  object. This object contains   464 records  and 108 
fields.  I  want to use column no. from 7 to 108 (see below) to create 
100 maps and want to save save as png  files. I tried following code but 
it  showed  error.  I used this code successfully for sppoint data 
frame. May be I am  missing something here.
Thanks
Zia

 for (i in 1:ns){
+ png(filename = paste(sim,i,.png,sep=), width=400, height=500, 
units=px,
+ pointsize=12, bg=transparent, res=NA)
+ print(spplot(sim[i], main = list(label=paste(Realization,i),cex=1.5),
+ par.settings=list(axis.line=list(col=NA)),at=at,
+ col.regions=rgb.palette (21),
+ colorkey=list(space=right,width=2,at=1:22,labels=list(cex=1.2,at=1:22,
+ labels=c(, 
104080120160,,,  
200, )))
+ ))
+ dev.off()
+ }
Error in print(spplot(sim[i], main = list(label = paste(Realization,  :
   error in evaluating the argument 'x' in selecting a method for 
function 'print': Error in .local(obj, ...) :
   length of col.regions should match number of factor levels
 

  summary(sim)
Object of class SpatialPolygonsDataFrame
Coordinates:
 min   max
x 2542038.2 3016189.2
y  322929.6  985549.3
Is projected: TRUE
proj4string : [+init=epsg:3106]
Data attributes:


  names(sim)
   [1] ID DIV_NAME   THANA_ID   x-cent  y-cent
   [6] DIST_NAME  THANA_NAME sim1   sim2   sim3
  [11] sim4   sim5   sim6   sim7   sim8
  [16] sim9   sim10  sim11  sim12  sim13
  [21] sim14  sim15  sim16  sim17  sim18
  [26] sim19  sim20  sim21  sim22  sim23
  [31] sim24  sim25  sim26  sim27  sim28
  [36] sim29  sim30  sim31  sim32  sim33
  [41] sim34  sim35  sim36  sim37  sim38
  [46] sim39  sim40  sim41  sim42  sim43
  [51] sim44  sim45  sim46  sim47  sim48
  [56] sim49  sim50  sim51  sim52  sim53
  [61] sim54  sim55  sim56  sim57  sim58
  [66] sim59  sim60  sim61  sim62  sim63
  [71] sim64  sim65  sim66  sim67  sim68
  [76] sim69  sim70  sim71  sim72  sim73
  [81] sim74  sim75  sim76  sim77  sim78
  [86] sim79  sim80  sim81  sim82  sim83
  [91] sim84  sim85  sim86  sim87  sim88
  [96] sim89  sim90  sim91  sim92  sim93
[101] sim94  sim95  sim96  sim97  sim98
[106] sim99  sim100

-- 
-
Zia Ahmed, PhD
Research Associate
Department of Crop and Soil Sciences
1002 Bradfield Hall, Cornell University
Ithaca, NY 14853-4203
t. 607.255.9387
f. 607.255.3207
email z...@cornell.edu



[[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] Thanks! Re: Help: create an animated GIF file from PNG files

[Zia Ahmed]

Hi Yihui,
Thank you so much for your help. It works using following function:
Zia

saveGIF(
for (i in 1:ns){
print(spplot(sim_ns[,i], main = list(label=paste(SIM,i),cex=1.5),
colorkey = list(space=right,height=1, width=1.5,labels=list(cex=1.5)),
par.settings=list(axis.line=list(col=NA)),at=at.pred,
col.regions=rgb.palette(100)))},
height = 600, width = 450, interval = .5, outdir = getwd())
}

On 1/25/2012 5:39 PM, Yihui Xie wrote:
 You can try the animation package with the function saveGIF(); see the
 help page for more info. Basically you just put the code in

 saveGIF({
whatever code that produces plots
 })

 and you do not need to take care of the issues like opening the
 graphical devices; the process is automatic.

 Regards,
 Yihui
 --
 Yihui Xiexieyi...@gmail.com
 Phone: 515-294-2465 Web: http://yihui.name
 Department of Statistics, Iowa State University
 2215 Snedecor Hall, Ames, IA



 On Wed, Jan 25, 2012 at 3:44 PM, Zia Ahmedz...@cornell.edu  wrote:
 I am trying to create an animated GIF file from a set of PNG files
 created from number realization.  Itried with animation ( depend -
 ImageMagick) and caTools packages, but I did not successful to make the
 GIS file.  Is  there any way to do this?
 Thanks
 Zia

 # GIF animation

 library(sp)
 library(gstat)
 library(lattice)
 library(RColorBrewer)
 library(maptools)
 library(maps)
 library(ggplot2)
 library(lubridate)
 library(animation)
 library(fields) # for tim.colors
 library(caTools) # for write.gif
 library(raster)

 data(meuse)
 data(meuse.grid)
 coordinates(meuse)- ~x+y
 coordinates(meuse.grid)- ~x+y
 gridded(meuse.grid)- TRUE

 lzn.vgm- variogram(log(zinc)~1, meuse)
 lzn.fit- fit.variogram(lzn.vgm, model = vgm(1, Sph, 900, 1))

 ns=4

 sim_ns = krige(log(zinc)~1, meuse, meuse.grid, model = lzn.fit,
 nmax = 30, nsim = ns)

 # Save figures as PNG files

 rgb.palette- colorRampPalette(c(lightyellow,green,orange,brown),
 space = rgb)

 setwd(C:/sim)

 for (i in 1:ns){
 png(filename = paste(sim,i,.png,sep=), width=400, height=500,
 units=px,
 pointsize=12, bg=transparent, res=NA)
 print(spplot(sim_ns[,i], main = list(label=paste(SIS,i),cex=1.5),
 par.settings=list(axis.line=list(col=NA)),
 col.regions=rgb.palette(22)))
 dev.off()
 }

 # Creat a gif.file
 system(convert -delay 40 *.png sim.gif) #
 Source(http://ryouready.wordpress.com/2010/11/21/animate-gif-images-in-r-imagemagick/)
 Invalid Parameter - 40
 Warning message:
 running command 'convert -delay 40 sim*.png plot.gif', invisible = F,
 wait=T' had status 4



 -
 Zia Ahmed, PhD
 Research Associate
 Department of Crop and Soil Sciences
 1002 Bradfield Hall, Cornell University
 Ithaca, NY 14853-4203
 t. 607.255.9387
 f. 607.255.3207
 email z...@cornell.edu



 [[alternative HTML version deleted]]

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

-- 
-
Zia Ahmed, PhD
Research Associate
Department of Crop and Soil Sciences
1002 Bradfield Hall, Cornell University
Ithaca, NY 14853-4203
t. 607.255.9387
f. 607.255.3207
email z...@cornell.edu



[[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] Help: create an animated GIF file from PNG files

[Zia Ahmed]

I am trying to create an animated GIF file from a set of PNG files 
created from number realization.  Itried with animation ( depend -  
ImageMagick) and caTools packages, but I did not successful to make the 
GIS file.  Is  there any way to do this?
Thanks
Zia

# GIF animation

library(sp)
library(gstat)
library(lattice)
library(RColorBrewer)
library(maptools)
library(maps)
library(ggplot2)
library(lubridate)
library(animation)
library(fields) # for tim.colors
library(caTools) # for write.gif
library(raster)

data(meuse)
data(meuse.grid)
coordinates(meuse) - ~x+y
coordinates(meuse.grid) - ~x+y
gridded(meuse.grid) - TRUE

lzn.vgm - variogram(log(zinc)~1, meuse)
lzn.fit - fit.variogram(lzn.vgm, model = vgm(1, Sph, 900, 1))

ns=4

sim_ns = krige(log(zinc)~1, meuse, meuse.grid, model = lzn.fit,
nmax = 30, nsim = ns)

# Save figures as PNG files

rgb.palette - colorRampPalette(c(lightyellow,green,orange,brown),
space = rgb)

setwd(C:/sim)

for (i in 1:ns){
png(filename = paste(sim,i,.png,sep=), width=400, height=500, 
units=px,
pointsize=12, bg=transparent, res=NA)
print(spplot(sim_ns[,i], main = list(label=paste(SIS,i),cex=1.5),
par.settings=list(axis.line=list(col=NA)),
col.regions=rgb.palette(22)))
dev.off()
}

# Creat a gif.file
system(convert -delay 40 *.png sim.gif) # 
Source(http://ryouready.wordpress.com/2010/11/21/animate-gif-images-in-r-imagemagick/)
Invalid Parameter - 40
Warning message:
running command 'convert -delay 40 sim*.png plot.gif', invisible = F, 
wait=T' had status 4



-
Zia Ahmed, PhD
Research Associate
Department of Crop and Soil Sciences
1002 Bradfield Hall, Cornell University
Ithaca, NY 14853-4203
t. 607.255.9387
f. 607.255.3207
email z...@cornell.edu



[[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] Calculating 95th percentile within polygons

[Zia Ahmed]

 I am trying to calculating  95th percentile within polygons  from a of 
set realizations - something like zonal statistics.

How do I calculate 95 th percentile for each polygon over all realizations.
Thanks
Zia

For example:

data(meuse)
data(meuse.grid)
coordinates(meuse) - ~x+y
coordinates(meuse.grid) - ~x+y

# Simulation
nsim=10
x - krige(log(zinc)~1, meuse, meuse.grid, model = vgm(.59, Sph, 874, 
.04), nmax=10, nsim=nsim)

over(sr, x[,1:4], fn = mean)

 over(sr, x[,1:4], fn = mean)
   sim1 sim2 sim3 sim4
r1 5.858169 5.792870 5.855246 5.868499
r2 5.588570 5.452744 5.596648 5.516289
r3 5.798087 5.860750 5.784194 5.848194
r4   NA   NA   NA   NA

# 95 th percentile at prediction grid:
x-as.data.frame(x)
y95-apply(x[3:nsim],1,stats::quantile,probs = 0.95,na.rm=TRUE) # 95 th 
percentile at each prediction grid




On 10/18/2011 3:30 PM, Edzer Pebesma wrote:

require(sp)
r1 = cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409,
180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676,
332618, 332413, 332349))
r2 = cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437,
179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683,
331133, 331623, 332152, 332357, 332373))
r3 = cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875,
179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004,
329783, 329665, 329720, 329933, 330478, 331062, 331086))
r4 = cbind(c(180304, 180403,179632,179420,180304),
c(332791, 333204, 333635, 333058, 332791))

sr1=Polygons(list(Polygon(r1)),r1)
sr2=Polygons(list(Polygon(r2)),r2)
sr3=Polygons(list(Polygon(r3)),r3)
sr4=Polygons(list(Polygon(r4)),r4)
sr=SpatialPolygons(list(sr1,sr2,sr3,sr4))
srdf=SpatialPolygonsDataFrame(sr, data.frame(cbind(1:4,5:2), 
row.names=c(r1,r2,r3,r4)))


data(meuse)
coordinates(meuse) = ~x+y

plot(meuse)
polygon(r1)
polygon(r2)
polygon(r3)
polygon(r4)
# retrieve mean heavy metal concentrations per polygon:
# attribute means over each polygon, NA for empty
over(sr, meuse[,1:4], fn = mean)

# return the number of points in each polygon:
sapply(over(sr, geometry(meuse), returnList = TRUE), length) 
attachment: zua3.vcf___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Calculating 95th percentile within polygons

[Zia Ahmed]

 Your are correct. Now, what I am  going to do: first calculate 95th 
percentile  of all realizations at prediction grid and then calculate 
the mean for each  polygon.  Is it possible to calculated weighted mean 
of the this  95th percentile? weight may be  area of polygon  or any 
value, for example population density  for each polygon.

Thanks
Zia

On 10/18/2011 5:12 PM, Edzer Pebesma wrote:

Zia,

You may want to rethink the question. Each realization has a 95 
percentile within a particular polygon. Over the set of realizations 
of some aggregated value for a polygon, you can take a 95 percentile.


These are two different things. The first is a spatial aggregation, 
the second an aggregation over the (sampled) probability distribution.


On 10/18/2011 11:07 PM, Zia Ahmed wrote:

I am trying to calculating 95th percentile within polygons from a of set
realizations - something like zonal statistics.
How do I calculate 95 th percentile for each polygon over all 
realizations.

Thanks
Zia

For example:

data(meuse)
data(meuse.grid)
coordinates(meuse) - ~x+y
coordinates(meuse.grid) - ~x+y

# Simulation
nsim=10
x - krige(log(zinc)~1, meuse, meuse.grid, model = vgm(.59, Sph, 874,
.04), nmax=10, nsim=nsim)
over(sr, x[,1:4], fn = mean)

 over(sr, x[,1:4], fn = mean)
sim1 sim2 sim3 sim4
r1 5.858169 5.792870 5.855246 5.868499
r2 5.588570 5.452744 5.596648 5.516289
r3 5.798087 5.860750 5.784194 5.848194
r4 NA NA NA NA

# 95 th percentile at prediction grid:
x-as.data.frame(x)
y95-apply(x[3:nsim],1,stats::quantile,probs = 0.95,na.rm=TRUE) # 95 th
percentile at each prediction grid



On 10/18/2011 3:30 PM, Edzer Pebesma wrote:

require(sp)
r1 = cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409,
180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676,
332618, 332413, 332349))
r2 = cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437,
179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683,
331133, 331623, 332152, 332357, 332373))
r3 = cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875,
179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004,
329783, 329665, 329720, 329933, 330478, 331062, 331086))
r4 = cbind(c(180304, 180403,179632,179420,180304),
c(332791, 333204, 333635, 333058, 332791))

sr1=Polygons(list(Polygon(r1)),r1)
sr2=Polygons(list(Polygon(r2)),r2)
sr3=Polygons(list(Polygon(r3)),r3)
sr4=Polygons(list(Polygon(r4)),r4)
sr=SpatialPolygons(list(sr1,sr2,sr3,sr4))
srdf=SpatialPolygonsDataFrame(sr, data.frame(cbind(1:4,5:2),
row.names=c(r1,r2,r3,r4)))

data(meuse)
coordinates(meuse) = ~x+y

plot(meuse)
polygon(r1)
polygon(r2)
polygon(r3)
polygon(r4)
# retrieve mean heavy metal concentrations per polygon:
# attribute means over each polygon, NA for empty
over(sr, meuse[,1:4], fn = mean)

# return the number of points in each polygon:
sapply(over(sr, geometry(meuse), returnList = TRUE), length)


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

Re: [R-sig-Geo] merging several shapefiles into one

[Zia Ahmed]

 I used following code to merge a series of  shape files with similar 
file structure. We can try this, but you need to do some modification of 
the code.

Zia

###

# Require packages: rgdal and maptool
#-

library(rgdal)
library(maptools)

# Get list of all shape files from county directories/sptial/
#

files - list.files(pattern=soilmu_a_ia.*.shp$, 
recursive=TRUE,full.names=TRUE)

uid-1

# Get polygons from first file (_001)
#-

poly.data- readOGR(files[1],gsub(^.*/(.*).shp$, \\1, files[1]))
n - length(slot(poly.data, polygons))
poly.data - spChFIDs(poly.data, as.character(uid:(uid+n-1)))
uid - uid + n

# mapunit polygoan: combin remaining  polygons with first polygoan
#-

for (i in 2:length(files)) {
temp.data - readOGR(files[i], gsub(^.*/(.*).shp$, \\1,files[i]))
n - length(slot(temp.data, polygons))
temp.data - spChFIDs(temp.data, as.character(uid:(uid+n-1)))
uid - uid + n
poly.data - spRbind(poly.data,temp.data)
}

names(poly.data)
proj4string(poly.data)

On 5/25/2011 3:22 AM, Pierre Roudier wrote:

Hi Nicolas,

I do not have a R session handy, but if I remember correctly, the
function gUnion() in the rgeos package would do what you want.

Hope this helps,

Pierre

2011/5/25 Nicolas Degalliernicolas.degall...@ird.fr:

Hi All,

I have read different shapefiles that are contiguous and that have different column names 
with the function readShapePoly.

example :

noumea- readShapePoly(noumea.file, IDvar = NOUMEA2_ID, proj4string = 
crs_rgnc)

and

dumbea- readShapePoly(dumbea.file, IDvar = DUMBEA2_, proj4string = crs_rgnc)

Now I would like to create a new shapefile with all the polygones in the noumea and 
dumbea spatial polygone data frames.

I tried :
spRbind(noumea, dumbea)

I get this message error :

Error in spRbind(as(obj, SpatialPolygons), as(x, SpatialPolygons)) :
  non-unique polygon IDs

I have tried with the function rbind()
I have tried to set noumea$NOUMEA2_ID and dumbea$DUMBEA2_ with unique values 
but I still get the same message error as if the function spRbind was not using 
those two columns as polygons ID...

Does anyone know how I can fix this problem ?

Thanks a lot

Nicolas


Nicolas Degallier

IRD UMR182
Laboratoire d'Océanographie et du Climat, Expérimentation et Approches 
Numériques (LOCEAN)
Tour 45-55, 4e ét., case 100, 4 place Jussieu
75252  Paris Cedex 5  France
tél: (33) 01 44 27 51 57
fax: (33) 01 44 27 38 05

E-mail:nicolas.degall...@ird.fr
pdf reprints:
http://www.locean-ipsl.upmc.fr/~ndelod/production/

Skype:xuxaxu

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




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

[R-sig-Geo] Problem with Kriging with in strata

[Zia Ahmed]

 Hi,
I am trying to do kriging with in strata. I  use the example in book 
Applied Spatial Data Analysis with R page:217. It works fine with 
meuse data. But when I am  trying to use my data it shows some error 
massage.

Do I missing something here?
Thanks in advance
Zia

 Note: I am using following codes:


coordinates(tala) - ~ x + y
gridded(tala.grid) = ~x+y
 names(tala.grid)
[1] lt.a  lt.b  elev  FDID  wfe_was   
was   was.ltelev.sqrt

# lt.a= flooding land type 1 and lt.b=flooding land type 2

tala$lt.a - tala.grid$lt.a[overlay(tala.grid, tala)]
 names(tala)

[1] YID  FDID elev was  sas  lt.a lt.b


x1 - krige(sas ~ 1, tala[tala$lt.a == 0, ],tala.grid[tala.grid$lt.a == 
0, ], model = vgm(35,Sph, 2000, 15), nmin = 8, nmax = 40, maxdist = 1000)


Error in tala[tala$lt.a == 0, ] : NAs not permitted in row index Error 
in krige(sas ~ 1, tala[tala$lt.a == 0, ], tala.grid[tala.grid$lt.a ==  :
  error in evaluating the argument 'locations' in selecting a method 
for function 'krige'


x2- krige(sas~1, tala[tala$lt.a == 1, ],tala.grid[tala.grid$lt.a == 1, 
], model = vgm(30,Sph, 2000), nmin = 8, nmax = 40, maxdist = 1000)


Error in tala[tala$lt.a == 1, ] : NAs not permitted in row index Error 
in krige(sas ~ 1, tala[tala$lt.a == 1, ], tala.grid[tala.grid$lt.a ==  :
  error in evaluating the argument 'locations' in selecting a method 
for function 'krige'



On 1/22/2011 4:51 PM, Edzer Pebesma wrote:


On 01/22/2011 11:18 AM, Patrick Giraudoux wrote:

Hi,

Does anybody knows if it is possible to fix user defined maximum and
minimum values in the key of a xyplot (or spplot) graphics ? I have a
set of seperate drawings and I would like to keep the same color scheme
between fixed values for all of them (which looks like being not
possible, autokey normally adapting the range of colors to the
particular range of values in each plot.

If your object is of class SpatialPointsDataFrame, it may suffice to
pass the argument cuts with the vector of cuts. If it is
grid/polygons/lines, use the at component of the colorkey argument
to specify breaks; see ?levelplot in package lattice.


Patrick

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