date:20170315

Re: [R-sig-Geo] Maintain SRID with st_write to postgis db

2017-03-15 Thread chris english

Michael,

I have the very same failure,
> library(sf)
Linking to GEOS 3.4.2, GDAL 2.1.3, proj.4 4.9.1
(pre upgrade attempt, which still works just fine)

checking for geos_c.h... yes
checking geos: linking with -lgeos_c... yes
checking for lwgeom_version in -llwgeom... no
configure: error: in `/tmp/RtmpsCmf6B/devtoolsd68e68107a/edzer-sfr-d620f3b':
configure: error: lwgeom test failed (--without-readline to disable)

I could nearly copy your install output though on 16.04, and the same
tangle of it depends on x but it won't be installed & etc.
But, if you have PostGIS, you have liblwgeom as it is fairly specific to
PostGIS. Perhaps, rather than a devtools::github_install,
a clone github lo a local directory and config, make, install might work,
but I'm just imagining.

then testing:
config.log --snip --
configure:3993: checking for lwgeom_version in -llwgeom
configure:4018: gcc -std=gnu99 -o conftest -g -O2 -fstack-protector
--param=ssp-
buffer-size=4 -Wformat -Werror=format-security -D_FORTIFY_SOURCE=2 -g
 -I/usr/lo
cal/include   -I/usr/local/include -I/usr/local/include
-Wl,-Bsymbolic-functions
 -Wl,-z,relro conftest.c -llwgeom  -lproj  -L/usr/local/lib -lproj
-L/usr/loca
l/lib -lgdal -L/usr/local/lib -lgeos_c -lproj >&5
/tmp/ccDBUh1h.o: In function `main':
/home/chris/sfr/sfr/conftest.c:34: undefined reference to `lwgeom_version'
collect2: error: ld returned 1 exit status
configure:4018: $? = 1
configure: failed program was:
| /* confdefs.h */
| #define PACKAGE_NAME ""
| #define PACKAGE_TARNAME ""
| #define PACKAGE_VERSION ""
| #define PACKAGE_STRING ""
| #define PACKAGE_BUGREPORT ""
| #define PACKAGE_URL ""
| #define STDC_HEADERS 1
| #define HAVE_SYS_TYPES_H 1
| #define HAVE_SYS_STAT_H 1
| #define HAVE_STDLIB_H 1
| #define HAVE_STRING_H 1
| #define HAVE_MEMORY_H 1
| #define HAVE_STRINGS_H 1
| #define HAVE_INTTYPES_H 1
| #define HAVE_STDINT_H 1
| #define HAVE_UNISTD_H 1
| #define HAVE_GDAL_H 1
| #define HAVE_PROJ_API_H 1
| #define HAVE_LIBPROJ 1
| #define HAVE_GEOS_C_H 1
| /* end confdefs.h.  */
|
| /* Override any GCC internal prototype to avoid an error.
|Use char because int might match the return type of a GCC
|builtin and then its argument prototype would still apply.  */
| #ifdef __cplusplus
| extern "C"
| #endif
| char lwgeom_version ();
| int
| main ()
| {
| return lwgeom_version ();
|   ;
|   return 0;
| }
configure:4027: result: no
configure:4036: error: in `/home/chris/sfr/sfr':
configure:4038: error: lwgeom test failed (--without-readline to disable)
See `config.log' for more details

so, it may or may not be the presence or absence of liblwgeom but simply
an undefined reference as suggested above:

/home/chris/sfr/sfr/conftest.c:34: undefined reference to `lwgeom_version' ,

but I'm unclear how the version was being requested (well, I can kind of
guess)
but failing on undefined reference suggests it is not the version per se
that may
or may not be present (though is because you are using PostGIS), but how
the version was requested. I intuit.

