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: 47 min 28 sec ago

Re: spatial autocorrelation in GAM residuals for large data set

Wed, 08/21/2019 - 10:15
On Wed, 21 Aug 2019, Elizabeth Webb wrote:

> Thank you, Roger for your help.   A quick follow-up:
>
> What do you mean when you say "Use one of the approaches described in
> the tutorial and you may be OK, but you should not trust the outcome of
> Moran's I on residuals without using an appropriate variant." ?  Or in
> other words, what is an appropriate variant in this context?

From Cliff & Ord (1973), we know that using the standard version of
Moran's I on OLS residuals is flawed, and they therefore proposed an
alternative taking into account the fact that the fitted values - needed
to calculate the residuals - depend on the fitted model, see the
differences that appear between running spdep::moran.test(...,
randomisation=FALSE) and lm.morantest(mod, ...) as one adds in covariates
in addition to the intercept. With just the intercept, they are identical,
as covariates are added, they diverge.

As of now, appropriate tests are not available for models other than OLS,
so one cannot take results as more than indicative. There has been work
moving from just lm() to glm() because the RHS may be seen as linear in
the covariates, but we don't have control of the distribution of the
residuals.

The problem is not widely discussed because it is intractable.

Had your pixels formed a rectangle, it might have been possible to use a
Moran eigenvector approach, as such approaches may use analytical
eigenvectors, but I am not aware of such work; proponents of Moran
eigenvectors claim that their use with larger data sets is possible; I
don't know what size data works in the spmoran package using ME.

Hope this helps,

Roger

>
> Elizabeth
>
> _______________________________________From: Roger Bivand <[hidden email]>
> Sent: Tuesday, August 20, 2019 4:43 PM
> To: Elizabeth Webb
> Cc: [hidden email]
> Subject: Re: [R-sig-Geo] spatial autocorrelation in GAM residuals for large data set
>
> On Tue, 20 Aug 2019, Elizabeth Webb wrote:
>
>> Hello,
>>
>> I have a large data set (~100k rows) containing observations at points
>> (MODIS pixels) across the northern hemisphere.  I have created a GAM
>> using the bam command in mgcv and I would like to check the model
>> residuals for spatial autocorrelation.
>>
>> One idea is to use the DHARMa package
>> (https://urldefense.proofpoint.com/v2/url?u=https-3A__cran.r-2Dproject.org_web_packages_DHARMa_vignettes_DHARMa.html-23spatial-2Dautocorrelation&d=DwIFAw&c=sJ6xIWYx-zLMB3EPkvcnVg&r=fIXSZTeOvV9o221vPRYLSw&m=IMILOSwwjZJkpsYiuj66ywrnluSI19Bn8ozD5p-NZks&s=lD9g96_TN9t-znQvV1M9V8CH2tgKAYHWKcTS8osjBSc&e= ).
>> The code looks something like this:
>>
>>     simulationOutput  <-   simulateResiduals(fittedModel = mymodel) # point at which R runs into memory problems
>>     testSpatialAutocorrelation(simulationOutput = simulationOutput, x =  data$latitude, y= data$longitude)
>>
>> However, this runs into memory problems.
>>
>> Another idea is to use the following code, after this tutorial
>> (https://urldefense.proofpoint.com/v2/url?u=http-3A__www.flutterbys.com.au_stats_tut_tut8.4a.html&d=DwIFAw&c=sJ6xIWYx-zLMB3EPkvcnVg&r=fIXSZTeOvV9o221vPRYLSw&m=IMILOSwwjZJkpsYiuj66ywrnluSI19Bn8ozD5p-NZks&s=v9Op69bwvOVRi1ujar_P8-LHsJFDAN_aE25i1_m22U4&e= ):
>>     library(ape)
>>     library(fields)
>>     coords = cbind(data$longitude, data$latitude)
>>    w = rdist(coords)  # point at which R runs into memory problems
>>     Moran.I(x = residuals(mymodel), w = w)
>>
>> But this also runs into memory problems.  I have tried increasing the
>> amount of memory allotted to R, but that just means R works for longer
>> before timing out.
>
> I do hope that you read
>
> https://urldefense.proofpoint.com/v2/url?u=https-3A__cran.r-2Dproject.org_web_packages_ape_vignettes_MoranI.pdf&d=DwIFAw&c=sJ6xIWYx-zLMB3EPkvcnVg&r=fIXSZTeOvV9o221vPRYLSw&m=IMILOSwwjZJkpsYiuj66ywrnluSI19Bn8ozD5p-NZks&s=tp6IoNz-8y-MBnLZldMpb3wUT_faHdoGMFszkZQxYBU&e=
>
> first, because the approach used in ape has been revised.
>
> The main problem is that ape uses by default a square matrix, and it is
> uncertain whether sparse matrices are accepted. This means that completely
> unneeded computations are carried out - dense matrices should never be
> used unless there is a convincing scientific argument (see
> https://urldefense.proofpoint.com/v2/url?u=https-3A__edzer.github.io_UseR2019_part2.html-23exercise-2Dreview-2D1&d=DwIFAw&c=sJ6xIWYx-zLMB3EPkvcnVg&r=fIXSZTeOvV9o221vPRYLSw&m=IMILOSwwjZJkpsYiuj66ywrnluSI19Bn8ozD5p-NZks&s=PIyJQgsz9qD81VCZsyfJQGdO-Gh6iJNgF2xH9jATbhI&e=  for a
> development on why distances are wasteful when edge counts on a graph do
> the same thing sparsely).
>
> Use one of the approaches described in the tutorial and you may be OK, but
> you should not trust the outcome of Moran's I on residuals without using
> an appropriate variant. Say you can represent your GAM with a linear model
> with say spline terms, you can use Moran's I for regression residuals.
> Take care that the average number of neighbours is very small (6-10), and
> large numbers of observations should not be a problem.
>
> A larger problem is that Moran's I (also for residuals) also responds to
> other mis-specifications than spatial autocorrelation, in particular
> missing variables and spatial processes with a different scale from the
> units of observation chosen.
>
>
>>
>> So, two questions: (1) Is there a memory efficient way to check for
>> spatial autocorrelation using Moran's I in large data sets? or (2) Is
>> there another way to check for spatial autocorrelation (besides Moran's
>> I) that won't have such memory problems?
>
> 1) Yes, see above, do not use dense matrices
>
> 2) Consider a higher level MRF term in your GAM for aggregates of your
> observations if such aggregation makes sense for your data.
>
> Hope this clarifies,
>
> Roger
>
>>
>> Thanks in advance,
>>
>> Elizabeth
>>
>>
>>
>>
>>
>>
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://urldefense.proofpoint.com/v2/url?u=https-3A__stat.ethz.ch_mailman_listinfo_r-2Dsig-2Dgeo&d=DwIFAw&c=sJ6xIWYx-zLMB3EPkvcnVg&r=fIXSZTeOvV9o221vPRYLSw&m=IMILOSwwjZJkpsYiuj66ywrnluSI19Bn8ozD5p-NZks&s=qq-Vo4xAxTCARWawQQYjCYvxm_Hg_iYYW1_n-Xphyg4&e=
>>
>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: [hidden email]
> https://urldefense.proofpoint.com/v2/url?u=https-3A__orcid.org_0000-2D0003-2D2392-2D6140&d=DwIFAw&c=sJ6xIWYx-zLMB3EPkvcnVg&r=fIXSZTeOvV9o221vPRYLSw&m=IMILOSwwjZJkpsYiuj66ywrnluSI19Bn8ozD5p-NZks&s=nyvB8TRse_NA-lALeTG3k_KOIVVzaNLMuDAqPo3dyGI&e=
> https://urldefense.proofpoint.com/v2/url?u=https-3A__scholar.google.no_citations-3Fuser-3DAWeghB0AAAAJ-26hl-3Den&d=DwIFAw&c=sJ6xIWYx-zLMB3EPkvcnVg&r=fIXSZTeOvV9o221vPRYLSw&m=IMILOSwwjZJkpsYiuj66ywrnluSI19Bn8ozD5p-NZks&s=UBVoMNMtGxwGxOGDcIJHGAuFm8gqb8X3kqNkGpPVKS4&e=
> --
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Parallel overlay of undefined multi-layer rasters

Wed, 08/21/2019 - 09:16
Dear list

I have an undefined number of multi-layer rasters (TIFF files with 10
bands), and I need to find the maximum (or median, etc.) among the rasters.
The calculations should operate band by band, so the expected output is one
stack also with 10 bands.

I thought of overlay, like this:


library(raster)
r <- raster(ncol=10, nrow=10)
r1 <- init(r, fun=runif)
r2 <- init(r, fun=runif)
r3 <- init(r, fun=runif)
r4 <- init(r, fun=runif)
s1<-stack(r1,r2)
s2<-stack(r3,r4)
s3<-stack(r1,r4)
f<-function(...){max(...,na.rm = TRUE)}

# do.call works for undefined stacks, but not in parallel
list1<-list(s1,s2,fun=f)
list2<-list(s1,s2,s3,fun=f)
do.call(overlay, list1)
do.call(overlay, list2)


As for parallel processing, I thought of clusteR, but it accepts only one
Raster* object.
The problem is that for running overlay with multi-layer rasters, they
cannot be stacked, because if stacked, the function will apply to all bands
and hence output a single band, while in this case the function should
apply band by band.

An alternative to this, is to split the bands of all images, process the
bands individually in parallel inside a loop (i.e. band 1 of all images,
band 2 of all images, etc.),
and then combine the outputs in a multi-layer raster. However, this is
inefficient, specially with large rasters not hold in memory.

Is there an elegant way to do this?

Thanks
Hugo

        [[alternative HTML version deleted]]

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

Re: Simulating variables with predefined correlation and autocorrelation

Wed, 08/21/2019 - 05:12
Dear Amitha,

I'm not sure I understand well your question (but then, I don't know
this Durbin model). From what I quickly looked up it follows this
formulation:

y = \rho W y + \alpha 1_n + X \beta + W X \theta + \varepsilon

where y are the observations (thus, what you are trying to simulate if I
understood correctly), W, X and 1_n are fixed and known matrices, and
the greek letters are the unknown parameters.

For simulating y from this model, you should assign values to the greek
letters and compute:

y = (I_n - \rho W)^{-1} [\alpha 1_n + X \beta + W X \theta + \varepsilon]

Does this help?

ƒacu.-


On 8/21/19 8:52 AM, Amitha Puranik wrote:
> Hello everyone,
>
> I'm trying to simulate a set of variables by specifying correlations. For
> the simulated data, I also want to specify autocorrelations. Basically I am
> trying to simulate data assuming spatial durbin model (a model which
> accounts for autocorrelation in both Y and X). I came across the post
> on ‘*Simulating
> spatially autocorrelated data*’ in R-sig-geo (
> https://stat.ethz.ch/pipermail/r-sig-geo/2011-September/012728.html) where
> similar query was discussed, however the focus there was on generating
> autocorrelated dependent variable.
>
> Can anyone assist me on this? Thanks in advance.
>
> [[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: Rasterize point process using specific rule

Wed, 08/21/2019 - 04:06
If I understand your question correctly you can simply use the spatstat
function `discs()`. Given `syn.ppp` created from your code the call is
something like:

     D <- discs(syn.ppp, radii = marks(syn.ppp)/2, mask = TRUE, eps = 1)

For your dataset this covers almost the entire window if both the
coordinates and the marks are in meters as I assume, but that is up to
you to investigate.

/Ege

On 20/08/2019 15.12, ASANTOS via R-sig-Geo wrote:
> Dear members,
>
> I've like to create a raster using my event coordinates and information
> about the size in this coordinates as a rule for my raster creation.
> First, I make:
>
> library(raster)
> library(spatstat)
>
> #Points coordinates in UTM
> xp<-c(371278.588,371250.722,371272.618,371328.421,371349.974,
> 371311.95,371296.265,371406.46,371411.551,371329.041,371338.081,
> 371334.182,371333.756,371299.818,371254.374,371193.673,371172.836,
> 371173.803,371153.73,371165.051,371140.417,371168.279,371166.367,
> 371180.575,371132.664,371129.791,371232.919,371208.502,371289.462,
> 371207.595,371219.008,371139.921,371133.215,371061.467,371053.69,
> 371099.897,371108.782,371112.52,371114.241,371176.236,371159.185,
> 371159.291,371158.552,370978.252,371120.03,371116.993)
>
> yp<-c(8246507.94,8246493.176,8246465.974,8246464.413,8246403.465,
> 8246386.098,8246432.144,8246394.827,8246366.201,8246337.626,
> 8246311.125,8246300.039,8246299.594,8246298.072,8246379.351,
> 8246431.998,8246423.913,8246423.476,8246431.658,8246418.226,
> 8246400.161,8246396.891,8246394.225,8246400.391,8246370.244,
> 8246367.019,8246311.075,8246255.174,8246255.085,8246226.514,
> 8246215.847,8246337.316,8246330.197,8246311.197,8246304.183,
> 8246239.282,8246239.887,8246241.678,8246240.361,8246167.364,
> 8246171.581,8246171.803,8246169.807,8246293.57,8246183.194,8246189.926)
>
> #Now I have the size of each nest in meters (marked process)
>
> area<-c(117,30,4,341,15,160,35,280,108,168,63,143,2,48,182,42,
> 88,56,27,156,288,45,49,234,72,270,91,40,304,56,35,4,56.7,9,4.6,
> 105,133,135,23.92,190,12.9,15.2,192.78,104,255,24)
>
> # Make a contour for the window creation
> W <- convexhull.xy(xp,yp)
>
> #Create ppm object
> syn.ppp <- ppp(x=xp,y=yp,window=W, marks=area)
>
> # Plot the point process
> plot(syn.ppp)
>
> #Definition of raster resolution - 1 meter square
> ext <- as(extent(c(W$xrange,W$yrange)), "SpatialPolygons")
> syn.ras.res<-rasterToPoints(raster(ext, resolution = 1), spatial = TRUE)
>
> #Now Rasterize
> syn.ras<- rasterize(syn.ras.res@coords, raster(syn.ras.res), *?????*)
>
> I' ve like to create some kind of rule in *?????*, where pixel values in
> my final raster was 1 when inside de area (marks=area) and zero for
> outside de circular area given by points coordinates
> (syn.ras.res@coords) using 1 meter as resolution?
>
> Any ideas, please
>
> Thanks,
>
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Simulating variables with predefined correlation and autocorrelation

Wed, 08/21/2019 - 01:52
Hello everyone,

I'm trying to simulate a set of variables by specifying correlations. For
the simulated data, I also want to specify autocorrelations. Basically I am
trying to simulate data assuming spatial durbin model (a model which
accounts for autocorrelation in both Y and X). I came across the post
on ‘*Simulating
spatially autocorrelated data*’ in R-sig-geo (
https://stat.ethz.ch/pipermail/r-sig-geo/2011-September/012728.html) where
similar query was discussed, however the focus there was on generating
autocorrelated dependent variable.

Can anyone assist me on this? Thanks in advance.

        [[alternative HTML version deleted]]

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

Re: spatial autocorrelation in GAM residuals for large data set

Tue, 08/20/2019 - 20:21
Thank you, Roger for your help.   A quick follow-up:

What do you mean when you say "Use one of the approaches described in the tutorial and you may be OK, but you should not trust the outcome of Moran's I on residuals without using an appropriate variant." ?  Or in other words, what is an appropriate variant in this context?

Elizabeth

_______________________________________From: Roger Bivand <[hidden email]>
Sent: Tuesday, August 20, 2019 4:43 PM
To: Elizabeth Webb
Cc: [hidden email]
Subject: Re: [R-sig-Geo] spatial autocorrelation in GAM residuals for large data set

On Tue, 20 Aug 2019, Elizabeth Webb wrote:

> Hello,
>
> I have a large data set (~100k rows) containing observations at points
> (MODIS pixels) across the northern hemisphere.  I have created a GAM
> using the bam command in mgcv and I would like to check the model
> residuals for spatial autocorrelation.
>
> One idea is to use the DHARMa package
> (https://urldefense.proofpoint.com/v2/url?u=https-3A__cran.r-2Dproject.org_web_packages_DHARMa_vignettes_DHARMa.html-23spatial-2Dautocorrelation&d=DwIFAw&c=sJ6xIWYx-zLMB3EPkvcnVg&r=fIXSZTeOvV9o221vPRYLSw&m=IMILOSwwjZJkpsYiuj66ywrnluSI19Bn8ozD5p-NZks&s=lD9g96_TN9t-znQvV1M9V8CH2tgKAYHWKcTS8osjBSc&e= ).
> The code looks something like this:
>
>     simulationOutput  <-   simulateResiduals(fittedModel = mymodel) # point at which R runs into memory problems
>     testSpatialAutocorrelation(simulationOutput = simulationOutput, x =  data$latitude, y= data$longitude)
>
> However, this runs into memory problems.
>
> Another idea is to use the following code, after this tutorial
> (https://urldefense.proofpoint.com/v2/url?u=http-3A__www.flutterbys.com.au_stats_tut_tut8.4a.html&d=DwIFAw&c=sJ6xIWYx-zLMB3EPkvcnVg&r=fIXSZTeOvV9o221vPRYLSw&m=IMILOSwwjZJkpsYiuj66ywrnluSI19Bn8ozD5p-NZks&s=v9Op69bwvOVRi1ujar_P8-LHsJFDAN_aE25i1_m22U4&e= ):
>     library(ape)
>     library(fields)
>     coords = cbind(data$longitude, data$latitude)
>    w = rdist(coords)  # point at which R runs into memory problems
>     Moran.I(x = residuals(mymodel), w = w)
>
> But this also runs into memory problems.  I have tried increasing the
> amount of memory allotted to R, but that just means R works for longer
> before timing out.
I do hope that you read

https://urldefense.proofpoint.com/v2/url?u=https-3A__cran.r-2Dproject.org_web_packages_ape_vignettes_MoranI.pdf&d=DwIFAw&c=sJ6xIWYx-zLMB3EPkvcnVg&r=fIXSZTeOvV9o221vPRYLSw&m=IMILOSwwjZJkpsYiuj66ywrnluSI19Bn8ozD5p-NZks&s=tp6IoNz-8y-MBnLZldMpb3wUT_faHdoGMFszkZQxYBU&e=

first, because the approach used in ape has been revised.

The main problem is that ape uses by default a square matrix, and it is
uncertain whether sparse matrices are accepted. This means that completely
unneeded computations are carried out - dense matrices should never be
used unless there is a convincing scientific argument (see
https://urldefense.proofpoint.com/v2/url?u=https-3A__edzer.github.io_UseR2019_part2.html-23exercise-2Dreview-2D1&d=DwIFAw&c=sJ6xIWYx-zLMB3EPkvcnVg&r=fIXSZTeOvV9o221vPRYLSw&m=IMILOSwwjZJkpsYiuj66ywrnluSI19Bn8ozD5p-NZks&s=PIyJQgsz9qD81VCZsyfJQGdO-Gh6iJNgF2xH9jATbhI&e=  for a
development on why distances are wasteful when edge counts on a graph do
the same thing sparsely).

Use one of the approaches described in the tutorial and you may be OK, but
you should not trust the outcome of Moran's I on residuals without using
an appropriate variant. Say you can represent your GAM with a linear model
with say spline terms, you can use Moran's I for regression residuals.
Take care that the average number of neighbours is very small (6-10), and
large numbers of observations should not be a problem.

A larger problem is that Moran's I (also for residuals) also responds to
other mis-specifications than spatial autocorrelation, in particular
missing variables and spatial processes with a different scale from the
units of observation chosen.


>
> So, two questions: (1) Is there a memory efficient way to check for
> spatial autocorrelation using Moran's I in large data sets? or (2) Is
> there another way to check for spatial autocorrelation (besides Moran's
> I) that won't have such memory problems?

1) Yes, see above, do not use dense matrices

2) Consider a higher level MRF term in your GAM for aggregates of your
observations if such aggregation makes sense for your data.

Hope this clarifies,

Roger

>
> Thanks in advance,
>
> Elizabeth
>
>
>
>
>
>
>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://urldefense.proofpoint.com/v2/url?u=https-3A__stat.ethz.ch_mailman_listinfo_r-2Dsig-2Dgeo&d=DwIFAw&c=sJ6xIWYx-zLMB3EPkvcnVg&r=fIXSZTeOvV9o221vPRYLSw&m=IMILOSwwjZJkpsYiuj66ywrnluSI19Bn8ozD5p-NZks&s=qq-Vo4xAxTCARWawQQYjCYvxm_Hg_iYYW1_n-Xphyg4&e=
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://urldefense.proofpoint.com/v2/url?u=https-3A__orcid.org_0000-2D0003-2D2392-2D6140&d=DwIFAw&c=sJ6xIWYx-zLMB3EPkvcnVg&r=fIXSZTeOvV9o221vPRYLSw&m=IMILOSwwjZJkpsYiuj66ywrnluSI19Bn8ozD5p-NZks&s=nyvB8TRse_NA-lALeTG3k_KOIVVzaNLMuDAqPo3dyGI&e=
https://urldefense.proofpoint.com/v2/url?u=https-3A__scholar.google.no_citations-3Fuser-3DAWeghB0AAAAJ-26hl-3Den&d=DwIFAw&c=sJ6xIWYx-zLMB3EPkvcnVg&r=fIXSZTeOvV9o221vPRYLSw&m=IMILOSwwjZJkpsYiuj66ywrnluSI19Bn8ozD5p-NZks&s=UBVoMNMtGxwGxOGDcIJHGAuFm8gqb8X3kqNkGpPVKS4&e=

        [[alternative HTML version deleted]]

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

Re: rgdal::writeOGR with driver='ESRI Shapefile' converts Polygon object into a hole

Tue, 08/20/2019 - 18:06
Hi Phil and Roger,

To close the loop on this issue, I have confirmed that the mapshaper javascript library (which rmapshaper wraps) can produce invalid polygons (https://github.com/ateucher/rmapshaper/issues/89), and I will address this in a future release of rmapshaper (https://github.com/ateucher/rmapshaper/issues/90). I have also alerted the maintainer of mapshaper and he is going to add functionality to address this as well (https://github.com/mbloch/mapshaper/issues/352).

Thanks for the thorough investigation and bringing this to my attention!

Cheers,
Andy Teucher (rmapshaper package maintainer)

> Date: Tue, 20 Aug 2019 11:55:00 +0200
> From: Roger Bivand <[hidden email]>
> To: Phil Haines <[hidden email]>
> Cc: <[hidden email]>
> Subject: Re: [R-sig-Geo]  rgdal::writeOGR with driver='ESRI Shapefile'
> converts Polygon object into a hole
> Message-ID: <[hidden email]>
> Content-Type: text/plain; charset="us-ascii"; Format="flowed"
>
> Good, thanks very much!
>
> Roger
>
> On Mon, 19 Aug 2019, Phil Haines wrote:
>
>> Hi Roger,
>>
>> I have added an example to https://github.com/ateucher/rmapshaper/issues/89
>>
>> Hopefully it illustrates the behaviour I was seeing from
>> rmapshaper::ms_simplify() but please let me know if not.
>>
>> And thanks for the nudge towards GPKG, I will investigate.
>> Phil
>>
>> On Mon, 19 Aug 2019 at 15:48, Roger Bivand <[hidden email]> wrote:
>>>
>>> Hi Phil,
>>>
>>> On Sun, 18 Aug 2019, Phil Haines wrote:
>>>
>>>> Hi Roger,
>>>>
>>>> I originally encountered this issue while trying to reduce the size of
>>>> German postal code boundaries. I used rmapshaper::ms_simplify() which
>>>> introduced the polygons sharing a common edge. I definitely didn't
>>>> need them to be separate polygons, and would happily have merged them
>>>> had I been more familiar with gUnaryUnion().
>>>
>>> If the German postal code boundaries are open, could you please put (an
>>> affected subset) on an URL, and provide a short script to replicate the
>>> workflow. Then we can engage with the rmapshaper maintainer, and see
>>> whether the underlying problem in the workflow is in the js library or its
>>> R wrapper. I've opened: https://github.com/ateucher/rmapshaper/issues/89.
>>>
>>>>
>>>> I apologise if my question has resulted in unnecessary effort on your
>>>> part. I reported the behaviour because I found it surprising. However,
>>>> I now understand that this example is not a valid simple feature
>>>> geometry and that the solution is to take steps to ensure I have valid
>>>> geometries, and to repair if not.
>>>>
>>>
>>> No need to apologise!! Your reproducible example was excellent, without it
>>> you wouldn't have contributed to getting the internal shapelib code
>>> changed to make it less restrictive in future releases of GDAL,
>>> potentially benefitting lots of users (who really should drop ESRI
>>> shapefiles, and/or should check geometries for validity, but ... whose
>>> work you've helped save). By the way, ESRI would be very happy if
>>> shapefiles stopped being produced in current workflows, and that new files
>>> should rather be GPKG.
>>>
>>>> Additionally I must confess that I only use the ESRI Shapefile driver
>>>> because I don't know any better... My use cases are reading and
>>>> writing from R (usually to produce maps) and occasionally opening in
>>>> QGIS. I can happily switch to any format that supports this workflow.
>>>>
>>>
>>> GPKG is now very broadly supported, is SF-compliant, and has the big
>>> user-facing benefits of not having field name length restrictions, and
>>> many fewer encoding issues for field names and string values.
>>>
>>> Best wishes,
>>>
>>> Roger
>>>
>>>> Thank you very much for your time,
>>>> Phil
>>>>
>>>> On Sun, 18 Aug 2019 at 13:59, Roger Bivand <[hidden email]> wrote:
>>>>>
>>>>> The reason that the problem occurred is that a MULTIPOLYGON with two
>>>>> exterior rings becomes invalid if the exterior rings touch along an edge
>>>>> (this case). It is important to know the use case, to see whether:
>>>>>
>>>>> library(rgeos)
>>>>> writeWKT(Ps1)
>>>>> gIsValid(Ps1, reason=TRUE)
>>>>> Ps1a <- gUnaryUnion(gBuffer(Ps1, width=0))
>>>>> gIsValid(Ps1a, reason=TRUE)
>>>>> writeWKT(Ps1a)
>>>>>
>>>>> or equivalently in sf for sf objects, should be applied before trying to
>>>>> write the object out to file with this driver.
>>>>>
>>>>> However, because drivers that are compliant with the simple features
>>>>> standard (which bans exterior rings sharing edges) have been permissive
>>>>> and do round-trip this invalid object, a relaxation in the OGR ESRI
>>>>> shapefile driver has been provided and may be included in a future
>>>>> release.
>>>>>
>>>>> We need to know (see these issues):
>>>>>
>>>>> https://github.com/r-spatial/sf/issues/1130
>>>>> https://github.com/OSGeo/gdal/issues/1787
>>>>>
>>>>> why it was desirable to write out this object using this driver? Could an
>>>>> alternative driver have been used, or is ESRI shapefile the only format
>>>>> used in the workflow?
>>>>>
>>>>> If it has to be this driver, could the workflow be changed to repair
>>>>> degenerate cases before writing? If using sp classes, rgeos may be used to
>>>>> test for and probably repair such geometries before they reach
>>>>> rgdal::writeOGR() for this driver. Adding code to sf and rgdal to trap
>>>>> degenerate cases does encumber all users with valid geometries with
>>>>> the time wasted on extra checking, so building checks into sf and rgdal is
>>>>> not desirable.
>>>>>
>>>>> I hope it is possible to find out more about the use case quickly, to pass
>>>>> on to GDAL developers to help motivate a relaxation in their current
>>>>> policy with regard to this driver, and to encourage them to include the
>>>>> fix branch in a future release of GDAL.
>>>>>
>>>>> Roger
>>>>>
>>>>> On Sat, 17 Aug 2019, Roger Bivand wrote:
>>>>>
>>>>>> Please follow up both here and on:
>>>>>>
>>>>>> https://github.com/r-spatial/sf/issues/1130
>>>>>>
>>>>>> as the problem is also seen in the sf package using the same GDAL ESRI
>>>>>> Shapefile driver.
>>>>>>
>>>>>> Roger
>>>>>>
>>>>>> On Fri, 16 Aug 2019, Roger Bivand wrote:
>>>>>>
>>>>>>> On Fri, 16 Aug 2019, Roger Bivand wrote:
>>>>>>>
>>>>>>>>  On Tue, 13 Aug 2019, Phil Haines wrote:
>>>>>>>>
>>>>>>>>>   Dear list,
>>>>>>>>>
>>>>>>>>>   I have a single Polygons object containing multiple Polygon objects
>>>>>>>>>   that share a common border. When I output this using writeOGR() one of
>>>>>>>>>   the Polygon objects becomes a hole, as the following example shows.
>>>>>>>>>
>>>>>>>>>   Create a Polygons object containing two adjoining Polygon objects
>>>>>>>>>>   library(rgdal)
>>>>>>>>>>   r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
>>>>>>>>>>   r2 <- r1; r2[,1] <- r2[,1]+1
>>>>>>>>>>   Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
>>>>>>>>>>   SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
>>>>>>>>>>   data.frame(Example=c("Minimal")))
>>>>>>>>>
>>>>>>>>>   Perform a write/readOGR() cycle
>>>>>>>>>>   fn <- tempfile()
>>>>>>>>>>   writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
>>>>>>>>>>   SPDF2 <- readOGR(dsn=fn,layer='test')
>>>>>>>>>
>>>>>>>>>   Second Polygon object is now a hole
>>>>>>>>>>   sapply(SPDF2@polygons[[1]]@Polygons,slot,"hole")
>>>>>>>>>   [1] FALSE  TRUE
>>>>>>>>>
>>>>>>>>>   I see from the sp documentation that "Polygon objects belonging to a
>>>>>>>>>   Polygons object should either not overlap one-other, or should be
>>>>>>>>>   fully included" but I am not sure how this relates to bordering
>>>>>>>>>   Polygon objects. I would welcome any advice as to whether what I am
>>>>>>>>>   asking of writeOGR is reasonable?
>>>>>>>>
>>>>>>>>  The problem is with the 'ESRI Shapefile' representation and driver:
>>>>>>>>
>>>>>>>>  library(rgdal)
>>>>>>>>  r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
>>>>>>>>  r2 <- r1; r2[,1] <- r2[,1]+1
>>>>>>>>  Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
>>>>>>>>  SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
>>>>>>>>  data.frame(Example=c("Minimal")))
>>>>>>>>  sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "hole")
>>>>>>>>  sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>>>>>>>
>>>>>>>>  # which constructs a MULTIPOLYGON object:
>>>>>>>>
>>>>>>>>  rgeos::writeWKT(SPDF)
>>>>>>>>  library(sf)
>>>>>>>>  st_as_text(st_geometry(st_as_sf(SPDF)))
>>>>>>>>
>>>>>>>> #    The 'ESRI Shapefile' driver is not Simple-Feature compliant (it
>>>>>>>> #    pre-dates it), so the failure occurs by seeing the second exterior
>>>>>>>> #    ring
>>>>>>>> #    as an interior ring
>>>>>>>>
>>>>>>>>  fn <- tempfile()
>>>>>>>>  writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
>>>>>>>>  SPDF2 <- readOGR(dsn=fn, layer='test')
>>>>>>>>  sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "hole")
>>>>>>>>  sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>>>>>>>  rgeos::writeWKT(SPDF2)
>>>>>>>>  st_as_text(st_geometry(st_as_sf(SPDF2)))
>>>>>>>>
>>>>>>>>  # This happens with sf too, using the same GDAL driver:
>>>>>>>>
>>>>>>>>  sf2 <- st_read(dsn=fn, layer='test')
>>>>>>>>  st_geometry(sf2)
>>>>>>>>  library(sf)
>>>>>>>>  st_as_text(st_geometry(st_as_sf(sf2)))
>>>>>>>>  rgeos::writeWKT(as(sf2, "Spatial"))
>>>>>>>>
>>>>>>>>  # Adding the comment fix doesn't help:
>>>>>>>>
>>>>>>>>  comment(slot(SPDF, "polygons")[[1]])
>>>>>>>>  SPDF_c <- rgeos::createSPComment(SPDF)
>>>>>>>>  comment(slot(SPDF_c, "polygons")[[1]])
>>>>>>>>  writeOGR(SPDF_c, fn, layer='test_c', driver='ESRI Shapefile',
>>>>>>>>    verbose=TRUE)
>>>>>>>>
>>>>>>>> #    reports
>>>>>>>> #    Object initially classed as: wkbPolygon
>>>>>>>> #    SFS comments in Polygons objects
>>>>>>>> #    Object reclassed as: wkbMultiPolygon
>>>>>>>>
>>>>>>>>  SPDF2_c <- readOGR(dsn=fn, layer='test_c')
>>>>>>>>  sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot, "hole")
>>>>>>>>  sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot,
>>>>>>>>  "ringDir")
>>>>>>>>  rgeos::writeWKT(SPDF2_c)
>>>>>>>>  st_as_text(st_geometry(st_as_sf(SPDF2_c)))
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>  # If the input object is written out with the GeoPackage driver:
>>>>>>>>
>>>>>>>>  fn1 <- tempfile(fileext=".gpkg")
>>>>>>>>  writeOGR(SPDF, fn1, layer="test", driver='GPKG')
>>>>>>>>  sf2a <- st_read(dsn=fn1,layer='test')
>>>>>>>>  st_coordinates(st_geometry(sf2a))
>>>>>>>>  SPDF2a <- readOGR(dsn=fn1)
>>>>>>>>  sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "hole")
>>>>>>>>  sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>>>>>>>  rgeos::writeWKT(SPDF2a)
>>>>>>>>  st_as_text(st_geometry(st_as_sf(SPDF2a)))
>>>>>>>>
>>>>>>>>  # the issue is resolved. If we separate the exterior rings:
>>>>>>>>
>>>>>>>>  r2a <- r1; r2a[,1] <- r2a[,1]+1.00001
>>>>>>>>  Ps1a = Polygons(list(Polygon(r1),Polygon(r2a)),ID=1)
>>>>>>>>  SPDFa = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1a)),
>>>>>>>>  data.frame(Example=c("Minimal")))
>>>>>>>>  fna <- tempfile()
>>>>>>>>  writeOGR(SPDFa, fna, layer='test', driver='ESRI Shapefile')
>>>>>>>>  SPDF2_a <- readOGR(dsn=fna, layer='test')
>>>>>>>>  sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot, "hole")
>>>>>>>>  sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot,
>>>>>>>>  "ringDir")
>>>>>>>>  rgeos::writeWKT(SPDF2_a)
>>>>>>>>  st_as_text(st_geometry(st_as_sf(SPDF2_a)))
>>>>>>>>
>>>>>>>>  # we are OK as the two exterior rings do not touch.
>>>>>>>>
>>>>>>>>  # does using sf make a difference?
>>>>>>>>
>>>>>>>>  fn_s <- tempfile(fileext=".shp")
>>>>>>>>  st_write(st_as_sf(SPDF), dsn=fn_s)
>>>>>>>>  sf_in <- st_read(fn_s)
>>>>>>>>  st_as_text(st_geometry(st_as_sf(sf_in)))
>>>>>>>>
>>>>>>>>  # No
>>>>>>>>
>>>>>>>>  fn_s <- tempfile(fileext=".shp")
>>>>>>>>  st_write(st_as_sf(SPDF_c), dsn=fn_s)
>>>>>>>>  sf_in_c <- st_read(fn_s)
>>>>>>>>  st_as_text(st_geometry(st_as_sf(sf_in_c)))
>>>>>>>>
>>>>>>>>  # nor with the pretend-SF-compliant comment set either.
>>>>>>>>
>>>>>>>>  So the weakness is in the "ESRI Shapefile" write driver, or possibly in
>>>>>>>>  the OGRGeometryFactory::organizePolygons() function in GDAL used in
>>>>>>>>  OGR_write() (a C++ function) called by writeOGR(). If sf::st_write()
>>>>>>>>  also
>>>>>>>>  calls OGRGeometryFactory::organizePolygons(), we'd maybe consider that
>>>>>>>>  it
>>>>>>>>  has a weakness for the "ESRI Shapefile" driver, but which does not
>>>>>>>>  affect
>>>>>>>>  SF-compliant drivers.
>>>>>>>
>>>>>>> Without the comment set, OGRGeometryFactory::organizePolygons() is used;
>>>>>>> with it set, OGRGeometryFactory::organizePolygons() is not used, because
>>>>>>> the object is declared as two exterior rings. In both cases, we have the
>>>>>>> output object written out and read back in incorrectly with the ESRI
>>>>>>> shapefile driver, but SF-compliant drivers round-trip (in the test
>>>>>>> GeoJSON), correctly.
>>>>>>>
>>>>>>> It is likely that the changes made in 2015 to accommodate GeoJSON led to
>>>>>>> this possible regression for the ESRI Shapefile driver. I'm adding this
>>>>>>> geometry to tests/tripup.R and data files; without code changes the hole
>>>>>>> slot is wrong and the ring direction is changes to match, so the
>>>>>>> coordinates change too.
>>>>>>>
>>>>>>> Reading using the deprecated maptools::readShapeSpatial() also gets a hole
>>>>>>> in the second external ring. However, writing with deprecated
>>>>>>> maptools::  writeSpatialShape() yields a shapefile that when read with
>>>>>>> maptools::  readShapeSpatial() gives the correct two exterior ring, no hole
>>>>>>> object. When read with sf::st_read() and rgdal::readOGR(), the object is
>>>>>>> also correct. So the problem definitely lies in rgdal::writeOGR(), and
>>>>>>> sf::st_write() - roundtripping with sf::st_write() and sf::st_read()
>>>>>>> degrades from MULTIPOLYGON to POLYGON with the ESRI Shapefile driver.
>>>>>>>
>>>>>>> Roger
>>>>>>>
>>>>>>>>
>>>>>>>>  This is probably related to a similar but inverse problem with the
>>>>>>>>  SF-compliant GeoJSON driver in 2015:
>>>>>>>>
>>>>>>>>  https://stat.ethz.ch/pipermail/r-sig-geo/2015-October/023609.html
>>>>>>>>
>>>>>>>>  continued the next month in:
>>>>>>>>
>>>>>>>>  https://stat.ethz.ch/pipermail/r-sig-geo/2015-November/023656.html
>>>>>>>>
>>>>>>>>  The details are in this SVN diff
>>>>>>>>
>>>>>>>>  https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=555&r2=571
>>>>>>>>
>>>>>>>>  up to the end og the list thread, and this one from then until now:
>>>>>>>>
>>>>>>>>  https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=571&r2=733
>>>>>>>>
>>>>>>>>  Summary: could you change drivers, or is it really necessary to fix an
>>>>>>>>  EOL
>>>>>>>>  problem? What is your use case?
>>>>>>>>
>>>>>>>>>
>>>>>>>>>   Thanks in advance for your time,
>>>>>>>>
>>>>>>>>  Thanks for a complete example,
>>>>>>>>
>>>>>>>>  Roger
>>>>>>>>
>>>>>>>>>   Phil
>>>>>>>>>
>>>>>>>>>>   sessionInfo()
>>>>>>>>>   R version 3.5.1 (2018-07-02)
>>>>>>>>>   Platform: x86_64-w64-mingw32/x64 (64-bit)
>>>>>>>>>   Running under: Windows >= 8 x64 (build 9200)
>>>>>>>>>
>>>>>>>>>   Matrix products: default
>>>>>>>>>
>>>>>>>>>   locale:
>>>>>>>>>   [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United
>>>>>>>>>   Kingdom.1252    LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
>>>>>>>>>                           LC_TIME=English_United Kingdom.1252
>>>>>>>>>
>>>>>>>>>   attached base packages:
>>>>>>>>>   [1] stats     graphics  grDevices utils     datasets  methods   base
>>>>>>>>>
>>>>>>>>>   other attached packages:
>>>>>>>>>   [1] rgdal_1.4-4 sp_1.3-1
>>>>>>>>>
>>>>>>>>>   loaded via a namespace (and not attached):
>>>>>>>>>   [1] compiler_3.5.1  tools_3.5.1     yaml_2.2.0      grid_3.5.1
>>>>>>>>>   lattice_0.20-35
>>>>>>>>>
>>>>>>>>>   _______________________________________________
>>>>>>>>>   R-sig-Geo mailing list
>>>>>>>>>   [hidden email]
>>>>>>>>>   https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>> --
>>>>> Roger Bivand
>>>>> Department of Economics, Norwegian School of Economics,
>>>>> Helleveien 30, N-5045 Bergen, Norway.
>>>>> voice: +47 55 95 93 55; e-mail: [hidden email]
>>>>> https://orcid.org/0000-0003-2392-6140
>>>>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>>>
>>>
>>> --
>>> Roger Bivand
>>> Department of Economics, Norwegian School of Economics,
>>> Helleveien 30, N-5045 Bergen, Norway.
>>> voice: +47 55 95 93 55; e-mail: [hidden email]
>>> https://orcid.org/0000-0003-2392-6140
>>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>
>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: [hidden email]
> https://orcid.org/0000-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: spatial autocorrelation in GAM residuals for large data set

Tue, 08/20/2019 - 15:43
On Tue, 20 Aug 2019, Elizabeth Webb wrote:

> Hello,
>
> I have a large data set (~100k rows) containing observations at points
> (MODIS pixels) across the northern hemisphere.  I have created a GAM
> using the bam command in mgcv and I would like to check the model
> residuals for spatial autocorrelation.
>
> One idea is to use the DHARMa package
> (https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html#spatial-autocorrelation).
> The code looks something like this:
>
>     simulationOutput  <-   simulateResiduals(fittedModel = mymodel) # point at which R runs into memory problems
>     testSpatialAutocorrelation(simulationOutput = simulationOutput, x =  data$latitude, y= data$longitude)
>
> However, this runs into memory problems.
>
> Another idea is to use the following code, after this tutorial
> (http://www.flutterbys.com.au/stats/tut/tut8.4a.html):
>     library(ape)
>     library(fields)
>     coords = cbind(data$longitude, data$latitude)     
>    w = rdist(coords)  # point at which R runs into memory problems
>     Moran.I(x = residuals(mymodel), w = w)
>
> But this also runs into memory problems.  I have tried increasing the
> amount of memory allotted to R, but that just means R works for longer
> before timing out. I do hope that you read

https://cran.r-project.org/web/packages/ape/vignettes/MoranI.pdf

first, because the approach used in ape has been revised.

The main problem is that ape uses by default a square matrix, and it is
uncertain whether sparse matrices are accepted. This means that completely
unneeded computations are carried out - dense matrices should never be
used unless there is a convincing scientific argument (see
https://edzer.github.io/UseR2019/part2.html#exercise-review-1 for a
development on why distances are wasteful when edge counts on a graph do
the same thing sparsely).

Use one of the approaches described in the tutorial and you may be OK, but
you should not trust the outcome of Moran's I on residuals without using
an appropriate variant. Say you can represent your GAM with a linear model
with say spline terms, you can use Moran's I for regression residuals.
Take care that the average number of neighbours is very small (6-10), and
large numbers of observations should not be a problem.

A larger problem is that Moran's I (also for residuals) also responds to
other mis-specifications than spatial autocorrelation, in particular
missing variables and spatial processes with a different scale from the
units of observation chosen.


>
> So, two questions: (1) Is there a memory efficient way to check for
> spatial autocorrelation using Moran's I in large data sets? or (2) Is
> there another way to check for spatial autocorrelation (besides Moran's
> I) that won't have such memory problems?

1) Yes, see above, do not use dense matrices

2) Consider a higher level MRF term in your GAM for aggregates of your
observations if such aggregation makes sense for your data.

Hope this clarifies,

Roger

>
> Thanks in advance,
>
> Elizabeth
>
>
>
>
>
>
>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Re: Rasterize point process using specific rule

Tue, 08/20/2019 - 10:26
Dear Alexandre,
I'm not familiar with point processes, but as far as I understand the
task well, it can be easily solved using sp+raster or sf+raster.

library(sf)
library(raster)
points_sf <- st_as_sf(data.frame(xp, yp, area), coords = c("xp", "yp"),
crs = 32629)
circles_sf <- st_buffer(x = points, dist = sqrt(points_sf$area / pi)) #
if it is really the area, or dist = points_sf$area if it is the radius
empty_raster <- raster(points_sf, resolution = 1)
filled_raster <- rasterize(x = circles_sf, y = empty_raster, field = 1,
fun = function(x, na.rm) 1, background = 0)
plot(filled_raster)

HTH,
Ákos Bede-Fazekas
Hungarian Academy of Sciences


2019.08.20. 15:12 keltezéssel, ASANTOS via R-sig-Geo írta:
> Dear members,
>
> I've like to create a raster using my event coordinates and information
> about the size in this coordinates as a rule for my raster creation.
> First, I make:
>
> library(raster)
> library(spatstat)
>
> #Points coordinates in UTM
> xp<-c(371278.588,371250.722,371272.618,371328.421,371349.974,
> 371311.95,371296.265,371406.46,371411.551,371329.041,371338.081,
> 371334.182,371333.756,371299.818,371254.374,371193.673,371172.836,
> 371173.803,371153.73,371165.051,371140.417,371168.279,371166.367,
> 371180.575,371132.664,371129.791,371232.919,371208.502,371289.462,
> 371207.595,371219.008,371139.921,371133.215,371061.467,371053.69,
> 371099.897,371108.782,371112.52,371114.241,371176.236,371159.185,
> 371159.291,371158.552,370978.252,371120.03,371116.993)
>
> yp<-c(8246507.94,8246493.176,8246465.974,8246464.413,8246403.465,
> 8246386.098,8246432.144,8246394.827,8246366.201,8246337.626,
> 8246311.125,8246300.039,8246299.594,8246298.072,8246379.351,
> 8246431.998,8246423.913,8246423.476,8246431.658,8246418.226,
> 8246400.161,8246396.891,8246394.225,8246400.391,8246370.244,
> 8246367.019,8246311.075,8246255.174,8246255.085,8246226.514,
> 8246215.847,8246337.316,8246330.197,8246311.197,8246304.183,
> 8246239.282,8246239.887,8246241.678,8246240.361,8246167.364,
> 8246171.581,8246171.803,8246169.807,8246293.57,8246183.194,8246189.926)
>
> #Now I have the size of each nest in meters (marked process)
>
> area<-c(117,30,4,341,15,160,35,280,108,168,63,143,2,48,182,42,
> 88,56,27,156,288,45,49,234,72,270,91,40,304,56,35,4,56.7,9,4.6,
> 105,133,135,23.92,190,12.9,15.2,192.78,104,255,24)
>
> # Make a contour for the window creation
> W <- convexhull.xy(xp,yp)
>
> #Create ppm object
> syn.ppp <- ppp(x=xp,y=yp,window=W, marks=area)
>
> # Plot the point process
> plot(syn.ppp)
>
> #Definition of raster resolution - 1 meter square
> ext <- as(extent(c(W$xrange,W$yrange)), "SpatialPolygons")
> syn.ras.res<-rasterToPoints(raster(ext, resolution = 1), spatial = TRUE)
>
> #Now Rasterize
> syn.ras<- rasterize(syn.ras.res@coords, raster(syn.ras.res), *?????*)
>
> I' ve like to create some kind of rule in *?????*, where pixel values in
> my final raster was 1 when inside de area (marks=area) and zero for
> outside de circular area given by points coordinates
> (syn.ras.res@coords) using 1 meter as resolution?
>
> Any ideas, please
>
> Thanks,
>
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

spatial autocorrelation in GAM residuals for large data set

Tue, 08/20/2019 - 08:46
Hello,

I have a large data set (~100k rows) containing observations at points (MODIS pixels) across the northern hemisphere.  I have created a GAM using the bam command in mgcv and I would like to check the model residuals for spatial autocorrelation.  

One idea is to use the DHARMa package (https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html#spatial-autocorrelation).  The code looks something like this:

    simulationOutput  <-   simulateResiduals(fittedModel = mymodel) # point at which R runs into memory problems
    testSpatialAutocorrelation(simulationOutput = simulationOutput, x =  data$latitude, y= data$longitude)

However, this runs into memory problems. 

Another idea is to use the following code, after this tutorial (http://www.flutterbys.com.au/stats/tut/tut8.4a.html):
    library(ape)
    library(fields)
    coords = cbind(data$longitude, data$latitude)     
    w = rdist(coords)  # point at which R runs into memory problems
    Moran.I(x = residuals(mymodel), w = w)

But this also runs into memory problems.  I have tried increasing the amount of memory allotted to R, but that just means R works for longer before timing out.  

So, two questions: (1) Is there a memory efficient way to check for spatial autocorrelation using Moran's I in large data sets? or (2) Is there another way to check for spatial autocorrelation (besides Moran's I) that won't have such memory problems?

Thanks in advance,

Elizabeth







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

Rasterize point process using specific rule

Tue, 08/20/2019 - 08:12
Dear members,

I've like to create a raster using my event coordinates and information
about the size in this coordinates as a rule for my raster creation.
First, I make:

library(raster)
library(spatstat)

#Points coordinates in UTM
xp<-c(371278.588,371250.722,371272.618,371328.421,371349.974,
371311.95,371296.265,371406.46,371411.551,371329.041,371338.081,
371334.182,371333.756,371299.818,371254.374,371193.673,371172.836,
371173.803,371153.73,371165.051,371140.417,371168.279,371166.367,
371180.575,371132.664,371129.791,371232.919,371208.502,371289.462,
371207.595,371219.008,371139.921,371133.215,371061.467,371053.69,
371099.897,371108.782,371112.52,371114.241,371176.236,371159.185,
371159.291,371158.552,370978.252,371120.03,371116.993)

yp<-c(8246507.94,8246493.176,8246465.974,8246464.413,8246403.465,
8246386.098,8246432.144,8246394.827,8246366.201,8246337.626,
8246311.125,8246300.039,8246299.594,8246298.072,8246379.351,
8246431.998,8246423.913,8246423.476,8246431.658,8246418.226,
8246400.161,8246396.891,8246394.225,8246400.391,8246370.244,
8246367.019,8246311.075,8246255.174,8246255.085,8246226.514,
8246215.847,8246337.316,8246330.197,8246311.197,8246304.183,
8246239.282,8246239.887,8246241.678,8246240.361,8246167.364,
8246171.581,8246171.803,8246169.807,8246293.57,8246183.194,8246189.926)

#Now I have the size of each nest in meters (marked process)

area<-c(117,30,4,341,15,160,35,280,108,168,63,143,2,48,182,42,
88,56,27,156,288,45,49,234,72,270,91,40,304,56,35,4,56.7,9,4.6,
105,133,135,23.92,190,12.9,15.2,192.78,104,255,24)

# Make a contour for the window creation
W <- convexhull.xy(xp,yp)

#Create ppm object
syn.ppp <- ppp(x=xp,y=yp,window=W, marks=area)

# Plot the point process
plot(syn.ppp)

#Definition of raster resolution - 1 meter square
ext <- as(extent(c(W$xrange,W$yrange)), "SpatialPolygons")
syn.ras.res<-rasterToPoints(raster(ext, resolution = 1), spatial = TRUE)

#Now Rasterize
syn.ras<- rasterize(syn.ras.res@coords, raster(syn.ras.res), *?????*)

I' ve like to create some kind of rule in *?????*, where pixel values in
my final raster was 1 when inside de area (marks=area) and zero for
outside de circular area given by points coordinates
(syn.ras.res@coords) using 1 meter as resolution?

Any ideas, please

Thanks,

--
======================================================================
Alexandre dos Santos
Prote??o Florestal
IFMT - Instituto Federal de Educa??o, Ci?ncia e Tecnologia de Mato Grosso
Campus C?ceres
Caixa Postal 244
Avenida dos Ramires, s/n
Bairro: Distrito Industrial
C?ceres - MT                      CEP: 78.200-000
Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO)

         [hidden email]
Lattes: http://lattes.cnpq.br/1360403201088680
OrcID: orcid.org/0000-0001-8232-6722
Researchgate: www.researchgate.net/profile/Alexandre_Santos10
LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635
Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/
======================================================================


        [[alternative HTML version deleted]]

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

Re: rgdal::writeOGR with driver='ESRI Shapefile' converts Polygon object into a hole

Tue, 08/20/2019 - 04:55
Good, thanks very much!

Roger

On Mon, 19 Aug 2019, Phil Haines wrote:

> Hi Roger,
>
> I have added an example to https://github.com/ateucher/rmapshaper/issues/89
>
> Hopefully it illustrates the behaviour I was seeing from
> rmapshaper::ms_simplify() but please let me know if not.
>
> And thanks for the nudge towards GPKG, I will investigate.
> Phil
>
> On Mon, 19 Aug 2019 at 15:48, Roger Bivand <[hidden email]> wrote:
>>
>> Hi Phil,
>>
>> On Sun, 18 Aug 2019, Phil Haines wrote:
>>
>>> Hi Roger,
>>>
>>> I originally encountered this issue while trying to reduce the size of
>>> German postal code boundaries. I used rmapshaper::ms_simplify() which
>>> introduced the polygons sharing a common edge. I definitely didn't
>>> need them to be separate polygons, and would happily have merged them
>>> had I been more familiar with gUnaryUnion().
>>
>> If the German postal code boundaries are open, could you please put (an
>> affected subset) on an URL, and provide a short script to replicate the
>> workflow. Then we can engage with the rmapshaper maintainer, and see
>> whether the underlying problem in the workflow is in the js library or its
>> R wrapper. I've opened: https://github.com/ateucher/rmapshaper/issues/89.
>>
>>>
>>> I apologise if my question has resulted in unnecessary effort on your
>>> part. I reported the behaviour because I found it surprising. However,
>>> I now understand that this example is not a valid simple feature
>>> geometry and that the solution is to take steps to ensure I have valid
>>> geometries, and to repair if not.
>>>
>>
>> No need to apologise!! Your reproducible example was excellent, without it
>> you wouldn't have contributed to getting the internal shapelib code
>> changed to make it less restrictive in future releases of GDAL,
>> potentially benefitting lots of users (who really should drop ESRI
>> shapefiles, and/or should check geometries for validity, but ... whose
>> work you've helped save). By the way, ESRI would be very happy if
>> shapefiles stopped being produced in current workflows, and that new files
>> should rather be GPKG.
>>
>>> Additionally I must confess that I only use the ESRI Shapefile driver
>>> because I don't know any better... My use cases are reading and
>>> writing from R (usually to produce maps) and occasionally opening in
>>> QGIS. I can happily switch to any format that supports this workflow.
>>>
>>
>> GPKG is now very broadly supported, is SF-compliant, and has the big
>> user-facing benefits of not having field name length restrictions, and
>> many fewer encoding issues for field names and string values.
>>
>> Best wishes,
>>
>> Roger
>>
>>> Thank you very much for your time,
>>> Phil
>>>
>>> On Sun, 18 Aug 2019 at 13:59, Roger Bivand <[hidden email]> wrote:
>>>>
>>>> The reason that the problem occurred is that a MULTIPOLYGON with two
>>>> exterior rings becomes invalid if the exterior rings touch along an edge
>>>> (this case). It is important to know the use case, to see whether:
>>>>
>>>> library(rgeos)
>>>> writeWKT(Ps1)
>>>> gIsValid(Ps1, reason=TRUE)
>>>> Ps1a <- gUnaryUnion(gBuffer(Ps1, width=0))
>>>> gIsValid(Ps1a, reason=TRUE)
>>>> writeWKT(Ps1a)
>>>>
>>>> or equivalently in sf for sf objects, should be applied before trying to
>>>> write the object out to file with this driver.
>>>>
>>>> However, because drivers that are compliant with the simple features
>>>> standard (which bans exterior rings sharing edges) have been permissive
>>>> and do round-trip this invalid object, a relaxation in the OGR ESRI
>>>> shapefile driver has been provided and may be included in a future
>>>> release.
>>>>
>>>> We need to know (see these issues):
>>>>
>>>> https://github.com/r-spatial/sf/issues/1130
>>>> https://github.com/OSGeo/gdal/issues/1787
>>>>
>>>> why it was desirable to write out this object using this driver? Could an
>>>> alternative driver have been used, or is ESRI shapefile the only format
>>>> used in the workflow?
>>>>
>>>> If it has to be this driver, could the workflow be changed to repair
>>>> degenerate cases before writing? If using sp classes, rgeos may be used to
>>>> test for and probably repair such geometries before they reach
>>>> rgdal::writeOGR() for this driver. Adding code to sf and rgdal to trap
>>>> degenerate cases does encumber all users with valid geometries with
>>>> the time wasted on extra checking, so building checks into sf and rgdal is
>>>> not desirable.
>>>>
>>>> I hope it is possible to find out more about the use case quickly, to pass
>>>> on to GDAL developers to help motivate a relaxation in their current
>>>> policy with regard to this driver, and to encourage them to include the
>>>> fix branch in a future release of GDAL.
>>>>
>>>> Roger
>>>>
>>>> On Sat, 17 Aug 2019, Roger Bivand wrote:
>>>>
>>>>> Please follow up both here and on:
>>>>>
>>>>> https://github.com/r-spatial/sf/issues/1130
>>>>>
>>>>> as the problem is also seen in the sf package using the same GDAL ESRI
>>>>> Shapefile driver.
>>>>>
>>>>> Roger
>>>>>
>>>>> On Fri, 16 Aug 2019, Roger Bivand wrote:
>>>>>
>>>>>>  On Fri, 16 Aug 2019, Roger Bivand wrote:
>>>>>>
>>>>>>>   On Tue, 13 Aug 2019, Phil Haines wrote:
>>>>>>>
>>>>>>>>    Dear list,
>>>>>>>>
>>>>>>>>    I have a single Polygons object containing multiple Polygon objects
>>>>>>>>    that share a common border. When I output this using writeOGR() one of
>>>>>>>>    the Polygon objects becomes a hole, as the following example shows.
>>>>>>>>
>>>>>>>>    Create a Polygons object containing two adjoining Polygon objects
>>>>>>>>>    library(rgdal)
>>>>>>>>>    r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
>>>>>>>>>    r2 <- r1; r2[,1] <- r2[,1]+1
>>>>>>>>>    Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
>>>>>>>>>    SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
>>>>>>>>>    data.frame(Example=c("Minimal")))
>>>>>>>>
>>>>>>>>    Perform a write/readOGR() cycle
>>>>>>>>>    fn <- tempfile()
>>>>>>>>>    writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
>>>>>>>>>    SPDF2 <- readOGR(dsn=fn,layer='test')
>>>>>>>>
>>>>>>>>    Second Polygon object is now a hole
>>>>>>>>>    sapply(SPDF2@polygons[[1]]@Polygons,slot,"hole")
>>>>>>>>    [1] FALSE  TRUE
>>>>>>>>
>>>>>>>>    I see from the sp documentation that "Polygon objects belonging to a
>>>>>>>>    Polygons object should either not overlap one-other, or should be
>>>>>>>>    fully included" but I am not sure how this relates to bordering
>>>>>>>>    Polygon objects. I would welcome any advice as to whether what I am
>>>>>>>>    asking of writeOGR is reasonable?
>>>>>>>
>>>>>>>   The problem is with the 'ESRI Shapefile' representation and driver:
>>>>>>>
>>>>>>>   library(rgdal)
>>>>>>>   r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
>>>>>>>   r2 <- r1; r2[,1] <- r2[,1]+1
>>>>>>>   Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
>>>>>>>   SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
>>>>>>>   data.frame(Example=c("Minimal")))
>>>>>>>   sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "hole")
>>>>>>>   sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>>>>>>
>>>>>>>   # which constructs a MULTIPOLYGON object:
>>>>>>>
>>>>>>>   rgeos::writeWKT(SPDF)
>>>>>>>   library(sf)
>>>>>>>   st_as_text(st_geometry(st_as_sf(SPDF)))
>>>>>>>
>>>>>>> #    The 'ESRI Shapefile' driver is not Simple-Feature compliant (it
>>>>>>> #    pre-dates it), so the failure occurs by seeing the second exterior
>>>>>>> #    ring
>>>>>>> #    as an interior ring
>>>>>>>
>>>>>>>   fn <- tempfile()
>>>>>>>   writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
>>>>>>>   SPDF2 <- readOGR(dsn=fn, layer='test')
>>>>>>>   sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "hole")
>>>>>>>   sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>>>>>>   rgeos::writeWKT(SPDF2)
>>>>>>>   st_as_text(st_geometry(st_as_sf(SPDF2)))
>>>>>>>
>>>>>>>   # This happens with sf too, using the same GDAL driver:
>>>>>>>
>>>>>>>   sf2 <- st_read(dsn=fn, layer='test')
>>>>>>>   st_geometry(sf2)
>>>>>>>   library(sf)
>>>>>>>   st_as_text(st_geometry(st_as_sf(sf2)))
>>>>>>>   rgeos::writeWKT(as(sf2, "Spatial"))
>>>>>>>
>>>>>>>   # Adding the comment fix doesn't help:
>>>>>>>
>>>>>>>   comment(slot(SPDF, "polygons")[[1]])
>>>>>>>   SPDF_c <- rgeos::createSPComment(SPDF)
>>>>>>>   comment(slot(SPDF_c, "polygons")[[1]])
>>>>>>>   writeOGR(SPDF_c, fn, layer='test_c', driver='ESRI Shapefile',
>>>>>>>     verbose=TRUE)
>>>>>>>
>>>>>>> #    reports
>>>>>>> #    Object initially classed as: wkbPolygon
>>>>>>> #    SFS comments in Polygons objects
>>>>>>> #    Object reclassed as: wkbMultiPolygon
>>>>>>>
>>>>>>>   SPDF2_c <- readOGR(dsn=fn, layer='test_c')
>>>>>>>   sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot, "hole")
>>>>>>>   sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot,
>>>>>>>   "ringDir")
>>>>>>>   rgeos::writeWKT(SPDF2_c)
>>>>>>>   st_as_text(st_geometry(st_as_sf(SPDF2_c)))
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>   # If the input object is written out with the GeoPackage driver:
>>>>>>>
>>>>>>>   fn1 <- tempfile(fileext=".gpkg")
>>>>>>>   writeOGR(SPDF, fn1, layer="test", driver='GPKG')
>>>>>>>   sf2a <- st_read(dsn=fn1,layer='test')
>>>>>>>   st_coordinates(st_geometry(sf2a))
>>>>>>>   SPDF2a <- readOGR(dsn=fn1)
>>>>>>>   sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "hole")
>>>>>>>   sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>>>>>>   rgeos::writeWKT(SPDF2a)
>>>>>>>   st_as_text(st_geometry(st_as_sf(SPDF2a)))
>>>>>>>
>>>>>>>   # the issue is resolved. If we separate the exterior rings:
>>>>>>>
>>>>>>>   r2a <- r1; r2a[,1] <- r2a[,1]+1.00001
>>>>>>>   Ps1a = Polygons(list(Polygon(r1),Polygon(r2a)),ID=1)
>>>>>>>   SPDFa = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1a)),
>>>>>>>   data.frame(Example=c("Minimal")))
>>>>>>>   fna <- tempfile()
>>>>>>>   writeOGR(SPDFa, fna, layer='test', driver='ESRI Shapefile')
>>>>>>>   SPDF2_a <- readOGR(dsn=fna, layer='test')
>>>>>>>   sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot, "hole")
>>>>>>>   sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot,
>>>>>>>   "ringDir")
>>>>>>>   rgeos::writeWKT(SPDF2_a)
>>>>>>>   st_as_text(st_geometry(st_as_sf(SPDF2_a)))
>>>>>>>
>>>>>>>   # we are OK as the two exterior rings do not touch.
>>>>>>>
>>>>>>>   # does using sf make a difference?
>>>>>>>
>>>>>>>   fn_s <- tempfile(fileext=".shp")
>>>>>>>   st_write(st_as_sf(SPDF), dsn=fn_s)
>>>>>>>   sf_in <- st_read(fn_s)
>>>>>>>   st_as_text(st_geometry(st_as_sf(sf_in)))
>>>>>>>
>>>>>>>   # No
>>>>>>>
>>>>>>>   fn_s <- tempfile(fileext=".shp")
>>>>>>>   st_write(st_as_sf(SPDF_c), dsn=fn_s)
>>>>>>>   sf_in_c <- st_read(fn_s)
>>>>>>>   st_as_text(st_geometry(st_as_sf(sf_in_c)))
>>>>>>>
>>>>>>>   # nor with the pretend-SF-compliant comment set either.
>>>>>>>
>>>>>>>   So the weakness is in the "ESRI Shapefile" write driver, or possibly in
>>>>>>>   the OGRGeometryFactory::organizePolygons() function in GDAL used in
>>>>>>>   OGR_write() (a C++ function) called by writeOGR(). If sf::st_write()
>>>>>>>   also
>>>>>>>   calls OGRGeometryFactory::organizePolygons(), we'd maybe consider that
>>>>>>>   it
>>>>>>>   has a weakness for the "ESRI Shapefile" driver, but which does not
>>>>>>>   affect
>>>>>>>   SF-compliant drivers.
>>>>>>
>>>>>>  Without the comment set, OGRGeometryFactory::organizePolygons() is used;
>>>>>>  with it set, OGRGeometryFactory::organizePolygons() is not used, because
>>>>>>  the object is declared as two exterior rings. In both cases, we have the
>>>>>>  output object written out and read back in incorrectly with the ESRI
>>>>>>  shapefile driver, but SF-compliant drivers round-trip (in the test
>>>>>>  GeoJSON), correctly.
>>>>>>
>>>>>>  It is likely that the changes made in 2015 to accommodate GeoJSON led to
>>>>>>  this possible regression for the ESRI Shapefile driver. I'm adding this
>>>>>>  geometry to tests/tripup.R and data files; without code changes the hole
>>>>>>  slot is wrong and the ring direction is changes to match, so the
>>>>>>  coordinates change too.
>>>>>>
>>>>>>  Reading using the deprecated maptools::readShapeSpatial() also gets a hole
>>>>>>  in the second external ring. However, writing with deprecated
>>>>>> maptools::  writeSpatialShape() yields a shapefile that when read with
>>>>>> maptools::  readShapeSpatial() gives the correct two exterior ring, no hole
>>>>>>  object. When read with sf::st_read() and rgdal::readOGR(), the object is
>>>>>>  also correct. So the problem definitely lies in rgdal::writeOGR(), and
>>>>>>  sf::st_write() - roundtripping with sf::st_write() and sf::st_read()
>>>>>>  degrades from MULTIPOLYGON to POLYGON with the ESRI Shapefile driver.
>>>>>>
>>>>>>  Roger
>>>>>>
>>>>>>>
>>>>>>>   This is probably related to a similar but inverse problem with the
>>>>>>>   SF-compliant GeoJSON driver in 2015:
>>>>>>>
>>>>>>>   https://stat.ethz.ch/pipermail/r-sig-geo/2015-October/023609.html
>>>>>>>
>>>>>>>   continued the next month in:
>>>>>>>
>>>>>>>   https://stat.ethz.ch/pipermail/r-sig-geo/2015-November/023656.html
>>>>>>>
>>>>>>>   The details are in this SVN diff
>>>>>>>
>>>>>>>   https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=555&r2=571
>>>>>>>
>>>>>>>   up to the end og the list thread, and this one from then until now:
>>>>>>>
>>>>>>>   https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=571&r2=733
>>>>>>>
>>>>>>>   Summary: could you change drivers, or is it really necessary to fix an
>>>>>>>   EOL
>>>>>>>   problem? What is your use case?
>>>>>>>
>>>>>>>>
>>>>>>>>    Thanks in advance for your time,
>>>>>>>
>>>>>>>   Thanks for a complete example,
>>>>>>>
>>>>>>>   Roger
>>>>>>>
>>>>>>>>    Phil
>>>>>>>>
>>>>>>>>>    sessionInfo()
>>>>>>>>    R version 3.5.1 (2018-07-02)
>>>>>>>>    Platform: x86_64-w64-mingw32/x64 (64-bit)
>>>>>>>>    Running under: Windows >= 8 x64 (build 9200)
>>>>>>>>
>>>>>>>>    Matrix products: default
>>>>>>>>
>>>>>>>>    locale:
>>>>>>>>    [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United
>>>>>>>>    Kingdom.1252    LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
>>>>>>>>                            LC_TIME=English_United Kingdom.1252
>>>>>>>>
>>>>>>>>    attached base packages:
>>>>>>>>    [1] stats     graphics  grDevices utils     datasets  methods   base
>>>>>>>>
>>>>>>>>    other attached packages:
>>>>>>>>    [1] rgdal_1.4-4 sp_1.3-1
>>>>>>>>
>>>>>>>>    loaded via a namespace (and not attached):
>>>>>>>>    [1] compiler_3.5.1  tools_3.5.1     yaml_2.2.0      grid_3.5.1
>>>>>>>>    lattice_0.20-35
>>>>>>>>
>>>>>>>>    _______________________________________________
>>>>>>>>    R-sig-Geo mailing list
>>>>>>>>    [hidden email]
>>>>>>>>    https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>
>>>> --
>>>> Roger Bivand
>>>> Department of Economics, Norwegian School of Economics,
>>>> Helleveien 30, N-5045 Bergen, Norway.
>>>> voice: +47 55 95 93 55; e-mail: [hidden email]
>>>> https://orcid.org/0000-0003-2392-6140
>>>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>>
>>
>> --
>> Roger Bivand
>> Department of Economics, Norwegian School of Economics,
>> Helleveien 30, N-5045 Bergen, Norway.
>> voice: +47 55 95 93 55; e-mail: [hidden email]
>> https://orcid.org/0000-0003-2392-6140
>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Re: rgdal::writeOGR with driver='ESRI Shapefile' converts Polygon object into a hole

Mon, 08/19/2019 - 11:50
Hi Roger,

I have added an example to https://github.com/ateucher/rmapshaper/issues/89

Hopefully it illustrates the behaviour I was seeing from
rmapshaper::ms_simplify() but please let me know if not.

And thanks for the nudge towards GPKG, I will investigate.
Phil

On Mon, 19 Aug 2019 at 15:48, Roger Bivand <[hidden email]> wrote:
>
> Hi Phil,
>
> On Sun, 18 Aug 2019, Phil Haines wrote:
>
> > Hi Roger,
> >
> > I originally encountered this issue while trying to reduce the size of
> > German postal code boundaries. I used rmapshaper::ms_simplify() which
> > introduced the polygons sharing a common edge. I definitely didn't
> > need them to be separate polygons, and would happily have merged them
> > had I been more familiar with gUnaryUnion().
>
> If the German postal code boundaries are open, could you please put (an
> affected subset) on an URL, and provide a short script to replicate the
> workflow. Then we can engage with the rmapshaper maintainer, and see
> whether the underlying problem in the workflow is in the js library or its
> R wrapper. I've opened: https://github.com/ateucher/rmapshaper/issues/89.
>
> >
> > I apologise if my question has resulted in unnecessary effort on your
> > part. I reported the behaviour because I found it surprising. However,
> > I now understand that this example is not a valid simple feature
> > geometry and that the solution is to take steps to ensure I have valid
> > geometries, and to repair if not.
> >
>
> No need to apologise!! Your reproducible example was excellent, without it
> you wouldn't have contributed to getting the internal shapelib code
> changed to make it less restrictive in future releases of GDAL,
> potentially benefitting lots of users (who really should drop ESRI
> shapefiles, and/or should check geometries for validity, but ... whose
> work you've helped save). By the way, ESRI would be very happy if
> shapefiles stopped being produced in current workflows, and that new files
> should rather be GPKG.
>
> > Additionally I must confess that I only use the ESRI Shapefile driver
> > because I don't know any better... My use cases are reading and
> > writing from R (usually to produce maps) and occasionally opening in
> > QGIS. I can happily switch to any format that supports this workflow.
> >
>
> GPKG is now very broadly supported, is SF-compliant, and has the big
> user-facing benefits of not having field name length restrictions, and
> many fewer encoding issues for field names and string values.
>
> Best wishes,
>
> Roger
>
> > Thank you very much for your time,
> > Phil
> >
> > On Sun, 18 Aug 2019 at 13:59, Roger Bivand <[hidden email]> wrote:
> >>
> >> The reason that the problem occurred is that a MULTIPOLYGON with two
> >> exterior rings becomes invalid if the exterior rings touch along an edge
> >> (this case). It is important to know the use case, to see whether:
> >>
> >> library(rgeos)
> >> writeWKT(Ps1)
> >> gIsValid(Ps1, reason=TRUE)
> >> Ps1a <- gUnaryUnion(gBuffer(Ps1, width=0))
> >> gIsValid(Ps1a, reason=TRUE)
> >> writeWKT(Ps1a)
> >>
> >> or equivalently in sf for sf objects, should be applied before trying to
> >> write the object out to file with this driver.
> >>
> >> However, because drivers that are compliant with the simple features
> >> standard (which bans exterior rings sharing edges) have been permissive
> >> and do round-trip this invalid object, a relaxation in the OGR ESRI
> >> shapefile driver has been provided and may be included in a future
> >> release.
> >>
> >> We need to know (see these issues):
> >>
> >> https://github.com/r-spatial/sf/issues/1130
> >> https://github.com/OSGeo/gdal/issues/1787
> >>
> >> why it was desirable to write out this object using this driver? Could an
> >> alternative driver have been used, or is ESRI shapefile the only format
> >> used in the workflow?
> >>
> >> If it has to be this driver, could the workflow be changed to repair
> >> degenerate cases before writing? If using sp classes, rgeos may be used to
> >> test for and probably repair such geometries before they reach
> >> rgdal::writeOGR() for this driver. Adding code to sf and rgdal to trap
> >> degenerate cases does encumber all users with valid geometries with
> >> the time wasted on extra checking, so building checks into sf and rgdal is
> >> not desirable.
> >>
> >> I hope it is possible to find out more about the use case quickly, to pass
> >> on to GDAL developers to help motivate a relaxation in their current
> >> policy with regard to this driver, and to encourage them to include the
> >> fix branch in a future release of GDAL.
> >>
> >> Roger
> >>
> >> On Sat, 17 Aug 2019, Roger Bivand wrote:
> >>
> >>> Please follow up both here and on:
> >>>
> >>> https://github.com/r-spatial/sf/issues/1130
> >>>
> >>> as the problem is also seen in the sf package using the same GDAL ESRI
> >>> Shapefile driver.
> >>>
> >>> Roger
> >>>
> >>> On Fri, 16 Aug 2019, Roger Bivand wrote:
> >>>
> >>>>  On Fri, 16 Aug 2019, Roger Bivand wrote:
> >>>>
> >>>>>   On Tue, 13 Aug 2019, Phil Haines wrote:
> >>>>>
> >>>>>>    Dear list,
> >>>>>>
> >>>>>>    I have a single Polygons object containing multiple Polygon objects
> >>>>>>    that share a common border. When I output this using writeOGR() one of
> >>>>>>    the Polygon objects becomes a hole, as the following example shows.
> >>>>>>
> >>>>>>    Create a Polygons object containing two adjoining Polygon objects
> >>>>>>>    library(rgdal)
> >>>>>>>    r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
> >>>>>>>    r2 <- r1; r2[,1] <- r2[,1]+1
> >>>>>>>    Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
> >>>>>>>    SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
> >>>>>>>    data.frame(Example=c("Minimal")))
> >>>>>>
> >>>>>>    Perform a write/readOGR() cycle
> >>>>>>>    fn <- tempfile()
> >>>>>>>    writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
> >>>>>>>    SPDF2 <- readOGR(dsn=fn,layer='test')
> >>>>>>
> >>>>>>    Second Polygon object is now a hole
> >>>>>>>    sapply(SPDF2@polygons[[1]]@Polygons,slot,"hole")
> >>>>>>    [1] FALSE  TRUE
> >>>>>>
> >>>>>>    I see from the sp documentation that "Polygon objects belonging to a
> >>>>>>    Polygons object should either not overlap one-other, or should be
> >>>>>>    fully included" but I am not sure how this relates to bordering
> >>>>>>    Polygon objects. I would welcome any advice as to whether what I am
> >>>>>>    asking of writeOGR is reasonable?
> >>>>>
> >>>>>   The problem is with the 'ESRI Shapefile' representation and driver:
> >>>>>
> >>>>>   library(rgdal)
> >>>>>   r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
> >>>>>   r2 <- r1; r2[,1] <- r2[,1]+1
> >>>>>   Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
> >>>>>   SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
> >>>>>   data.frame(Example=c("Minimal")))
> >>>>>   sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "hole")
> >>>>>   sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "ringDir")
> >>>>>
> >>>>>   # which constructs a MULTIPOLYGON object:
> >>>>>
> >>>>>   rgeos::writeWKT(SPDF)
> >>>>>   library(sf)
> >>>>>   st_as_text(st_geometry(st_as_sf(SPDF)))
> >>>>>
> >>>>> #    The 'ESRI Shapefile' driver is not Simple-Feature compliant (it
> >>>>> #    pre-dates it), so the failure occurs by seeing the second exterior
> >>>>> #    ring
> >>>>> #    as an interior ring
> >>>>>
> >>>>>   fn <- tempfile()
> >>>>>   writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
> >>>>>   SPDF2 <- readOGR(dsn=fn, layer='test')
> >>>>>   sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "hole")
> >>>>>   sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "ringDir")
> >>>>>   rgeos::writeWKT(SPDF2)
> >>>>>   st_as_text(st_geometry(st_as_sf(SPDF2)))
> >>>>>
> >>>>>   # This happens with sf too, using the same GDAL driver:
> >>>>>
> >>>>>   sf2 <- st_read(dsn=fn, layer='test')
> >>>>>   st_geometry(sf2)
> >>>>>   library(sf)
> >>>>>   st_as_text(st_geometry(st_as_sf(sf2)))
> >>>>>   rgeos::writeWKT(as(sf2, "Spatial"))
> >>>>>
> >>>>>   # Adding the comment fix doesn't help:
> >>>>>
> >>>>>   comment(slot(SPDF, "polygons")[[1]])
> >>>>>   SPDF_c <- rgeos::createSPComment(SPDF)
> >>>>>   comment(slot(SPDF_c, "polygons")[[1]])
> >>>>>   writeOGR(SPDF_c, fn, layer='test_c', driver='ESRI Shapefile',
> >>>>>     verbose=TRUE)
> >>>>>
> >>>>> #    reports
> >>>>> #    Object initially classed as: wkbPolygon
> >>>>> #    SFS comments in Polygons objects
> >>>>> #    Object reclassed as: wkbMultiPolygon
> >>>>>
> >>>>>   SPDF2_c <- readOGR(dsn=fn, layer='test_c')
> >>>>>   sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot, "hole")
> >>>>>   sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot,
> >>>>>   "ringDir")
> >>>>>   rgeos::writeWKT(SPDF2_c)
> >>>>>   st_as_text(st_geometry(st_as_sf(SPDF2_c)))
> >>>>>
> >>>>>
> >>>>>
> >>>>>   # If the input object is written out with the GeoPackage driver:
> >>>>>
> >>>>>   fn1 <- tempfile(fileext=".gpkg")
> >>>>>   writeOGR(SPDF, fn1, layer="test", driver='GPKG')
> >>>>>   sf2a <- st_read(dsn=fn1,layer='test')
> >>>>>   st_coordinates(st_geometry(sf2a))
> >>>>>   SPDF2a <- readOGR(dsn=fn1)
> >>>>>   sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "hole")
> >>>>>   sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "ringDir")
> >>>>>   rgeos::writeWKT(SPDF2a)
> >>>>>   st_as_text(st_geometry(st_as_sf(SPDF2a)))
> >>>>>
> >>>>>   # the issue is resolved. If we separate the exterior rings:
> >>>>>
> >>>>>   r2a <- r1; r2a[,1] <- r2a[,1]+1.00001
> >>>>>   Ps1a = Polygons(list(Polygon(r1),Polygon(r2a)),ID=1)
> >>>>>   SPDFa = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1a)),
> >>>>>   data.frame(Example=c("Minimal")))
> >>>>>   fna <- tempfile()
> >>>>>   writeOGR(SPDFa, fna, layer='test', driver='ESRI Shapefile')
> >>>>>   SPDF2_a <- readOGR(dsn=fna, layer='test')
> >>>>>   sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot, "hole")
> >>>>>   sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot,
> >>>>>   "ringDir")
> >>>>>   rgeos::writeWKT(SPDF2_a)
> >>>>>   st_as_text(st_geometry(st_as_sf(SPDF2_a)))
> >>>>>
> >>>>>   # we are OK as the two exterior rings do not touch.
> >>>>>
> >>>>>   # does using sf make a difference?
> >>>>>
> >>>>>   fn_s <- tempfile(fileext=".shp")
> >>>>>   st_write(st_as_sf(SPDF), dsn=fn_s)
> >>>>>   sf_in <- st_read(fn_s)
> >>>>>   st_as_text(st_geometry(st_as_sf(sf_in)))
> >>>>>
> >>>>>   # No
> >>>>>
> >>>>>   fn_s <- tempfile(fileext=".shp")
> >>>>>   st_write(st_as_sf(SPDF_c), dsn=fn_s)
> >>>>>   sf_in_c <- st_read(fn_s)
> >>>>>   st_as_text(st_geometry(st_as_sf(sf_in_c)))
> >>>>>
> >>>>>   # nor with the pretend-SF-compliant comment set either.
> >>>>>
> >>>>>   So the weakness is in the "ESRI Shapefile" write driver, or possibly in
> >>>>>   the OGRGeometryFactory::organizePolygons() function in GDAL used in
> >>>>>   OGR_write() (a C++ function) called by writeOGR(). If sf::st_write()
> >>>>>   also
> >>>>>   calls OGRGeometryFactory::organizePolygons(), we'd maybe consider that
> >>>>>   it
> >>>>>   has a weakness for the "ESRI Shapefile" driver, but which does not
> >>>>>   affect
> >>>>>   SF-compliant drivers.
> >>>>
> >>>>  Without the comment set, OGRGeometryFactory::organizePolygons() is used;
> >>>>  with it set, OGRGeometryFactory::organizePolygons() is not used, because
> >>>>  the object is declared as two exterior rings. In both cases, we have the
> >>>>  output object written out and read back in incorrectly with the ESRI
> >>>>  shapefile driver, but SF-compliant drivers round-trip (in the test
> >>>>  GeoJSON), correctly.
> >>>>
> >>>>  It is likely that the changes made in 2015 to accommodate GeoJSON led to
> >>>>  this possible regression for the ESRI Shapefile driver. I'm adding this
> >>>>  geometry to tests/tripup.R and data files; without code changes the hole
> >>>>  slot is wrong and the ring direction is changes to match, so the
> >>>>  coordinates change too.
> >>>>
> >>>>  Reading using the deprecated maptools::readShapeSpatial() also gets a hole
> >>>>  in the second external ring. However, writing with deprecated
> >>>> maptools::  writeSpatialShape() yields a shapefile that when read with
> >>>> maptools::  readShapeSpatial() gives the correct two exterior ring, no hole
> >>>>  object. When read with sf::st_read() and rgdal::readOGR(), the object is
> >>>>  also correct. So the problem definitely lies in rgdal::writeOGR(), and
> >>>>  sf::st_write() - roundtripping with sf::st_write() and sf::st_read()
> >>>>  degrades from MULTIPOLYGON to POLYGON with the ESRI Shapefile driver.
> >>>>
> >>>>  Roger
> >>>>
> >>>>>
> >>>>>   This is probably related to a similar but inverse problem with the
> >>>>>   SF-compliant GeoJSON driver in 2015:
> >>>>>
> >>>>>   https://stat.ethz.ch/pipermail/r-sig-geo/2015-October/023609.html
> >>>>>
> >>>>>   continued the next month in:
> >>>>>
> >>>>>   https://stat.ethz.ch/pipermail/r-sig-geo/2015-November/023656.html
> >>>>>
> >>>>>   The details are in this SVN diff
> >>>>>
> >>>>>   https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=555&r2=571
> >>>>>
> >>>>>   up to the end og the list thread, and this one from then until now:
> >>>>>
> >>>>>   https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=571&r2=733
> >>>>>
> >>>>>   Summary: could you change drivers, or is it really necessary to fix an
> >>>>>   EOL
> >>>>>   problem? What is your use case?
> >>>>>
> >>>>>>
> >>>>>>    Thanks in advance for your time,
> >>>>>
> >>>>>   Thanks for a complete example,
> >>>>>
> >>>>>   Roger
> >>>>>
> >>>>>>    Phil
> >>>>>>
> >>>>>>>    sessionInfo()
> >>>>>>    R version 3.5.1 (2018-07-02)
> >>>>>>    Platform: x86_64-w64-mingw32/x64 (64-bit)
> >>>>>>    Running under: Windows >= 8 x64 (build 9200)
> >>>>>>
> >>>>>>    Matrix products: default
> >>>>>>
> >>>>>>    locale:
> >>>>>>    [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United
> >>>>>>    Kingdom.1252    LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
> >>>>>>                            LC_TIME=English_United Kingdom.1252
> >>>>>>
> >>>>>>    attached base packages:
> >>>>>>    [1] stats     graphics  grDevices utils     datasets  methods   base
> >>>>>>
> >>>>>>    other attached packages:
> >>>>>>    [1] rgdal_1.4-4 sp_1.3-1
> >>>>>>
> >>>>>>    loaded via a namespace (and not attached):
> >>>>>>    [1] compiler_3.5.1  tools_3.5.1     yaml_2.2.0      grid_3.5.1
> >>>>>>    lattice_0.20-35
> >>>>>>
> >>>>>>    _______________________________________________
> >>>>>>    R-sig-Geo mailing list
> >>>>>>    [hidden email]
> >>>>>>    https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>>>>>
> >>>>>
> >>>>>
> >>>>
> >>>>
> >>>
> >>>
> >>
> >> --
> >> Roger Bivand
> >> Department of Economics, Norwegian School of Economics,
> >> Helleveien 30, N-5045 Bergen, Norway.
> >> voice: +47 55 95 93 55; e-mail: [hidden email]
> >> https://orcid.org/0000-0003-2392-6140
> >> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
> >
>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: [hidden email]
> https://orcid.org/0000-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: rgdal::writeOGR with driver='ESRI Shapefile' converts Polygon object into a hole

Mon, 08/19/2019 - 09:48
Hi Phil,

On Sun, 18 Aug 2019, Phil Haines wrote:

> Hi Roger,
>
> I originally encountered this issue while trying to reduce the size of
> German postal code boundaries. I used rmapshaper::ms_simplify() which
> introduced the polygons sharing a common edge. I definitely didn't
> need them to be separate polygons, and would happily have merged them
> had I been more familiar with gUnaryUnion().

If the German postal code boundaries are open, could you please put (an
affected subset) on an URL, and provide a short script to replicate the
workflow. Then we can engage with the rmapshaper maintainer, and see
whether the underlying problem in the workflow is in the js library or its
R wrapper. I've opened: https://github.com/ateucher/rmapshaper/issues/89.

>
> I apologise if my question has resulted in unnecessary effort on your
> part. I reported the behaviour because I found it surprising. However,
> I now understand that this example is not a valid simple feature
> geometry and that the solution is to take steps to ensure I have valid
> geometries, and to repair if not.
>

No need to apologise!! Your reproducible example was excellent, without it
you wouldn't have contributed to getting the internal shapelib code
changed to make it less restrictive in future releases of GDAL,
potentially benefitting lots of users (who really should drop ESRI
shapefiles, and/or should check geometries for validity, but ... whose
work you've helped save). By the way, ESRI would be very happy if
shapefiles stopped being produced in current workflows, and that new files
should rather be GPKG.

> Additionally I must confess that I only use the ESRI Shapefile driver
> because I don't know any better... My use cases are reading and
> writing from R (usually to produce maps) and occasionally opening in
> QGIS. I can happily switch to any format that supports this workflow.
>

GPKG is now very broadly supported, is SF-compliant, and has the big
user-facing benefits of not having field name length restrictions, and
many fewer encoding issues for field names and string values.

Best wishes,

Roger

> Thank you very much for your time,
> Phil
>
> On Sun, 18 Aug 2019 at 13:59, Roger Bivand <[hidden email]> wrote:
>>
>> The reason that the problem occurred is that a MULTIPOLYGON with two
>> exterior rings becomes invalid if the exterior rings touch along an edge
>> (this case). It is important to know the use case, to see whether:
>>
>> library(rgeos)
>> writeWKT(Ps1)
>> gIsValid(Ps1, reason=TRUE)
>> Ps1a <- gUnaryUnion(gBuffer(Ps1, width=0))
>> gIsValid(Ps1a, reason=TRUE)
>> writeWKT(Ps1a)
>>
>> or equivalently in sf for sf objects, should be applied before trying to
>> write the object out to file with this driver.
>>
>> However, because drivers that are compliant with the simple features
>> standard (which bans exterior rings sharing edges) have been permissive
>> and do round-trip this invalid object, a relaxation in the OGR ESRI
>> shapefile driver has been provided and may be included in a future
>> release.
>>
>> We need to know (see these issues):
>>
>> https://github.com/r-spatial/sf/issues/1130
>> https://github.com/OSGeo/gdal/issues/1787
>>
>> why it was desirable to write out this object using this driver? Could an
>> alternative driver have been used, or is ESRI shapefile the only format
>> used in the workflow?
>>
>> If it has to be this driver, could the workflow be changed to repair
>> degenerate cases before writing? If using sp classes, rgeos may be used to
>> test for and probably repair such geometries before they reach
>> rgdal::writeOGR() for this driver. Adding code to sf and rgdal to trap
>> degenerate cases does encumber all users with valid geometries with
>> the time wasted on extra checking, so building checks into sf and rgdal is
>> not desirable.
>>
>> I hope it is possible to find out more about the use case quickly, to pass
>> on to GDAL developers to help motivate a relaxation in their current
>> policy with regard to this driver, and to encourage them to include the
>> fix branch in a future release of GDAL.
>>
>> Roger
>>
>> On Sat, 17 Aug 2019, Roger Bivand wrote:
>>
>>> Please follow up both here and on:
>>>
>>> https://github.com/r-spatial/sf/issues/1130
>>>
>>> as the problem is also seen in the sf package using the same GDAL ESRI
>>> Shapefile driver.
>>>
>>> Roger
>>>
>>> On Fri, 16 Aug 2019, Roger Bivand wrote:
>>>
>>>>  On Fri, 16 Aug 2019, Roger Bivand wrote:
>>>>
>>>>>   On Tue, 13 Aug 2019, Phil Haines wrote:
>>>>>
>>>>>>    Dear list,
>>>>>>
>>>>>>    I have a single Polygons object containing multiple Polygon objects
>>>>>>    that share a common border. When I output this using writeOGR() one of
>>>>>>    the Polygon objects becomes a hole, as the following example shows.
>>>>>>
>>>>>>    Create a Polygons object containing two adjoining Polygon objects
>>>>>>>    library(rgdal)
>>>>>>>    r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
>>>>>>>    r2 <- r1; r2[,1] <- r2[,1]+1
>>>>>>>    Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
>>>>>>>    SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
>>>>>>>    data.frame(Example=c("Minimal")))
>>>>>>
>>>>>>    Perform a write/readOGR() cycle
>>>>>>>    fn <- tempfile()
>>>>>>>    writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
>>>>>>>    SPDF2 <- readOGR(dsn=fn,layer='test')
>>>>>>
>>>>>>    Second Polygon object is now a hole
>>>>>>>    sapply(SPDF2@polygons[[1]]@Polygons,slot,"hole")
>>>>>>    [1] FALSE  TRUE
>>>>>>
>>>>>>    I see from the sp documentation that "Polygon objects belonging to a
>>>>>>    Polygons object should either not overlap one-other, or should be
>>>>>>    fully included" but I am not sure how this relates to bordering
>>>>>>    Polygon objects. I would welcome any advice as to whether what I am
>>>>>>    asking of writeOGR is reasonable?
>>>>>
>>>>>   The problem is with the 'ESRI Shapefile' representation and driver:
>>>>>
>>>>>   library(rgdal)
>>>>>   r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
>>>>>   r2 <- r1; r2[,1] <- r2[,1]+1
>>>>>   Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
>>>>>   SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
>>>>>   data.frame(Example=c("Minimal")))
>>>>>   sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "hole")
>>>>>   sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>>>>
>>>>>   # which constructs a MULTIPOLYGON object:
>>>>>
>>>>>   rgeos::writeWKT(SPDF)
>>>>>   library(sf)
>>>>>   st_as_text(st_geometry(st_as_sf(SPDF)))
>>>>>
>>>>> #    The 'ESRI Shapefile' driver is not Simple-Feature compliant (it
>>>>> #    pre-dates it), so the failure occurs by seeing the second exterior
>>>>> #    ring
>>>>> #    as an interior ring
>>>>>
>>>>>   fn <- tempfile()
>>>>>   writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
>>>>>   SPDF2 <- readOGR(dsn=fn, layer='test')
>>>>>   sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "hole")
>>>>>   sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>>>>   rgeos::writeWKT(SPDF2)
>>>>>   st_as_text(st_geometry(st_as_sf(SPDF2)))
>>>>>
>>>>>   # This happens with sf too, using the same GDAL driver:
>>>>>
>>>>>   sf2 <- st_read(dsn=fn, layer='test')
>>>>>   st_geometry(sf2)
>>>>>   library(sf)
>>>>>   st_as_text(st_geometry(st_as_sf(sf2)))
>>>>>   rgeos::writeWKT(as(sf2, "Spatial"))
>>>>>
>>>>>   # Adding the comment fix doesn't help:
>>>>>
>>>>>   comment(slot(SPDF, "polygons")[[1]])
>>>>>   SPDF_c <- rgeos::createSPComment(SPDF)
>>>>>   comment(slot(SPDF_c, "polygons")[[1]])
>>>>>   writeOGR(SPDF_c, fn, layer='test_c', driver='ESRI Shapefile',
>>>>>     verbose=TRUE)
>>>>>
>>>>> #    reports
>>>>> #    Object initially classed as: wkbPolygon
>>>>> #    SFS comments in Polygons objects
>>>>> #    Object reclassed as: wkbMultiPolygon
>>>>>
>>>>>   SPDF2_c <- readOGR(dsn=fn, layer='test_c')
>>>>>   sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot, "hole")
>>>>>   sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot,
>>>>>   "ringDir")
>>>>>   rgeos::writeWKT(SPDF2_c)
>>>>>   st_as_text(st_geometry(st_as_sf(SPDF2_c)))
>>>>>
>>>>>
>>>>>
>>>>>   # If the input object is written out with the GeoPackage driver:
>>>>>
>>>>>   fn1 <- tempfile(fileext=".gpkg")
>>>>>   writeOGR(SPDF, fn1, layer="test", driver='GPKG')
>>>>>   sf2a <- st_read(dsn=fn1,layer='test')
>>>>>   st_coordinates(st_geometry(sf2a))
>>>>>   SPDF2a <- readOGR(dsn=fn1)
>>>>>   sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "hole")
>>>>>   sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>>>>   rgeos::writeWKT(SPDF2a)
>>>>>   st_as_text(st_geometry(st_as_sf(SPDF2a)))
>>>>>
>>>>>   # the issue is resolved. If we separate the exterior rings:
>>>>>
>>>>>   r2a <- r1; r2a[,1] <- r2a[,1]+1.00001
>>>>>   Ps1a = Polygons(list(Polygon(r1),Polygon(r2a)),ID=1)
>>>>>   SPDFa = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1a)),
>>>>>   data.frame(Example=c("Minimal")))
>>>>>   fna <- tempfile()
>>>>>   writeOGR(SPDFa, fna, layer='test', driver='ESRI Shapefile')
>>>>>   SPDF2_a <- readOGR(dsn=fna, layer='test')
>>>>>   sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot, "hole")
>>>>>   sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot,
>>>>>   "ringDir")
>>>>>   rgeos::writeWKT(SPDF2_a)
>>>>>   st_as_text(st_geometry(st_as_sf(SPDF2_a)))
>>>>>
>>>>>   # we are OK as the two exterior rings do not touch.
>>>>>
>>>>>   # does using sf make a difference?
>>>>>
>>>>>   fn_s <- tempfile(fileext=".shp")
>>>>>   st_write(st_as_sf(SPDF), dsn=fn_s)
>>>>>   sf_in <- st_read(fn_s)
>>>>>   st_as_text(st_geometry(st_as_sf(sf_in)))
>>>>>
>>>>>   # No
>>>>>
>>>>>   fn_s <- tempfile(fileext=".shp")
>>>>>   st_write(st_as_sf(SPDF_c), dsn=fn_s)
>>>>>   sf_in_c <- st_read(fn_s)
>>>>>   st_as_text(st_geometry(st_as_sf(sf_in_c)))
>>>>>
>>>>>   # nor with the pretend-SF-compliant comment set either.
>>>>>
>>>>>   So the weakness is in the "ESRI Shapefile" write driver, or possibly in
>>>>>   the OGRGeometryFactory::organizePolygons() function in GDAL used in
>>>>>   OGR_write() (a C++ function) called by writeOGR(). If sf::st_write()
>>>>>   also
>>>>>   calls OGRGeometryFactory::organizePolygons(), we'd maybe consider that
>>>>>   it
>>>>>   has a weakness for the "ESRI Shapefile" driver, but which does not
>>>>>   affect
>>>>>   SF-compliant drivers.
>>>>
>>>>  Without the comment set, OGRGeometryFactory::organizePolygons() is used;
>>>>  with it set, OGRGeometryFactory::organizePolygons() is not used, because
>>>>  the object is declared as two exterior rings. In both cases, we have the
>>>>  output object written out and read back in incorrectly with the ESRI
>>>>  shapefile driver, but SF-compliant drivers round-trip (in the test
>>>>  GeoJSON), correctly.
>>>>
>>>>  It is likely that the changes made in 2015 to accommodate GeoJSON led to
>>>>  this possible regression for the ESRI Shapefile driver. I'm adding this
>>>>  geometry to tests/tripup.R and data files; without code changes the hole
>>>>  slot is wrong and the ring direction is changes to match, so the
>>>>  coordinates change too.
>>>>
>>>>  Reading using the deprecated maptools::readShapeSpatial() also gets a hole
>>>>  in the second external ring. However, writing with deprecated
>>>> maptools::  writeSpatialShape() yields a shapefile that when read with
>>>> maptools::  readShapeSpatial() gives the correct two exterior ring, no hole
>>>>  object. When read with sf::st_read() and rgdal::readOGR(), the object is
>>>>  also correct. So the problem definitely lies in rgdal::writeOGR(), and
>>>>  sf::st_write() - roundtripping with sf::st_write() and sf::st_read()
>>>>  degrades from MULTIPOLYGON to POLYGON with the ESRI Shapefile driver.
>>>>
>>>>  Roger
>>>>
>>>>>
>>>>>   This is probably related to a similar but inverse problem with the
>>>>>   SF-compliant GeoJSON driver in 2015:
>>>>>
>>>>>   https://stat.ethz.ch/pipermail/r-sig-geo/2015-October/023609.html
>>>>>
>>>>>   continued the next month in:
>>>>>
>>>>>   https://stat.ethz.ch/pipermail/r-sig-geo/2015-November/023656.html
>>>>>
>>>>>   The details are in this SVN diff
>>>>>
>>>>>   https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=555&r2=571
>>>>>
>>>>>   up to the end og the list thread, and this one from then until now:
>>>>>
>>>>>   https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=571&r2=733
>>>>>
>>>>>   Summary: could you change drivers, or is it really necessary to fix an
>>>>>   EOL
>>>>>   problem? What is your use case?
>>>>>
>>>>>>
>>>>>>    Thanks in advance for your time,
>>>>>
>>>>>   Thanks for a complete example,
>>>>>
>>>>>   Roger
>>>>>
>>>>>>    Phil
>>>>>>
>>>>>>>    sessionInfo()
>>>>>>    R version 3.5.1 (2018-07-02)
>>>>>>    Platform: x86_64-w64-mingw32/x64 (64-bit)
>>>>>>    Running under: Windows >= 8 x64 (build 9200)
>>>>>>
>>>>>>    Matrix products: default
>>>>>>
>>>>>>    locale:
>>>>>>    [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United
>>>>>>    Kingdom.1252    LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
>>>>>>                            LC_TIME=English_United Kingdom.1252
>>>>>>
>>>>>>    attached base packages:
>>>>>>    [1] stats     graphics  grDevices utils     datasets  methods   base
>>>>>>
>>>>>>    other attached packages:
>>>>>>    [1] rgdal_1.4-4 sp_1.3-1
>>>>>>
>>>>>>    loaded via a namespace (and not attached):
>>>>>>    [1] compiler_3.5.1  tools_3.5.1     yaml_2.2.0      grid_3.5.1
>>>>>>    lattice_0.20-35
>>>>>>
>>>>>>    _______________________________________________
>>>>>>    R-sig-Geo mailing list
>>>>>>    [hidden email]
>>>>>>    https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>
>>>
>>
>> --
>> Roger Bivand
>> Department of Economics, Norwegian School of Economics,
>> Helleveien 30, N-5045 Bergen, Norway.
>> voice: +47 55 95 93 55; e-mail: [hidden email]
>> https://orcid.org/0000-0003-2392-6140
>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Re: rgdal::writeOGR with driver='ESRI Shapefile' converts Polygon object into a hole

Sun, 08/18/2019 - 15:33
Hi Roger,

I originally encountered this issue while trying to reduce the size of
German postal code boundaries. I used rmapshaper::ms_simplify() which
introduced the polygons sharing a common edge. I definitely didn't
need them to be separate polygons, and would happily have merged them
had I been more familiar with gUnaryUnion().

I apologise if my question has resulted in unnecessary effort on your
part. I reported the behaviour because I found it surprising. However,
I now understand that this example is not a valid simple feature
geometry and that the solution is to take steps to ensure I have valid
geometries, and to repair if not.

Additionally I must confess that I only use the ESRI Shapefile driver
because I don't know any better... My use cases are reading and
writing from R (usually to produce maps) and occasionally opening in
QGIS. I can happily switch to any format that supports this workflow.

Thank you very much for your time,
Phil

On Sun, 18 Aug 2019 at 13:59, Roger Bivand <[hidden email]> wrote:
>
> The reason that the problem occurred is that a MULTIPOLYGON with two
> exterior rings becomes invalid if the exterior rings touch along an edge
> (this case). It is important to know the use case, to see whether:
>
> library(rgeos)
> writeWKT(Ps1)
> gIsValid(Ps1, reason=TRUE)
> Ps1a <- gUnaryUnion(gBuffer(Ps1, width=0))
> gIsValid(Ps1a, reason=TRUE)
> writeWKT(Ps1a)
>
> or equivalently in sf for sf objects, should be applied before trying to
> write the object out to file with this driver.
>
> However, because drivers that are compliant with the simple features
> standard (which bans exterior rings sharing edges) have been permissive
> and do round-trip this invalid object, a relaxation in the OGR ESRI
> shapefile driver has been provided and may be included in a future
> release.
>
> We need to know (see these issues):
>
> https://github.com/r-spatial/sf/issues/1130
> https://github.com/OSGeo/gdal/issues/1787
>
> why it was desirable to write out this object using this driver? Could an
> alternative driver have been used, or is ESRI shapefile the only format
> used in the workflow?
>
> If it has to be this driver, could the workflow be changed to repair
> degenerate cases before writing? If using sp classes, rgeos may be used to
> test for and probably repair such geometries before they reach
> rgdal::writeOGR() for this driver. Adding code to sf and rgdal to trap
> degenerate cases does encumber all users with valid geometries with
> the time wasted on extra checking, so building checks into sf and rgdal is
> not desirable.
>
> I hope it is possible to find out more about the use case quickly, to pass
> on to GDAL developers to help motivate a relaxation in their current
> policy with regard to this driver, and to encourage them to include the
> fix branch in a future release of GDAL.
>
> Roger
>
> On Sat, 17 Aug 2019, Roger Bivand wrote:
>
> > Please follow up both here and on:
> >
> > https://github.com/r-spatial/sf/issues/1130
> >
> > as the problem is also seen in the sf package using the same GDAL ESRI
> > Shapefile driver.
> >
> > Roger
> >
> > On Fri, 16 Aug 2019, Roger Bivand wrote:
> >
> >>  On Fri, 16 Aug 2019, Roger Bivand wrote:
> >>
> >>>   On Tue, 13 Aug 2019, Phil Haines wrote:
> >>>
> >>>>    Dear list,
> >>>>
> >>>>    I have a single Polygons object containing multiple Polygon objects
> >>>>    that share a common border. When I output this using writeOGR() one of
> >>>>    the Polygon objects becomes a hole, as the following example shows.
> >>>>
> >>>>    Create a Polygons object containing two adjoining Polygon objects
> >>>>>    library(rgdal)
> >>>>>    r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
> >>>>>    r2 <- r1; r2[,1] <- r2[,1]+1
> >>>>>    Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
> >>>>>    SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
> >>>>>    data.frame(Example=c("Minimal")))
> >>>>
> >>>>    Perform a write/readOGR() cycle
> >>>>>    fn <- tempfile()
> >>>>>    writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
> >>>>>    SPDF2 <- readOGR(dsn=fn,layer='test')
> >>>>
> >>>>    Second Polygon object is now a hole
> >>>>>    sapply(SPDF2@polygons[[1]]@Polygons,slot,"hole")
> >>>>    [1] FALSE  TRUE
> >>>>
> >>>>    I see from the sp documentation that "Polygon objects belonging to a
> >>>>    Polygons object should either not overlap one-other, or should be
> >>>>    fully included" but I am not sure how this relates to bordering
> >>>>    Polygon objects. I would welcome any advice as to whether what I am
> >>>>    asking of writeOGR is reasonable?
> >>>
> >>>   The problem is with the 'ESRI Shapefile' representation and driver:
> >>>
> >>>   library(rgdal)
> >>>   r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
> >>>   r2 <- r1; r2[,1] <- r2[,1]+1
> >>>   Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
> >>>   SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
> >>>   data.frame(Example=c("Minimal")))
> >>>   sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "hole")
> >>>   sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "ringDir")
> >>>
> >>>   # which constructs a MULTIPOLYGON object:
> >>>
> >>>   rgeos::writeWKT(SPDF)
> >>>   library(sf)
> >>>   st_as_text(st_geometry(st_as_sf(SPDF)))
> >>>
> >>> #    The 'ESRI Shapefile' driver is not Simple-Feature compliant (it
> >>> #    pre-dates it), so the failure occurs by seeing the second exterior
> >>> #    ring
> >>> #    as an interior ring
> >>>
> >>>   fn <- tempfile()
> >>>   writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
> >>>   SPDF2 <- readOGR(dsn=fn, layer='test')
> >>>   sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "hole")
> >>>   sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "ringDir")
> >>>   rgeos::writeWKT(SPDF2)
> >>>   st_as_text(st_geometry(st_as_sf(SPDF2)))
> >>>
> >>>   # This happens with sf too, using the same GDAL driver:
> >>>
> >>>   sf2 <- st_read(dsn=fn, layer='test')
> >>>   st_geometry(sf2)
> >>>   library(sf)
> >>>   st_as_text(st_geometry(st_as_sf(sf2)))
> >>>   rgeos::writeWKT(as(sf2, "Spatial"))
> >>>
> >>>   # Adding the comment fix doesn't help:
> >>>
> >>>   comment(slot(SPDF, "polygons")[[1]])
> >>>   SPDF_c <- rgeos::createSPComment(SPDF)
> >>>   comment(slot(SPDF_c, "polygons")[[1]])
> >>>   writeOGR(SPDF_c, fn, layer='test_c', driver='ESRI Shapefile',
> >>>     verbose=TRUE)
> >>>
> >>> #    reports
> >>> #    Object initially classed as: wkbPolygon
> >>> #    SFS comments in Polygons objects
> >>> #    Object reclassed as: wkbMultiPolygon
> >>>
> >>>   SPDF2_c <- readOGR(dsn=fn, layer='test_c')
> >>>   sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot, "hole")
> >>>   sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot,
> >>>   "ringDir")
> >>>   rgeos::writeWKT(SPDF2_c)
> >>>   st_as_text(st_geometry(st_as_sf(SPDF2_c)))
> >>>
> >>>
> >>>
> >>>   # If the input object is written out with the GeoPackage driver:
> >>>
> >>>   fn1 <- tempfile(fileext=".gpkg")
> >>>   writeOGR(SPDF, fn1, layer="test", driver='GPKG')
> >>>   sf2a <- st_read(dsn=fn1,layer='test')
> >>>   st_coordinates(st_geometry(sf2a))
> >>>   SPDF2a <- readOGR(dsn=fn1)
> >>>   sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "hole")
> >>>   sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "ringDir")
> >>>   rgeos::writeWKT(SPDF2a)
> >>>   st_as_text(st_geometry(st_as_sf(SPDF2a)))
> >>>
> >>>   # the issue is resolved. If we separate the exterior rings:
> >>>
> >>>   r2a <- r1; r2a[,1] <- r2a[,1]+1.00001
> >>>   Ps1a = Polygons(list(Polygon(r1),Polygon(r2a)),ID=1)
> >>>   SPDFa = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1a)),
> >>>   data.frame(Example=c("Minimal")))
> >>>   fna <- tempfile()
> >>>   writeOGR(SPDFa, fna, layer='test', driver='ESRI Shapefile')
> >>>   SPDF2_a <- readOGR(dsn=fna, layer='test')
> >>>   sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot, "hole")
> >>>   sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot,
> >>>   "ringDir")
> >>>   rgeos::writeWKT(SPDF2_a)
> >>>   st_as_text(st_geometry(st_as_sf(SPDF2_a)))
> >>>
> >>>   # we are OK as the two exterior rings do not touch.
> >>>
> >>>   # does using sf make a difference?
> >>>
> >>>   fn_s <- tempfile(fileext=".shp")
> >>>   st_write(st_as_sf(SPDF), dsn=fn_s)
> >>>   sf_in <- st_read(fn_s)
> >>>   st_as_text(st_geometry(st_as_sf(sf_in)))
> >>>
> >>>   # No
> >>>
> >>>   fn_s <- tempfile(fileext=".shp")
> >>>   st_write(st_as_sf(SPDF_c), dsn=fn_s)
> >>>   sf_in_c <- st_read(fn_s)
> >>>   st_as_text(st_geometry(st_as_sf(sf_in_c)))
> >>>
> >>>   # nor with the pretend-SF-compliant comment set either.
> >>>
> >>>   So the weakness is in the "ESRI Shapefile" write driver, or possibly in
> >>>   the OGRGeometryFactory::organizePolygons() function in GDAL used in
> >>>   OGR_write() (a C++ function) called by writeOGR(). If sf::st_write()
> >>>   also
> >>>   calls OGRGeometryFactory::organizePolygons(), we'd maybe consider that
> >>>   it
> >>>   has a weakness for the "ESRI Shapefile" driver, but which does not
> >>>   affect
> >>>   SF-compliant drivers.
> >>
> >>  Without the comment set, OGRGeometryFactory::organizePolygons() is used;
> >>  with it set, OGRGeometryFactory::organizePolygons() is not used, because
> >>  the object is declared as two exterior rings. In both cases, we have the
> >>  output object written out and read back in incorrectly with the ESRI
> >>  shapefile driver, but SF-compliant drivers round-trip (in the test
> >>  GeoJSON), correctly.
> >>
> >>  It is likely that the changes made in 2015 to accommodate GeoJSON led to
> >>  this possible regression for the ESRI Shapefile driver. I'm adding this
> >>  geometry to tests/tripup.R and data files; without code changes the hole
> >>  slot is wrong and the ring direction is changes to match, so the
> >>  coordinates change too.
> >>
> >>  Reading using the deprecated maptools::readShapeSpatial() also gets a hole
> >>  in the second external ring. However, writing with deprecated
> >> maptools::  writeSpatialShape() yields a shapefile that when read with
> >> maptools::  readShapeSpatial() gives the correct two exterior ring, no hole
> >>  object. When read with sf::st_read() and rgdal::readOGR(), the object is
> >>  also correct. So the problem definitely lies in rgdal::writeOGR(), and
> >>  sf::st_write() - roundtripping with sf::st_write() and sf::st_read()
> >>  degrades from MULTIPOLYGON to POLYGON with the ESRI Shapefile driver.
> >>
> >>  Roger
> >>
> >>>
> >>>   This is probably related to a similar but inverse problem with the
> >>>   SF-compliant GeoJSON driver in 2015:
> >>>
> >>>   https://stat.ethz.ch/pipermail/r-sig-geo/2015-October/023609.html
> >>>
> >>>   continued the next month in:
> >>>
> >>>   https://stat.ethz.ch/pipermail/r-sig-geo/2015-November/023656.html
> >>>
> >>>   The details are in this SVN diff
> >>>
> >>>   https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=555&r2=571
> >>>
> >>>   up to the end og the list thread, and this one from then until now:
> >>>
> >>>   https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=571&r2=733
> >>>
> >>>   Summary: could you change drivers, or is it really necessary to fix an
> >>>   EOL
> >>>   problem? What is your use case?
> >>>
> >>>>
> >>>>    Thanks in advance for your time,
> >>>
> >>>   Thanks for a complete example,
> >>>
> >>>   Roger
> >>>
> >>>>    Phil
> >>>>
> >>>>>    sessionInfo()
> >>>>    R version 3.5.1 (2018-07-02)
> >>>>    Platform: x86_64-w64-mingw32/x64 (64-bit)
> >>>>    Running under: Windows >= 8 x64 (build 9200)
> >>>>
> >>>>    Matrix products: default
> >>>>
> >>>>    locale:
> >>>>    [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United
> >>>>    Kingdom.1252    LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
> >>>>                            LC_TIME=English_United Kingdom.1252
> >>>>
> >>>>    attached base packages:
> >>>>    [1] stats     graphics  grDevices utils     datasets  methods   base
> >>>>
> >>>>    other attached packages:
> >>>>    [1] rgdal_1.4-4 sp_1.3-1
> >>>>
> >>>>    loaded via a namespace (and not attached):
> >>>>    [1] compiler_3.5.1  tools_3.5.1     yaml_2.2.0      grid_3.5.1
> >>>>    lattice_0.20-35
> >>>>
> >>>>    _______________________________________________
> >>>>    R-sig-Geo mailing list
> >>>>    [hidden email]
> >>>>    https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>>>
> >>>
> >>>
> >>
> >>
> >
> >
>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: [hidden email]
> https://orcid.org/0000-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: rgdal::writeOGR with driver='ESRI Shapefile' converts Polygon object into a hole

Sun, 08/18/2019 - 07:59
The reason that the problem occurred is that a MULTIPOLYGON with two
exterior rings becomes invalid if the exterior rings touch along an edge
(this case). It is important to know the use case, to see whether:

library(rgeos)
writeWKT(Ps1)
gIsValid(Ps1, reason=TRUE)
Ps1a <- gUnaryUnion(gBuffer(Ps1, width=0))
gIsValid(Ps1a, reason=TRUE)
writeWKT(Ps1a)

or equivalently in sf for sf objects, should be applied before trying to
write the object out to file with this driver.

However, because drivers that are compliant with the simple features
standard (which bans exterior rings sharing edges) have been permissive
and do round-trip this invalid object, a relaxation in the OGR ESRI
shapefile driver has been provided and may be included in a future
release.

We need to know (see these issues):

https://github.com/r-spatial/sf/issues/1130
https://github.com/OSGeo/gdal/issues/1787

why it was desirable to write out this object using this driver? Could an
alternative driver have been used, or is ESRI shapefile the only format
used in the workflow?

If it has to be this driver, could the workflow be changed to repair
degenerate cases before writing? If using sp classes, rgeos may be used to
test for and probably repair such geometries before they reach
rgdal::writeOGR() for this driver. Adding code to sf and rgdal to trap
degenerate cases does encumber all users with valid geometries with
the time wasted on extra checking, so building checks into sf and rgdal is
not desirable.

I hope it is possible to find out more about the use case quickly, to pass
on to GDAL developers to help motivate a relaxation in their current
policy with regard to this driver, and to encourage them to include the
fix branch in a future release of GDAL.

Roger

On Sat, 17 Aug 2019, Roger Bivand wrote:

> Please follow up both here and on:
>
> https://github.com/r-spatial/sf/issues/1130
>
> as the problem is also seen in the sf package using the same GDAL ESRI
> Shapefile driver.
>
> Roger
>
> On Fri, 16 Aug 2019, Roger Bivand wrote:
>
>>  On Fri, 16 Aug 2019, Roger Bivand wrote:
>>
>>>   On Tue, 13 Aug 2019, Phil Haines wrote:
>>>
>>>>    Dear list,
>>>>
>>>>    I have a single Polygons object containing multiple Polygon objects
>>>>    that share a common border. When I output this using writeOGR() one of
>>>>    the Polygon objects becomes a hole, as the following example shows.
>>>>
>>>>    Create a Polygons object containing two adjoining Polygon objects
>>>>>    library(rgdal)
>>>>>    r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
>>>>>    r2 <- r1; r2[,1] <- r2[,1]+1
>>>>>    Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
>>>>>    SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
>>>>>    data.frame(Example=c("Minimal")))
>>>>
>>>>    Perform a write/readOGR() cycle
>>>>>    fn <- tempfile()
>>>>>    writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
>>>>>    SPDF2 <- readOGR(dsn=fn,layer='test')
>>>>
>>>>    Second Polygon object is now a hole
>>>>>    sapply(SPDF2@polygons[[1]]@Polygons,slot,"hole")
>>>>    [1] FALSE  TRUE
>>>>
>>>>    I see from the sp documentation that "Polygon objects belonging to a
>>>>    Polygons object should either not overlap one-other, or should be
>>>>    fully included" but I am not sure how this relates to bordering
>>>>    Polygon objects. I would welcome any advice as to whether what I am
>>>>    asking of writeOGR is reasonable?
>>>
>>>   The problem is with the 'ESRI Shapefile' representation and driver:
>>>
>>>   library(rgdal)
>>>   r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
>>>   r2 <- r1; r2[,1] <- r2[,1]+1
>>>   Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
>>>   SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
>>>   data.frame(Example=c("Minimal")))
>>>   sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "hole")
>>>   sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>>
>>>   # which constructs a MULTIPOLYGON object:
>>>
>>>   rgeos::writeWKT(SPDF)
>>>   library(sf)
>>>   st_as_text(st_geometry(st_as_sf(SPDF)))
>>>
>>> #    The 'ESRI Shapefile' driver is not Simple-Feature compliant (it
>>> #    pre-dates it), so the failure occurs by seeing the second exterior
>>> #    ring
>>> #    as an interior ring
>>>
>>>   fn <- tempfile()
>>>   writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
>>>   SPDF2 <- readOGR(dsn=fn, layer='test')
>>>   sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "hole")
>>>   sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>>   rgeos::writeWKT(SPDF2)
>>>   st_as_text(st_geometry(st_as_sf(SPDF2)))
>>>
>>>   # This happens with sf too, using the same GDAL driver:
>>>
>>>   sf2 <- st_read(dsn=fn, layer='test')
>>>   st_geometry(sf2)
>>>   library(sf)
>>>   st_as_text(st_geometry(st_as_sf(sf2)))
>>>   rgeos::writeWKT(as(sf2, "Spatial"))
>>>
>>>   # Adding the comment fix doesn't help:
>>>
>>>   comment(slot(SPDF, "polygons")[[1]])
>>>   SPDF_c <- rgeos::createSPComment(SPDF)
>>>   comment(slot(SPDF_c, "polygons")[[1]])
>>>   writeOGR(SPDF_c, fn, layer='test_c', driver='ESRI Shapefile',
>>>     verbose=TRUE)
>>>
>>> #    reports
>>> #    Object initially classed as: wkbPolygon
>>> #    SFS comments in Polygons objects
>>> #    Object reclassed as: wkbMultiPolygon
>>>
>>>   SPDF2_c <- readOGR(dsn=fn, layer='test_c')
>>>   sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot, "hole")
>>>   sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot,
>>>   "ringDir")
>>>   rgeos::writeWKT(SPDF2_c)
>>>   st_as_text(st_geometry(st_as_sf(SPDF2_c)))
>>>
>>>
>>>
>>>   # If the input object is written out with the GeoPackage driver:
>>>
>>>   fn1 <- tempfile(fileext=".gpkg")
>>>   writeOGR(SPDF, fn1, layer="test", driver='GPKG')
>>>   sf2a <- st_read(dsn=fn1,layer='test')
>>>   st_coordinates(st_geometry(sf2a))
>>>   SPDF2a <- readOGR(dsn=fn1)
>>>   sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "hole")
>>>   sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>>   rgeos::writeWKT(SPDF2a)
>>>   st_as_text(st_geometry(st_as_sf(SPDF2a)))
>>>
>>>   # the issue is resolved. If we separate the exterior rings:
>>>
>>>   r2a <- r1; r2a[,1] <- r2a[,1]+1.00001
>>>   Ps1a = Polygons(list(Polygon(r1),Polygon(r2a)),ID=1)
>>>   SPDFa = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1a)),
>>>   data.frame(Example=c("Minimal")))
>>>   fna <- tempfile()
>>>   writeOGR(SPDFa, fna, layer='test', driver='ESRI Shapefile')
>>>   SPDF2_a <- readOGR(dsn=fna, layer='test')
>>>   sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot, "hole")
>>>   sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot,
>>>   "ringDir")
>>>   rgeos::writeWKT(SPDF2_a)
>>>   st_as_text(st_geometry(st_as_sf(SPDF2_a)))
>>>
>>>   # we are OK as the two exterior rings do not touch.
>>>
>>>   # does using sf make a difference?
>>>
>>>   fn_s <- tempfile(fileext=".shp")
>>>   st_write(st_as_sf(SPDF), dsn=fn_s)
>>>   sf_in <- st_read(fn_s)
>>>   st_as_text(st_geometry(st_as_sf(sf_in)))
>>>
>>>   # No
>>>
>>>   fn_s <- tempfile(fileext=".shp")
>>>   st_write(st_as_sf(SPDF_c), dsn=fn_s)
>>>   sf_in_c <- st_read(fn_s)
>>>   st_as_text(st_geometry(st_as_sf(sf_in_c)))
>>>
>>>   # nor with the pretend-SF-compliant comment set either.
>>>
>>>   So the weakness is in the "ESRI Shapefile" write driver, or possibly in
>>>   the OGRGeometryFactory::organizePolygons() function in GDAL used in
>>>   OGR_write() (a C++ function) called by writeOGR(). If sf::st_write()
>>>   also
>>>   calls OGRGeometryFactory::organizePolygons(), we'd maybe consider that
>>>   it
>>>   has a weakness for the "ESRI Shapefile" driver, but which does not
>>>   affect
>>>   SF-compliant drivers.
>>
>>  Without the comment set, OGRGeometryFactory::organizePolygons() is used;
>>  with it set, OGRGeometryFactory::organizePolygons() is not used, because
>>  the object is declared as two exterior rings. In both cases, we have the
>>  output object written out and read back in incorrectly with the ESRI
>>  shapefile driver, but SF-compliant drivers round-trip (in the test
>>  GeoJSON), correctly.
>>
>>  It is likely that the changes made in 2015 to accommodate GeoJSON led to
>>  this possible regression for the ESRI Shapefile driver. I'm adding this
>>  geometry to tests/tripup.R and data files; without code changes the hole
>>  slot is wrong and the ring direction is changes to match, so the
>>  coordinates change too.
>>
>>  Reading using the deprecated maptools::readShapeSpatial() also gets a hole
>>  in the second external ring. However, writing with deprecated
>> maptools::  writeSpatialShape() yields a shapefile that when read with
>> maptools::  readShapeSpatial() gives the correct two exterior ring, no hole
>>  object. When read with sf::st_read() and rgdal::readOGR(), the object is
>>  also correct. So the problem definitely lies in rgdal::writeOGR(), and
>>  sf::st_write() - roundtripping with sf::st_write() and sf::st_read()
>>  degrades from MULTIPOLYGON to POLYGON with the ESRI Shapefile driver.
>>
>>  Roger
>>
>>>
>>>   This is probably related to a similar but inverse problem with the
>>>   SF-compliant GeoJSON driver in 2015:
>>>
>>>   https://stat.ethz.ch/pipermail/r-sig-geo/2015-October/023609.html
>>>
>>>   continued the next month in:
>>>
>>>   https://stat.ethz.ch/pipermail/r-sig-geo/2015-November/023656.html
>>>
>>>   The details are in this SVN diff
>>>
>>>   https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=555&r2=571
>>>
>>>   up to the end og the list thread, and this one from then until now:
>>>
>>>   https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=571&r2=733
>>>
>>>   Summary: could you change drivers, or is it really necessary to fix an
>>>   EOL
>>>   problem? What is your use case?
>>>
>>>>
>>>>    Thanks in advance for your time,
>>>
>>>   Thanks for a complete example,
>>>
>>>   Roger
>>>
>>>>    Phil
>>>>
>>>>>    sessionInfo()
>>>>    R version 3.5.1 (2018-07-02)
>>>>    Platform: x86_64-w64-mingw32/x64 (64-bit)
>>>>    Running under: Windows >= 8 x64 (build 9200)
>>>>
>>>>    Matrix products: default
>>>>
>>>>    locale:
>>>>    [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United
>>>>    Kingdom.1252    LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
>>>>                            LC_TIME=English_United Kingdom.1252
>>>>
>>>>    attached base packages:
>>>>    [1] stats     graphics  grDevices utils     datasets  methods   base
>>>>
>>>>    other attached packages:
>>>>    [1] rgdal_1.4-4 sp_1.3-1
>>>>
>>>>    loaded via a namespace (and not attached):
>>>>    [1] compiler_3.5.1  tools_3.5.1     yaml_2.2.0      grid_3.5.1
>>>>    lattice_0.20-35
>>>>
>>>>    _______________________________________________
>>>>    R-sig-Geo mailing list
>>>>    [hidden email]
>>>>    https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>>
>>
>>
>
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Re: rgdal::writeOGR with driver='ESRI Shapefile' converts Polygon object into a hole

Sat, 08/17/2019 - 04:06
Please follow up both here and on:

https://github.com/r-spatial/sf/issues/1130

as the problem is also seen in the sf package using the same GDAL ESRI
Shapefile driver.

Roger

On Fri, 16 Aug 2019, Roger Bivand wrote:

> On Fri, 16 Aug 2019, Roger Bivand wrote:
>
>>  On Tue, 13 Aug 2019, Phil Haines wrote:
>>
>>>   Dear list,
>>>
>>>   I have a single Polygons object containing multiple Polygon objects
>>>   that share a common border. When I output this using writeOGR() one of
>>>   the Polygon objects becomes a hole, as the following example shows.
>>>
>>>   Create a Polygons object containing two adjoining Polygon objects
>>>>   library(rgdal)
>>>>   r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
>>>>   r2 <- r1; r2[,1] <- r2[,1]+1
>>>>   Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
>>>>   SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
>>>>   data.frame(Example=c("Minimal")))
>>>
>>>   Perform a write/readOGR() cycle
>>>>   fn <- tempfile()
>>>>   writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
>>>>   SPDF2 <- readOGR(dsn=fn,layer='test')
>>>
>>>   Second Polygon object is now a hole
>>>>   sapply(SPDF2@polygons[[1]]@Polygons,slot,"hole")
>>>   [1] FALSE  TRUE
>>>
>>>   I see from the sp documentation that "Polygon objects belonging to a
>>>   Polygons object should either not overlap one-other, or should be
>>>   fully included" but I am not sure how this relates to bordering
>>>   Polygon objects. I would welcome any advice as to whether what I am
>>>   asking of writeOGR is reasonable?
>>
>>  The problem is with the 'ESRI Shapefile' representation and driver:
>>
>>  library(rgdal)
>>  r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
>>  r2 <- r1; r2[,1] <- r2[,1]+1
>>  Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
>>  SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
>>  data.frame(Example=c("Minimal")))
>>  sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "hole")
>>  sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>
>>  # which constructs a MULTIPOLYGON object:
>>
>>  rgeos::writeWKT(SPDF)
>>  library(sf)
>>  st_as_text(st_geometry(st_as_sf(SPDF)))
>>
>> #   The 'ESRI Shapefile' driver is not Simple-Feature compliant (it
>> #   pre-dates it), so the failure occurs by seeing the second exterior ring
>> #   as an interior ring
>>
>>  fn <- tempfile()
>>  writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
>>  SPDF2 <- readOGR(dsn=fn, layer='test')
>>  sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "hole")
>>  sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>  rgeos::writeWKT(SPDF2)
>>  st_as_text(st_geometry(st_as_sf(SPDF2)))
>>
>>  # This happens with sf too, using the same GDAL driver:
>>
>>  sf2 <- st_read(dsn=fn, layer='test')
>>  st_geometry(sf2)
>>  library(sf)
>>  st_as_text(st_geometry(st_as_sf(sf2)))
>>  rgeos::writeWKT(as(sf2, "Spatial"))
>>
>>  # Adding the comment fix doesn't help:
>>
>>  comment(slot(SPDF, "polygons")[[1]])
>>  SPDF_c <- rgeos::createSPComment(SPDF)
>>  comment(slot(SPDF_c, "polygons")[[1]])
>>  writeOGR(SPDF_c, fn, layer='test_c', driver='ESRI Shapefile',
>>    verbose=TRUE)
>>
>> #   reports
>> #   Object initially classed as: wkbPolygon
>> #   SFS comments in Polygons objects
>> #   Object reclassed as: wkbMultiPolygon
>>
>>  SPDF2_c <- readOGR(dsn=fn, layer='test_c')
>>  sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot, "hole")
>>  sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>  rgeos::writeWKT(SPDF2_c)
>>  st_as_text(st_geometry(st_as_sf(SPDF2_c)))
>>
>>
>>
>>  # If the input object is written out with the GeoPackage driver:
>>
>>  fn1 <- tempfile(fileext=".gpkg")
>>  writeOGR(SPDF, fn1, layer="test", driver='GPKG')
>>  sf2a <- st_read(dsn=fn1,layer='test')
>>  st_coordinates(st_geometry(sf2a))
>>  SPDF2a <- readOGR(dsn=fn1)
>>  sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "hole")
>>  sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>  rgeos::writeWKT(SPDF2a)
>>  st_as_text(st_geometry(st_as_sf(SPDF2a)))
>>
>>  # the issue is resolved. If we separate the exterior rings:
>>
>>  r2a <- r1; r2a[,1] <- r2a[,1]+1.00001
>>  Ps1a = Polygons(list(Polygon(r1),Polygon(r2a)),ID=1)
>>  SPDFa = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1a)),
>>  data.frame(Example=c("Minimal")))
>>  fna <- tempfile()
>>  writeOGR(SPDFa, fna, layer='test', driver='ESRI Shapefile')
>>  SPDF2_a <- readOGR(dsn=fna, layer='test')
>>  sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot, "hole")
>>  sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot, "ringDir")
>>  rgeos::writeWKT(SPDF2_a)
>>  st_as_text(st_geometry(st_as_sf(SPDF2_a)))
>>
>>  # we are OK as the two exterior rings do not touch.
>>
>>  # does using sf make a difference?
>>
>>  fn_s <- tempfile(fileext=".shp")
>>  st_write(st_as_sf(SPDF), dsn=fn_s)
>>  sf_in <- st_read(fn_s)
>>  st_as_text(st_geometry(st_as_sf(sf_in)))
>>
>>  # No
>>
>>  fn_s <- tempfile(fileext=".shp")
>>  st_write(st_as_sf(SPDF_c), dsn=fn_s)
>>  sf_in_c <- st_read(fn_s)
>>  st_as_text(st_geometry(st_as_sf(sf_in_c)))
>>
>>  # nor with the pretend-SF-compliant comment set either.
>>
>>  So the weakness is in the "ESRI Shapefile" write driver, or possibly in
>>  the OGRGeometryFactory::organizePolygons() function in GDAL used in
>>  OGR_write() (a C++ function) called by writeOGR(). If sf::st_write() also
>>  calls OGRGeometryFactory::organizePolygons(), we'd maybe consider that it
>>  has a weakness for the "ESRI Shapefile" driver, but which does not affect
>>  SF-compliant drivers.
>
> Without the comment set, OGRGeometryFactory::organizePolygons() is used; with
> it set, OGRGeometryFactory::organizePolygons() is not used, because the
> object is declared as two exterior rings. In both cases, we have the output
> object written out and read back in incorrectly with the ESRI shapefile
> driver, but SF-compliant drivers round-trip (in the test GeoJSON), correctly.
>
> It is likely that the changes made in 2015 to accommodate GeoJSON led to this
> possible regression for the ESRI Shapefile driver. I'm adding this geometry
> to tests/tripup.R and data files; without code changes the hole slot is wrong
> and the ring direction is changes to match, so the coordinates change too.
>
> Reading using the deprecated maptools::readShapeSpatial() also gets a hole in
> the second external ring. However, writing with deprecated
> maptools:: writeSpatialShape() yields a shapefile that when read with
> maptools:: readShapeSpatial() gives the correct two exterior ring, no hole
> object. When read with sf::st_read() and rgdal::readOGR(), the object is also
> correct. So the problem definitely lies in rgdal::writeOGR(), and
> sf::st_write() - roundtripping with sf::st_write() and sf::st_read() degrades
> from MULTIPOLYGON to POLYGON with the ESRI Shapefile driver.
>
> Roger
>
>>
>>  This is probably related to a similar but inverse problem with the
>>  SF-compliant GeoJSON driver in 2015:
>>
>>  https://stat.ethz.ch/pipermail/r-sig-geo/2015-October/023609.html
>>
>>  continued the next month in:
>>
>>  https://stat.ethz.ch/pipermail/r-sig-geo/2015-November/023656.html
>>
>>  The details are in this SVN diff
>>
>>  https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=555&r2=571
>>
>>  up to the end og the list thread, and this one from then until now:
>>
>>  https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=571&r2=733
>>
>>  Summary: could you change drivers, or is it really necessary to fix an EOL
>>  problem? What is your use case?
>>
>>>
>>>   Thanks in advance for your time,
>>
>>  Thanks for a complete example,
>>
>>  Roger
>>
>>>   Phil
>>>
>>>>   sessionInfo()
>>>   R version 3.5.1 (2018-07-02)
>>>   Platform: x86_64-w64-mingw32/x64 (64-bit)
>>>   Running under: Windows >= 8 x64 (build 9200)
>>>
>>>   Matrix products: default
>>>
>>>   locale:
>>>   [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United
>>>   Kingdom.1252    LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
>>>                           LC_TIME=English_United Kingdom.1252
>>>
>>>   attached base packages:
>>>   [1] stats     graphics  grDevices utils     datasets  methods   base
>>>
>>>   other attached packages:
>>>   [1] rgdal_1.4-4 sp_1.3-1
>>>
>>>   loaded via a namespace (and not attached):
>>>   [1] compiler_3.5.1  tools_3.5.1     yaml_2.2.0      grid_3.5.1
>>>   lattice_0.20-35
>>>
>>>   _______________________________________________
>>>   R-sig-Geo mailing list
>>>   [hidden email]
>>>   https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>>
>
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Re: rgdal::writeOGR with driver='ESRI Shapefile' converts Polygon object into a hole

Fri, 08/16/2019 - 14:27
On Fri, 16 Aug 2019, Roger Bivand wrote:

> On Tue, 13 Aug 2019, Phil Haines wrote:
>
>>  Dear list,
>>
>>  I have a single Polygons object containing multiple Polygon objects
>>  that share a common border. When I output this using writeOGR() one of
>>  the Polygon objects becomes a hole, as the following example shows.
>>
>>  Create a Polygons object containing two adjoining Polygon objects
>>>  library(rgdal)
>>>  r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
>>>  r2 <- r1; r2[,1] <- r2[,1]+1
>>>  Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
>>>  SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
>>>  data.frame(Example=c("Minimal")))
>>
>>  Perform a write/readOGR() cycle
>>>  fn <- tempfile()
>>>  writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
>>>  SPDF2 <- readOGR(dsn=fn,layer='test')
>>
>>  Second Polygon object is now a hole
>>>  sapply(SPDF2@polygons[[1]]@Polygons,slot,"hole")
>>  [1] FALSE  TRUE
>>
>>  I see from the sp documentation that "Polygon objects belonging to a
>>  Polygons object should either not overlap one-other, or should be
>>  fully included" but I am not sure how this relates to bordering
>>  Polygon objects. I would welcome any advice as to whether what I am
>>  asking of writeOGR is reasonable?
>
> The problem is with the 'ESRI Shapefile' representation and driver:
>
> library(rgdal)
> r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
> r2 <- r1; r2[,1] <- r2[,1]+1
> Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
> SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
> data.frame(Example=c("Minimal")))
> sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "hole")
> sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "ringDir")
>
> # which constructs a MULTIPOLYGON object:
>
> rgeos::writeWKT(SPDF)
> library(sf)
> st_as_text(st_geometry(st_as_sf(SPDF)))
>
> #  The 'ESRI Shapefile' driver is not Simple-Feature compliant (it pre-dates
> #  it), so the failure occurs by seeing the second exterior ring
> #  as an interior ring
>
> fn <- tempfile()
> writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
> SPDF2 <- readOGR(dsn=fn, layer='test')
> sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "hole")
> sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "ringDir")
> rgeos::writeWKT(SPDF2)
> st_as_text(st_geometry(st_as_sf(SPDF2)))
>
> # This happens with sf too, using the same GDAL driver:
>
> sf2 <- st_read(dsn=fn, layer='test')
> st_geometry(sf2)
> library(sf)
> st_as_text(st_geometry(st_as_sf(sf2)))
> rgeos::writeWKT(as(sf2, "Spatial"))
>
> # Adding the comment fix doesn't help:
>
> comment(slot(SPDF, "polygons")[[1]])
> SPDF_c <- rgeos::createSPComment(SPDF)
> comment(slot(SPDF_c, "polygons")[[1]])
> writeOGR(SPDF_c, fn, layer='test_c', driver='ESRI Shapefile',
>   verbose=TRUE)
>
> #  reports
> #  Object initially classed as: wkbPolygon
> #  SFS comments in Polygons objects
> #  Object reclassed as: wkbMultiPolygon
>
> SPDF2_c <- readOGR(dsn=fn, layer='test_c')
> sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot, "hole")
> sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot, "ringDir")
> rgeos::writeWKT(SPDF2_c)
> st_as_text(st_geometry(st_as_sf(SPDF2_c)))
>
>
>
> # If the input object is written out with the GeoPackage driver:
>
> fn1 <- tempfile(fileext=".gpkg")
> writeOGR(SPDF, fn1, layer="test", driver='GPKG')
> sf2a <- st_read(dsn=fn1,layer='test')
> st_coordinates(st_geometry(sf2a))
> SPDF2a <- readOGR(dsn=fn1)
> sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "hole")
> sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "ringDir")
> rgeos::writeWKT(SPDF2a)
> st_as_text(st_geometry(st_as_sf(SPDF2a)))
>
> # the issue is resolved. If we separate the exterior rings:
>
> r2a <- r1; r2a[,1] <- r2a[,1]+1.00001
> Ps1a = Polygons(list(Polygon(r1),Polygon(r2a)),ID=1)
> SPDFa = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1a)),
>  data.frame(Example=c("Minimal")))
> fna <- tempfile()
> writeOGR(SPDFa, fna, layer='test', driver='ESRI Shapefile')
> SPDF2_a <- readOGR(dsn=fna, layer='test')
> sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot, "hole")
> sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot, "ringDir")
> rgeos::writeWKT(SPDF2_a)
> st_as_text(st_geometry(st_as_sf(SPDF2_a)))
>
> # we are OK as the two exterior rings do not touch.
>
> # does using sf make a difference?
>
> fn_s <- tempfile(fileext=".shp")
> st_write(st_as_sf(SPDF), dsn=fn_s)
> sf_in <- st_read(fn_s)
> st_as_text(st_geometry(st_as_sf(sf_in)))
>
> # No
>
> fn_s <- tempfile(fileext=".shp")
> st_write(st_as_sf(SPDF_c), dsn=fn_s)
> sf_in_c <- st_read(fn_s)
> st_as_text(st_geometry(st_as_sf(sf_in_c)))
>
> # nor with the pretend-SF-compliant comment set either.
>
> So the weakness is in the "ESRI Shapefile" write driver, or possibly in the
> OGRGeometryFactory::organizePolygons() function in GDAL used in OGR_write()
> (a C++ function) called by writeOGR(). If sf::st_write() also calls
> OGRGeometryFactory::organizePolygons(), we'd maybe consider that it has a
> weakness for the "ESRI Shapefile" driver, but which does not affect
> SF-compliant drivers.
Without the comment set, OGRGeometryFactory::organizePolygons() is used;
with it set, OGRGeometryFactory::organizePolygons() is not used, because
the object is declared as two exterior rings. In both cases, we have the
output object written out and read back in incorrectly with the ESRI
shapefile driver, but SF-compliant drivers round-trip (in the test
GeoJSON), correctly.