Ah, and reading  /* confdefs.h */, that it indeed ends with #define
HAVE_GEOS_C_H 1,
and indeed, conftest.c:34: would have an  undefined reference to
`lwgeom_version'.

And I've said as much as I understand.

Chris

On Wed, Mar 15, 2017 at 7:53 PM, Michael Treglia  wrote:

> Sorry - this follow-up isn't entirely an R question, so if best to take
> this off list, let me know:
>
> Installing the dev version from github fails for me (installing on ubuntu
> 14.04.5) - I've included the output from the install process below - seems
> to fail due to failed check for liblwgeom version.
>
> Looks like I don't have liblwgeom installed - and when I try (via sudo
> apt-get install liblwgeom-2.3-0), it indicates it requires libgeos-c1.
> Installing libgeos-c1, however, requires removal of a newer version of
> libgeos (libgeos-c1v5), which other FOSS GIS tools depend on (at least if
> my understanding is correct).  Is there a way around this?  Sorry if I'm
> just missing something, and thanks again for input.
> mike
>
>
>
> #Output of install from github
> > install_github("edzer/sfr")
> Downloading GitHub repo edzer/sfr@master
> from URL https://api.github.com/repos/edzer/sfr/zipball/master
> Installing sf
> '/usr/lib/R/bin/R' --no-site-file --no-environ --no-save --no-restore
> --quiet  \
>   CMD INSTALL '/tmp/RtmpKLmhyt/devtools9fd11eb9c3/edzer-sfr-d620f3b'  \
>   --library='/home/mlt244/R/x86_64-pc-linux-gnu-library/3.3'
> --install-tests
>
> * installing *source* package ‘sf’ ...
> configure: CC: gcc -std=gnu99
> configure: CXX: g++
> configure: :
> checking for gdal-config... /usr/bin/gdal-config
> checking gdal-config usability... yes
> configure: GDAL: 2.1.0
> checking GDAL version >= 2.0.0... yes
> checking for gcc... gcc -std=gnu99
> checking whether the C compiler works... yes
> checking for C compiler default output file name... a.out
> checking for suffix of executables...
> checking whether we are cross compiling... no
>

Re: [R-sig-Geo] Maintain SRID with st_write to postgis db

2017-03-15 Thread Michael Treglia

Sorry - this follow-up isn't entirely an R question, so if best to take
this off list, let me know:

Installing the dev version from github fails for me (installing on ubuntu
14.04.5) - I've included the output from the install process below - seems
to fail due to failed check for liblwgeom version.

Looks like I don't have liblwgeom installed - and when I try (via sudo
apt-get install liblwgeom-2.3-0), it indicates it requires libgeos-c1.
Installing libgeos-c1, however, requires removal of a newer version of
libgeos (libgeos-c1v5), which other FOSS GIS tools depend on (at least if
my understanding is correct).  Is there a way around this?  Sorry if I'm
just missing something, and thanks again for input.
mike

#Output of install from github
> install_github("edzer/sfr")
Downloading GitHub repo edzer/sfr@master
from URL https://api.github.com/repos/edzer/sfr/zipball/master
Installing sf
'/usr/lib/R/bin/R' --no-site-file --no-environ --no-save --no-restore
--quiet  \
  CMD INSTALL '/tmp/RtmpKLmhyt/devtools9fd11eb9c3/edzer-sfr-d620f3b'  \
  --library='/home/mlt244/R/x86_64-pc-linux-gnu-library/3.3'
--install-tests

* installing *source* package ‘sf’ ...
configure: CC: gcc -std=gnu99
configure: CXX: g++
configure: :
checking for gdal-config... /usr/bin/gdal-config
checking gdal-config usability... yes
configure: GDAL: 2.1.0
checking GDAL version >= 2.0.0... yes
checking for gcc... gcc -std=gnu99
checking whether the C compiler works... yes
checking for C compiler default output file name... a.out
checking for suffix of executables...
checking whether we are cross compiling... no
checking for suffix of object files... o
checking whether we are using the GNU C compiler... yes
checking whether gcc -std=gnu99 accepts -g... yes
checking for gcc -std=gnu99 option to accept ISO C89... none needed
checking how to run the C preprocessor... gcc -std=gnu99 -E
checking for grep that handles long lines and -e... /bin/grep
checking for egrep... /bin/grep -E
checking for ANSI C header files... yes
checking for sys/types.h... yes
checking for sys/stat.h... yes
checking for stdlib.h... yes
checking for string.h... yes
checking for memory.h... yes
checking for strings.h... yes
checking for inttypes.h... yes
checking for stdint.h... yes
checking for unistd.h... yes
checking gdal.h usability... yes
checking gdal.h presence... yes
checking for gdal.h... yes
checking GDAL: linking with --libs only... yes
checking GDAL: /usr/share/gdal/2.1/pcs.csv readable... yes
checking GDAL: checking whether PROJ.4 is available for linking:... yes
checking GDAL: checking whether PROJ.4 is available fur running:... yes
checking proj_api.h usability... yes
checking proj_api.h presence... yes
checking for proj_api.h... yes
checking for pj_init_plus in -lproj... yes
checking PROJ.4: epsg found and readable... yes
proj_conf_test.c: In function 'main':
proj_conf_test.c:5:5: error: unknown type name 'PAFile'
 PAFile fp;
 ^
proj_conf_test.c:8:5: warning: implicit declaration of function
'pj_open_lib' [-Wimplicit-function-declaration]
 fp = pj_open_lib(ctx, "conus", "rb");
 ^
proj_conf_test.c:9:12: warning: comparison between pointer and integer
[enabled by default]
 if (fp == NULL) exit(1);
^
proj_conf_test.c:10:5: warning: implicit declaration of function
'pj_ctx_fclose' [-Wimplicit-function-declaration]
 pj_ctx_fclose(ctx, fp);
 ^
./configure: line 3834: ./proj_conf_test: No such file or directory
checking PROJ.4: conus found and readable... yes
checking geos-config usability... yes
configure: GEOS: 3.5.0
checking GEOS version >= 3.3.0... yes
checking geos_c.h usability... yes
checking geos_c.h presence... yes
checking for geos_c.h... yes
checking geos: linking with -lgeos_c... yes
checking for lwgeom_version in -llwgeom... no
configure: error: in `/tmp/RtmpKLmhyt/devtools9fd11eb9c3/edzer-sfr-d620f3b':
configure: error: lwgeom test failed (--without-readline to disable)
See `config.log' for more details
ERROR: configuration failed for package ‘sf’
* removing ‘/home/mlt244/R/x86_64-pc-linux-gnu-library/3.3/sf’
* restoring previous ‘/home/mlt244/R/x86_64-pc-linux-gnu-library/3.3/sf’
Error: Command failed (1)

On Wed, Mar 15, 2017 at 6:14 PM, Michael Treglia  wrote:

> Wow - thanks so much for this quick fix, Edzer! I like your solution to
> having syntatically different but semantically identical proj4string. Will
> try this in a bit, and write back if I have any issues.
> Best,
> mike
>
> On Wed, Mar 15, 2017 at 5:58 PM, Edzer Pebesma <
> edzer.pebe...@uni-muenster.de> wrote:
>
>> Great question! Short answer: now solved in the dev version on github;
>> data are written directly to postgis having epsg 2263.
>>
>>
>> Long answer:
>>
>> I believe this was caused by gdal and proj.4 doing different things when
>> resolving epsg codes to proj4strings. sf uses proj.4 for this. When
>> writing a proj4string either gdal or postgis don't recognize the
>> proj4string that proj.4 returns on the epsg

Re: [R-sig-Geo] Maintain SRID with st_write to postgis db

2017-03-15 Thread Michael Treglia

Wow - thanks so much for this quick fix, Edzer! I like your solution to
having syntatically different but semantically identical proj4string. Will
try this in a bit, and write back if I have any issues.
Best,
mike

On Wed, Mar 15, 2017 at 5:58 PM, Edzer Pebesma <
edzer.pebe...@uni-muenster.de> wrote:

> Great question! Short answer: now solved in the dev version on github;
> data are written directly to postgis having epsg 2263.
>
>
> Long answer:
>
> I believe this was caused by gdal and proj.4 doing different things when
> resolving epsg codes to proj4strings. sf uses proj.4 for this. When
> writing a proj4string either gdal or postgis don't recognize the
> proj4string that proj.4 returns on the epsg code. Neither gdal nor
> postgis understand that syntactically different proj4strings may be
> semantically identical.
>
>
> After running your example, in postgis
>
> # select proj4text from spatial_ref_sys where SRID = 900914;
>
> gives:
>
>  +proj=lcc +lat_1=41.03 +lat_2=40.66
> +lat_0=40.16 +lon_0=-74 +x_0=30 +y_0=0 +ellps=GRS80
> +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs
>
> # select proj4text from spatial_ref_sys where SRID = 2263;
>
> gives:
>
>  +proj=lcc +lat_1=41.03 +lat_2=40.66
> +lat_0=40.16 +lon_0=-74 +x_0=30.01 +y_0=0
> +datum=NAD83 +units=us-ft +no_defs
>
> sf causes this:
>
> > st_crs(2263)
> $epsg
> [1] 2263
>
> $proj4string
> [1] "+proj=lcc +lat_1=41.03 +lat_2=40.66
> +lat_0=40.16 +lon_0=-74 +x_0=30.01 +y_0=0
> +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs"
>
> attr(,"class")
> [1] "crs"
>
> because it uses proj.4 directly to resolve SRID strings:
>
> /usr/share/proj/epsg has:
>
> # NAD83 / New York Long Island (ftUS)
> <2263> +proj=lcc +lat_1=41.03 +lat_2=40.66
> +lat_0=40.16 +lon_0=-74 +x_0=30.01 +y_0=0
> +datum=NAD83 +units=us-ft +no_defs  <>
>
>
> The change I made to sf (
> https://github.com/edzer/sfr/commit/d620f3b748d487bae910e3884c554f
> a9e3ce7090)
> now first asks gdal to resolve the epsg (in case it is not NA), and
> otherwise resolve the proj4string (if not NA), instead of ONLY trying to
> resolve a non-NA proj4string.
>
> On 15/03/17 20:50, Michael Treglia wrote:
> > Hi All,
> >
> > Been working to import and manipulate a CSV file with point data
> > (lat/long), and then export to a PostGIS db.
> >
> > Overall, successful, but one thing I'd like to fix - when I write out the
> > layer to postgis, the SRID is not maintained. The final EPSG/SRID should
> be
> > 2263, but when I check in PostGIS, it comes up as 900914.
> >
> > Below is my code and sessionInfo, and the data are from here:
> > https://data.cityofnewyork.us/Public-Safety/NYPD-Complaint-
> Data-Current-YTD/5uac-w243
> > (downloaded as plain old CSV)
> >
> > Anything I might be missing? Thanks in advance for giving a quick look!
> > Mike
> >
> >
> > ##Start Code
> >
> > #load packages
> > library(sf)
> > library(RPostgreSQL)
> >
> > #read data
> > crime_current <- read.csv("NYPD_Complaint_Data_Current_YTD.csv",
> > stringsAsFactors = FALSE)
> >
> > #format time columns for easier reading in postgres (I think), as keeping
> > as date format caused problems in export
> > crime_current$CMPLNT_FR_TIME <-
> > as.character(as.POSIXct(paste(crime_current$CMPLNT_FR_DT,
> > crime_current$CMPLNT_FR_TM), format="%m/%d/%Y\ %H:%M", tz=""))
> > crime_current$CMPLNT_TO_TIME <-
> > as.character(as.POSIXct(paste(crime_current$CMPLNT_TO_DT,
> > crime_current$CMPLNT_TO_TM), format="%m/%d/%Y\ %H:%M", tz=""))
> > crime_current$RPT_DT <- as.character(as.POSIXct(crime_current$RPT_DT,
> > format="%m/%d/%Y", tz=""))
> >
> > #convert to sf object
> > crime_current.sf <- st_as_sf(crime_current, coords = c("Longitude",
> > "Latitude"), crs = 4326)
> > #reproject to EPSG 2263
> > crime_current.sf <- st_transform(crime_current.sf, crs=2263)
> >
> > #write to postgres
> > st_write(crime_current.sf, "PG:dbname=mydb user=user host=xx.xx.xx.xx",
> > 'health_safety.crime_current')
> > ###End Code
> >
> >
> >
> >> sessionInfo()
> > R version 3.3.1 (2016-06-21)
> > Platform: x86_64-pc-linux-gnu (64-bit)
> > Running under: Ubuntu 14.04.5 LTS
> >
> > locale:
> >  [1] LC_CTYPE=en_US.UTF-8   LC_NUMERIC=C
> > LC_TIME=en_US.UTF-8LC_COLLATE=en_US.UTF-8
> >  [5] LC_MONETARY=en_US.UTF-8LC_MESSAGES=en_US.UTF-8
> >  LC_PAPER=en_US.UTF-8   LC_NAME=C
> >  [9] LC_ADDRESS=C   LC_TELEPHONE=C
> > LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
> >
> > attached base packages:
> > [1] stats graphics  grDevices utils datasets  methods   base
> >
> > other attached packages:
> > [1] sp_1.2-3  RPostgreSQL_0.4-1 DBI_0.6   sf_0.3-4
> >
> > loaded via a namespace (and not attached):
> > [1] tools_3.3.1 units_0.4-2 Rcpp_0.12.9 udunits2_0.13
> > grid_3.3.1  lattice_0.20-33
> >
> >   [[alternative HTML version

Re: [R-sig-Geo] Sentienl 2 gdal translate

2017-03-15 Thread Loïc Dutrieux

Hi,

I tried your code with some S2 images I had lying around, and it works
on my system. I modified it a bit to write the output layers to the same
directory and not to my working directory.

library(gdalUtils)
library(raster)

Image_Path<-'/path/to/images/'

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

for (file in S2_JP2_List) {
  out_file <- extension(file, 'tif')
  gdal_translate(src_dataset = file, dst_dataset = out_file, ot =
"UInt16", of = "GTiff")
}

gdal has a sentinel 2 driver, but it's not extremely old, therefore it's
possible that your installation doesn't have it; you can check with:

gdal-config --formats | grep sentinel2

Cheers,
Loïc



On 15/03/17 12:49, Miguel Castro Gómez wrote:
> Hi,
> 
> I have Sentinel 2 images (jp2) that I want to convert to GTiff using 
> gdal_translate in R (from gdalutils package). Here is the code I’m using
> 
> Image_Path<-“/path/to/wd“
> 
> S2_JP2_List <- list.files(Image_Path, full.names = TRUE, pattern = ".jp2$")
> 
> for (i in 1:length(S2_JP2_List)) {
> gdal_translate(src_dataset = S2_JP2_List[[i]], dst_dataset = paste("Band", i 
> , "tif"), ot = "UInt16", of = “GTiff")
> }
> 
> When running this code, the process starts without any error but its never 
> ending neither producing any output
> 
> I've tried to do it in gdal and the same code works, except that it is 
> necessary to skip  the default driver since it cannot manage large jp2 files. 
> 
> gdal_translate -of GTiff /path/to/input/S2_b1.jp2  
> /path/to/output/S2_b1converted.tif --config GDAL_SKIP JP2ECW
> 
> Any idea how could I run this in R?
> 
> Thanks 
> M
> 
> 
> 
>   [[alternative HTML version deleted]]
> 
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

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

Re: [R-sig-Geo] Maintain SRID with st_write to postgis db

2017-03-15 Thread Michael Treglia

Hi Chris,

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

Best,
Mike

On Wed, Mar 15, 2017 at 5:21 PM, chris english <
englishchristoph...@gmail.com> wrote:

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

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] Maintain SRID with st_write to postgis db

2017-03-15 Thread chris english

Michael,

this from the sf vignettes on read write to databases:

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

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

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

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

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

HTH,

Chris

On Wed, Mar 15, 2017 at 3:50 PM, Michael Treglia  wrote:

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

[[alternative HTML version deleted]]

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

[R-sig-Geo] Maintain SRID with st_write to postgis db

2017-03-15 Thread Michael Treglia

Hi All,

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

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

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

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


##Start Code

#load packages
library(sf)
library(RPostgreSQL)

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

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

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

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



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

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

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

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

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

[[alternative HTML version deleted]]

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

[R-sig-Geo] Sentienl 2 gdal translate

2017-03-15 Thread Miguel Castro Gómez

Hi,

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

Image_Path<-“/path/to/wd“

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

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

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

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

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

Any idea how could I run this in R?

Thanks 
M



[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] loop memory problem

2017-03-15 Thread marta azores

Hi again,

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

Any idea?

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

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

stopCluster(cl)
registerDoSEQ()

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

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] Same mask for two rasters, but obtaining different extents in R

2017-03-15 Thread Sarah Goslee

The most likely explanation seems to me that the rasters were not
aligned before cropping, so they are not aligned after cropping.

Have you checked the original extent and resolution?

Sarah

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

-- 
Sarah Goslee
http://www.functionaldiversity.org

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

Re: [R-sig-Geo] Same mask for two rasters, but obtaining different extents in R

2017-03-15 Thread Ben Tupper

Hi,

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

Cheers,
Ben


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

Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org

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

[R-sig-Geo] Same mask for two rasters, but obtaining different extents in R

2017-03-15 Thread Vijay Ramesh via R-sig-Geo

I loaded a shapefile in R, and cropped and masked it with the same mask,
but I get different extents. Any suggestions?

Code below:

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

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

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

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

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

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

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


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

[[alternative HTML version deleted]]

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

Re: [R-sig-Geo] Maintain SRID with st_write to postgis db

Re: [R-sig-Geo] Maintain SRID with st_write to postgis db

Re: [R-sig-Geo] Maintain SRID with st_write to postgis db

Re: [R-sig-Geo] Sentienl 2 gdal translate

Re: [R-sig-Geo] Maintain SRID with st_write to postgis db

Re: [R-sig-Geo] Maintain SRID with st_write to postgis db

[R-sig-Geo] Maintain SRID with st_write to postgis db

[R-sig-Geo] Sentienl 2 gdal translate

Re: [R-sig-Geo] loop memory problem

Re: [R-sig-Geo] Same mask for two rasters, but obtaining different extents in R

Re: [R-sig-Geo] Same mask for two rasters, but obtaining different extents in R

[R-sig-Geo] Same mask for two rasters, but obtaining different extents in R

12 matches

Site Navigation

Mail list logo

Footer information