It is likely that the changes made in 2015 to accommodate GeoJSON led to
this possible regression for the ESRI Shapefile driver. I'm adding this
geometry to tests/tripup.R and data files; without code changes the hole
slot is wrong and the ring direction is changes to match, so the
coordinates change too.

Reading using the deprecated maptools::readShapeSpatial() also gets a hole
in the second external ring. However, writing with deprecated
maptools::writeSpatialShape() yields a shapefile that when read with
maptools::readShapeSpatial() gives the correct two exterior ring, no hole
object. When read with sf::st_read() and rgdal::readOGR(), the object is
also correct. So the problem definitely lies in rgdal::writeOGR(), and
sf::st_write() - roundtripping with sf::st_write() and sf::st_read()
degrades from MULTIPOLYGON to POLYGON with the ESRI Shapefile driver.

Roger

>
> This is probably related to a similar but inverse problem with the
> SF-compliant GeoJSON driver in 2015:
>
> https://stat.ethz.ch/pipermail/r-sig-geo/2015-October/023609.html
>
> continued the next month in:
>
> https://stat.ethz.ch/pipermail/r-sig-geo/2015-November/023656.html
>
> The details are in this SVN diff
>
> https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=555&r2=571
>
> up to the end og the list thread, and this one from then until now:
>
> https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=571&r2=733
>
> Summary: could you change drivers, or is it really necessary to fix an EOL
> problem? What is your use case?
>
>>
>>  Thanks in advance for your time,
>
> Thanks for a complete example,
>
> Roger
>
>>  Phil
>>
>>>  sessionInfo()
>>  R version 3.5.1 (2018-07-02)
>>  Platform: x86_64-w64-mingw32/x64 (64-bit)
>>  Running under: Windows >= 8 x64 (build 9200)
>>
>>  Matrix products: default
>>
>>  locale:
>>  [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United
>>  Kingdom.1252    LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
>>                          LC_TIME=English_United Kingdom.1252
>>
>>  attached base packages:
>>  [1] stats     graphics  grDevices utils     datasets  methods   base
>>
>>  other attached packages:
>>  [1] rgdal_1.4-4 sp_1.3-1
>>
>>  loaded via a namespace (and not attached):
>>  [1] compiler_3.5.1  tools_3.5.1     yaml_2.2.0      grid_3.5.1
>>  lattice_0.20-35
>>
>>  _______________________________________________
>>  R-sig-Geo mailing list
>>  [hidden email]
>>  https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Re: rgdal::writeOGR with driver='ESRI Shapefile' converts Polygon object into a hole

Fri, 08/16/2019 - 07:22
On Tue, 13 Aug 2019, Phil Haines wrote:

> Dear list,
>
> I have a single Polygons object containing multiple Polygon objects
> that share a common border. When I output this using writeOGR() one of
> the Polygon objects becomes a hole, as the following example shows.
>
> Create a Polygons object containing two adjoining Polygon objects
>> library(rgdal)
>> r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
>> r2 <- r1; r2[,1] <- r2[,1]+1
>> Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
>> SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)), data.frame(Example=c("Minimal")))
>
> Perform a write/readOGR() cycle
>> fn <- tempfile()
>> writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
>> SPDF2 <- readOGR(dsn=fn,layer='test')
>
> Second Polygon object is now a hole
>> sapply(SPDF2@polygons[[1]]@Polygons,slot,"hole")
> [1] FALSE  TRUE
>
> I see from the sp documentation that "Polygon objects belonging to a
> Polygons object should either not overlap one-other, or should be
> fully included" but I am not sure how this relates to bordering
> Polygon objects. I would welcome any advice as to whether what I am
> asking of writeOGR is reasonable?
The problem is with the 'ESRI Shapefile' representation and driver:

library(rgdal)
r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
r2 <- r1; r2[,1] <- r2[,1]+1
Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)),
  data.frame(Example=c("Minimal")))
sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "hole")
sapply(slot(slot(SPDF, "polygons")[[1]], "Polygons"), slot, "ringDir")

# which constructs a MULTIPOLYGON object:

rgeos::writeWKT(SPDF)
library(sf)
st_as_text(st_geometry(st_as_sf(SPDF)))

# The 'ESRI Shapefile' driver is not Simple-Feature compliant (it
# pre-dates it), so the failure occurs by seeing the second exterior ring
# as an interior ring

fn <- tempfile()
writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
SPDF2 <- readOGR(dsn=fn, layer='test')
sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "hole")
sapply(slot(slot(SPDF2, "polygons")[[1]], "Polygons"), slot, "ringDir")
rgeos::writeWKT(SPDF2)
st_as_text(st_geometry(st_as_sf(SPDF2)))

# This happens with sf too, using the same GDAL driver:

sf2 <- st_read(dsn=fn, layer='test')
st_geometry(sf2)
library(sf)
st_as_text(st_geometry(st_as_sf(sf2)))
rgeos::writeWKT(as(sf2, "Spatial"))

# Adding the comment fix doesn't help:

comment(slot(SPDF, "polygons")[[1]])
SPDF_c <- rgeos::createSPComment(SPDF)
comment(slot(SPDF_c, "polygons")[[1]])
writeOGR(SPDF_c, fn, layer='test_c', driver='ESRI Shapefile',
   verbose=TRUE)

# reports
# Object initially classed as: wkbPolygon
# SFS comments in Polygons objects
# Object reclassed as: wkbMultiPolygon

SPDF2_c <- readOGR(dsn=fn, layer='test_c')
sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot, "hole")
sapply(slot(slot(SPDF2_c, "polygons")[[1]], "Polygons"), slot, "ringDir")
rgeos::writeWKT(SPDF2_c)
st_as_text(st_geometry(st_as_sf(SPDF2_c)))



# If the input object is written out with the GeoPackage driver:

fn1 <- tempfile(fileext=".gpkg")
writeOGR(SPDF, fn1, layer="test", driver='GPKG')
sf2a <- st_read(dsn=fn1,layer='test')
st_coordinates(st_geometry(sf2a))
SPDF2a <- readOGR(dsn=fn1)
sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "hole")
sapply(slot(slot(SPDF2a, "polygons")[[1]], "Polygons"), slot, "ringDir")
rgeos::writeWKT(SPDF2a)
st_as_text(st_geometry(st_as_sf(SPDF2a)))

# the issue is resolved. If we separate the exterior rings:

r2a <- r1; r2a[,1] <- r2a[,1]+1.00001
Ps1a = Polygons(list(Polygon(r1),Polygon(r2a)),ID=1)
SPDFa = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1a)),
   data.frame(Example=c("Minimal")))
fna <- tempfile()
writeOGR(SPDFa, fna, layer='test', driver='ESRI Shapefile')
SPDF2_a <- readOGR(dsn=fna, layer='test')
sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot, "hole")
sapply(slot(slot(SPDF2_a, "polygons")[[1]], "Polygons"), slot, "ringDir")
rgeos::writeWKT(SPDF2_a)
st_as_text(st_geometry(st_as_sf(SPDF2_a)))

# we are OK as the two exterior rings do not touch.

# does using sf make a difference?

fn_s <- tempfile(fileext=".shp")
st_write(st_as_sf(SPDF), dsn=fn_s)
sf_in <- st_read(fn_s)
st_as_text(st_geometry(st_as_sf(sf_in)))

# No

fn_s <- tempfile(fileext=".shp")
st_write(st_as_sf(SPDF_c), dsn=fn_s)
sf_in_c <- st_read(fn_s)
st_as_text(st_geometry(st_as_sf(sf_in_c)))

# nor with the pretend-SF-compliant comment set either.

So the weakness is in the "ESRI Shapefile" write driver, or possibly in
the OGRGeometryFactory::organizePolygons() function in GDAL used in
OGR_write() (a C++ function) called by writeOGR(). If sf::st_write() also
calls OGRGeometryFactory::organizePolygons(), we'd maybe consider that it
has a weakness for the "ESRI Shapefile" driver, but which does not affect
SF-compliant drivers.

This is probably related to a similar but inverse problem with the
SF-compliant GeoJSON driver in 2015:

https://stat.ethz.ch/pipermail/r-sig-geo/2015-October/023609.html

continued the next month in:

https://stat.ethz.ch/pipermail/r-sig-geo/2015-November/023656.html

The details are in this SVN diff

https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=555&r2=571

up to the end og the list thread, and this one from then until now:

https://r-forge.r-project.org/scm/viewvc.php/pkg/src/OGR_write.cpp?root=rgdal&r1=571&r2=733

Summary: could you change drivers, or is it really necessary to fix an EOL
problem? What is your use case?

>
> Thanks in advance for your time,

Thanks for a complete example,

Roger

> Phil
>
>> sessionInfo()
> R version 3.5.1 (2018-07-02)
> Platform: x86_64-w64-mingw32/x64 (64-bit)
> Running under: Windows >= 8 x64 (build 9200)
>
> Matrix products: default
>
> locale:
> [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United
> Kingdom.1252    LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
>                         LC_TIME=English_United Kingdom.1252
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] rgdal_1.4-4 sp_1.3-1
>
> loaded via a namespace (and not attached):
> [1] compiler_3.5.1  tools_3.5.1     yaml_2.2.0      grid_3.5.1
> lattice_0.20-35
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Forcing use of alpha-band in RGBA, as alpha for other bands in plotRGB

Thu, 08/15/2019 - 09:43
Can one ensure that the alpha-band (i.e. band 4) of an RGBA RasterBrick is applied as a mask when using plotRGB(my_brick)? 

I have a GeoTIFF with an alpha-band marking an irregularly-shaped mask, which displays correctly in QGIS (using band 4 as transparency, values either simply 0 for mask or 255 for visible), but I cannot get plotRGB to display the masked area as anything other than black (which simply means that RGB bands 1,2,3 all =0, and that the alpha mask is being ignored). 

summary(my_brick) acknowledges the presence of the 4th band, and plot(my_brick[[4]]) shows the correctly-shaped mask. 

I hope there is a blindingly obvious solution to this, but have singularly failed to find anything to help me via normal web searches or from the manual. 

(Reclassification is not ideal given that this particular raster is very large and eats memory). 

Thank you in advance, 
Toby


        [[alternative HTML version deleted]]

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

Pages