[R-sig-Geo] converting raster to netcdf .nc with keeping CRS information

[Ahmed El-Gabbas]

Hello all,


I am having a problem while converting a raster object as NetCDF (.nc)
file, with keeping the CRS information in the output file.

Here is a reproducible code:

require(raster)
require(ncdf4)
CurrTemp <- tempfile()
download.file(url =
"https://seaice.uni-bremen.de/data/amsre/asi_daygrid_swath/s6250/2003/feb/Antarctic/asi-s6250-20030214-v5.hdf";,
destfile = CurrTemp, mode = "wb", quiet = T)
r <- raster(CurrTemp)
r <- flip(r,2)
extent(r) <- c(-395, 395, -395, 435)
crs(r) <- "+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0
+y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs"
r# class   : RasterLayer # dimensions  : 1328, 1264, 1678592
(nrow, ncol, ncell)# resolution  : 6250, 6250  (x, y)# extent  :
-395, 395, -395, 435  (xmin, xmax, ymin, ymax)# coord.
ref. : +proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0
+a=6378273 +b=6356889.449 +units=m +no_defs # data source : in memory#
names   : layer # values  : 0, 100  (min, max)

So far the raster object reads well, with CRS information included.

However, when I try to save it as .nc file, R prints "coord. ref. : NA",
and the produced file does not contain the CRS information.

writeRaster(r, filename = "O:/Ahmed/r001.nc", varname="IceConc",
overwrite=TRUE, format="CDF",
xname="Longitude", yname="Latitude")# class   :
RasterLayer # dimensions  : 1328, 1264, 1678592  (nrow, ncol, ncell)#
resolution  : 6250, 6250  (x, y)# extent  : -395, 395,
-395, 435  (xmin, xmax, ymin, ymax)# coord. ref. : NA # data
source : O:/Ahmed/r001.nc # names   : IceConc # zvar:
IceConc

raster("O:/Ahmed/r001.nc")# class   : RasterLayer # dimensions  :
1328, 1264, 1678592  (nrow, ncol, ncell)# resolution  : 6250, 6250
(x, y)# extent  : -395, 395, -395, 435  (xmin,
xmax, ymin, ymax)# coord. ref. : NA # data source : O:/Ahmed/r001.nc #
names   : IceConc # zvar: IceConc

Any solution?


N.B. I also sent the same question at stackoverflow
,
without a solution so far.


Thanks

~
*Dr. Ahmed El-Gabbas,*
Ocean Acoustics Lab,
Alfred-Wegener-Institut
Helmholtz-Zentrum für Polar und Meeresforschung
[image: image.png]
My Website:  https://elgabbas.github.io
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] converting raster to netcdf .nc with keeping CRS information

[Ben Tupper]

Hi,

Have you tried using the 'prj' argument to writeRaster?  I don't know that it 
will be the solution, but according to ...

https://www.rdocumentation.org/packages/raster/versions/2.8-19/topics/writeRaster
 


It is documented among the '...' arguments.  Setting prj = TRUE will cause the 
crs to be written to the file.

Cheers,
Ben

> On Feb 12, 2019, at 8:34 AM, Ahmed El-Gabbas  wrote:
> 
> Hello all,
> 
> I am having a problem while converting a raster object as NetCDF (.nc) file, 
> with keeping the CRS information in the output file.
> Here is a reproducible code:
> require(raster)
> require(ncdf4)
> CurrTemp <- tempfile()
> download.file(url = 
> "https://seaice.uni-bremen.de/data/amsre/asi_daygrid_swath/s6250/2003/feb/Antarctic/asi-s6250-20030214-v5.hdf
>  
> ",
>  destfile = CurrTemp, mode = "wb", quiet = T)
> r <- raster(CurrTemp)
> r <- flip(r,2)
> extent(r) <- c(-395, 395, -395, 435)
> crs(r) <- "+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 
> +a=6378273 +b=6356889.449 +units=m +no_defs"
> r
> # class   : RasterLayer 
> # dimensions  : 1328, 1264, 1678592  (nrow, ncol, ncell)
> # resolution  : 6250, 6250  (x, y)
> # extent  : -395, 395, -395, 435  (xmin, xmax, ymin, ymax)
> # coord. ref. : +proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 
> +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs 
> # data source : in memory
> # names   : layer 
> # values  : 0, 100  (min, max)
> So far the raster object reads well, with CRS information included.
> However, when I try to save it as .nc file, R prints "coord. ref. : NA", and 
> the produced file does not contain the CRS information.
> writeRaster(r, filename = "O:/Ahmed/r001.nc ", 
> varname="IceConc",
> overwrite=TRUE, format="CDF",
> xname="Longitude", yname="Latitude")
> # class   : RasterLayer 
> # dimensions  : 1328, 1264, 1678592  (nrow, ncol, ncell)
> # resolution  : 6250, 6250  (x, y)
> # extent  : -395, 395, -395, 435  (xmin, xmax, ymin, ymax)
> # coord. ref. : NA 
> # data source : O:/Ahmed/r001.nc  
> # names   : IceConc 
> # zvar: IceConc
> 
> raster("O:/Ahmed/r001.nc ")
> # class   : RasterLayer 
> # dimensions  : 1328, 1264, 1678592  (nrow, ncol, ncell)
> # resolution  : 6250, 6250  (x, y)
> # extent  : -395, 395, -395, 435  (xmin, xmax, ymin, ymax)
> # coord. ref. : NA 
> # data source : O:/Ahmed/r001.nc  
> # names   : IceConc 
> # zvar: IceConc
> Any solution?
> 
> N.B. I also sent the same question at stackoverflow 
> ,
>  without a solution so far.
> 
> Thanks
> 
> ~
> Dr. Ahmed El-Gabbas,
> Ocean Acoustics Lab,
> Alfred-Wegener-Institut
> Helmholtz-Zentrum für Polar und Meeresforschung
> 
> My Website:  https://elgabbas.github.io 
> ___
> 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

Ecological Forecasting: https://eco.bigelow.org/






[[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] converting raster to netcdf .nc with keeping CRS information

[Ahmed El-Gabbas]

Thanks Ben for your suggestion.

Adding *prj = TRUE *argument caused this error:
   Error in ncdf4::ncvar_def(varname, varunit, list(xdim, ydim),
NAvalue(x),  :
   unused argument (prj = T)

Ahmed

On Tue, Feb 12, 2019 at 3:13 PM Ben Tupper  wrote:

> Hi,
>
> Have you tried using the 'prj' argument to writeRaster?  I don't know that
> it will be the solution, but according to ...
>
>
> https://www.rdocumentation.org/packages/raster/versions/2.8-19/topics/writeRaster
>
> It is documented among the '...' arguments.  Setting prj = TRUE will cause
> the crs to be written to the file.
>
> Cheers,
> Ben
>
> On Feb 12, 2019, at 8:34 AM, Ahmed El-Gabbas  wrote:
>
> Hello all,
>
> I am having a problem while converting a raster object as NetCDF (.nc)
> file, with keeping the CRS information in the output file.
> Here is a reproducible code:
>
> require(raster)
> require(ncdf4)
> CurrTemp <- tempfile()
> download.file(url = 
> "https://seaice.uni-bremen.de/data/amsre/asi_daygrid_swath/s6250/2003/feb/Antarctic/asi-s6250-20030214-v5.hdf";,
>  destfile = CurrTemp, mode = "wb", quiet = T)
> r <- raster(CurrTemp)
> r <- flip(r,2)
> extent(r) <- c(-395, 395, -395, 435)
> crs(r) <- "+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 
> +a=6378273 +b=6356889.449 +units=m +no_defs"
> r# class   : RasterLayer # dimensions  : 1328, 1264, 1678592  (nrow, 
> ncol, ncell)# resolution  : 6250, 6250  (x, y)# extent  : -395, 
> 395, -395, 435  (xmin, xmax, ymin, ymax)# coord. ref. : 
> +proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273 
> +b=6356889.449 +units=m +no_defs # data source : in memory# names   : 
> layer # values  : 0, 100  (min, max)
>
> So far the raster object reads well, with CRS information included.
> However, when I try to save it as .nc file, R prints "coord. ref. : NA",
> and the produced file does not contain the CRS information.
>
> writeRaster(r, filename = "O:/Ahmed/r001.nc", varname="IceConc",
> overwrite=TRUE, format="CDF",
> xname="Longitude", yname="Latitude")# class   : RasterLayer # 
> dimensions  : 1328, 1264, 1678592  (nrow, ncol, ncell)# resolution  : 6250, 
> 6250  (x, y)# extent  : -395, 395, -395, 435  (xmin, 
> xmax, ymin, ymax)# coord. ref. : NA # data source : O:/Ahmed/r001.nc # names  
>  : IceConc # zvar: IceConc
>
> raster("O:/Ahmed/r001.nc")# class   : RasterLayer # dimensions  : 1328, 
> 1264, 1678592  (nrow, ncol, ncell)# resolution  : 6250, 6250  (x, y)# extent  
> : -395, 395, -395, 435  (xmin, xmax, ymin, ymax)# coord. 
> ref. : NA # data source : O:/Ahmed/r001.nc # names   : IceConc # zvar 
>: IceConc
>
> Any solution?
>
> N.B. I also sent the same question at stackoverflow
> ,
> without a solution so far.
>
> Thanks
>
> ~
> *Dr. Ahmed El-Gabbas,*
> Ocean Acoustics Lab,
> Alfred-Wegener-Institut
> Helmholtz-Zentrum für Polar und Meeresforschung
> 
> My Website:  https://elgabbas.github.io
> ___
> 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
>
> Ecological Forecasting: https://eco.bigelow.org/
>
>
>
>
>
>

[[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] converting raster to netcdf .nc with keeping CRS information

[Ben Tupper]

Hi Ahmed,

I can confirm that on my system I get the same result as you.

Does your workflow confine you to using NetCDF?  I can successfully write/read 
to the native raster format.

r2 <- writeRaster(r, filename = "r001.grd", overwrite=TRUE)
r2 
#class   : RasterLayer 
#dimensions  : 1328, 1264, 1678592  (nrow, ncol, ncell)
#resolution  : 6250, 6250  (x, y)
#extent  : -395, 395, -395, 435  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 
+a=6378273 +b=6356889.449 +units=m +no_defs 
#data source : /home/btupper/r001.grd 
#names   : layer 
#values  : 0, 100  (min, max)


r3 <- raster("r001.grd")
r3 
#class   : RasterLayer 
#dimensions  : 1328, 1264, 1678592  (nrow, ncol, ncell)
#resolution  : 6250, 6250  (x, y)
#extent  : -395, 395, -395, 435  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 
+a=6378273 +b=6356889.449 +units=m +no_defs 
#data source : /home/btupper/r001.grd 
#names   : layer 
#values  : 0, 100  (min, max)

Ben

> sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so

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

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

other attached packages:
[1] raster_2.8-19 sp_1.3-1  ncdf4_1.16   

loaded via a namespace (and not attached):
[1] compiler_3.5.1   rgdal_1.3-6  Rcpp_1.0.0   codetools_0.2-15
[5] grid_3.5.1   lattice_0.20-35 





> On Feb 12, 2019, at 9:19 AM, Ahmed El-Gabbas  wrote:
> 
> Thanks Ben for your suggestion.
> 
> Adding *prj = TRUE *argument caused this error:
>   Error in ncdf4::ncvar_def(varname, varunit, list(xdim, ydim),
> NAvalue(x),  :
>   unused argument (prj = T)
> 
> Ahmed
> 
> On Tue, Feb 12, 2019 at 3:13 PM Ben Tupper  wrote:
> 
>> Hi,
>> 
>> Have you tried using the 'prj' argument to writeRaster?  I don't know that
>> it will be the solution, but according to ...
>> 
>> 
>> https://www.rdocumentation.org/packages/raster/versions/2.8-19/topics/writeRaster
>> 
>> It is documented among the '...' arguments.  Setting prj = TRUE will cause
>> the crs to be written to the file.
>> 
>> Cheers,
>> Ben
>> 
>> On Feb 12, 2019, at 8:34 AM, Ahmed El-Gabbas  wrote:
>> 
>> Hello all,
>> 
>> I am having a problem while converting a raster object as NetCDF (.nc)
>> file, with keeping the CRS information in the output file.
>> Here is a reproducible code:
>> 
>> require(raster)
>> require(ncdf4)
>> CurrTemp <- tempfile()
>> download.file(url = 
>> "https://seaice.uni-bremen.de/data/amsre/asi_daygrid_swath/s6250/2003/feb/Antarctic/asi-s6250-20030214-v5.hdf";,
>>  destfile = CurrTemp, mode = "wb", quiet = T)
>> r <- raster(CurrTemp)
>> r <- flip(r,2)
>> extent(r) <- c(-395, 395, -395, 435)
>> crs(r) <- "+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 
>> +a=6378273 +b=6356889.449 +units=m +no_defs"
>> r# class   : RasterLayer # dimensions  : 1328, 1264, 1678592  (nrow, 
>> ncol, ncell)# resolution  : 6250, 6250  (x, y)# extent  : -395, 
>> 395, -395, 435  (xmin, xmax, ymin, ymax)# coord. ref. : 
>> +proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273 
>> +b=6356889.449 +units=m +no_defs # data source : in memory# names : 
>> layer # values  : 0, 100  (min, max)
>> 
>> So far the raster object reads well, with CRS information included.
>> However, when I try to save it as .nc file, R prints "coord. ref. : NA",
>> and the produced file does not contain the CRS information.
>> 
>> writeRaster(r, filename = "O:/Ahmed/r001.nc", varname="IceConc",
>>overwrite=TRUE, format="CDF",
>>xname="Longitude", yname="Latitude")# class   : RasterLayer # 
>> dimensions  : 1328, 1264, 1678592  (nrow, ncol, ncell)# resolution  : 6250, 
>> 6250  (x, y)# extent  : -395, 395, -395, 435  (xmin, 
>> xmax, ymin, ymax)# coord. ref. : NA # data source : O:/Ahmed/r001.nc # names 
>>   : IceConc # zvar: IceConc
>> 
>> raster("O:/Ahmed/r001.nc")# class   : RasterLayer # dimensions  : 1328, 
>> 1264, 1678592  (nrow, ncol, ncell)# resolution  : 6250, 6250  (x, y)# extent 
>>  : -395, 395, -395, 435  (xmin, xmax, ymin, ymax)# 
>> coord. ref. : NA # data source : O:/Ahmed/r001.nc # names   : IceConc # 
>> zvar: IceConc
>> 
>> Any solution?
>> 
>> N.B. I also sent the same question at stackoverflow
>> 

Re: [R-sig-Geo] converting raster to netcdf .nc with keeping CRS information

[Michael Sumner]

I agree, I wouldn't try to bend NetCDF to work this way via raster. The
convention/s require parameters of the projection to be specified via
attributes and you could do that via RNetCDF or ncdf4, but it's a lot of
painful work.

It's not easy to create NetCDF from R, it's a simple as that - and
typically those files don't play well with projections - it's changed
recently but you'll see a lot of patchy partial support across different
softwares.

You can use rgdal or call out out to GDAL - but then round-tripping gives
warnings, because again raster has to unpack all the components of metadata
that GDAL created, and none of this is very systematic and clearly is
broken atm.

system('gdal_translate r001.nc r001_aux.nc -of NetCDF -a_srs
"+init=epsg:3412"')
## you could do this with  stars too
rgdal::writeGDAL(as(r, "SpatialGridDataFrame"), "r001_rgdal.nc", drivername
= "NetCDF")


But, beware as in neither case is raster guaranteed to get the right
interpretation (THIS IS WRONG):

raster("r001_aux.nc")
raster("r001_rgdal.nc")
class   : RasterLayer
dimensions  : 1328, 1264, 1678592  (nrow, ncol, ncell)
resolution  : 6250, 6250  (x, y)
extent  : -395, 395, -395, 435  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=stere +lon_0=0 +x_0=0 +y_0=0 +lat_0=-90 +pm=0
+a=6378273 +rf=298.279411123064
data source : r001_rgdal.nc
names   : GDAL.Band.Number.1
zvar: Band1

BTW, you called your X and Y "longitude" and "latitude" but this grid isn't
defined by those (eastings and northings is not a better name either, but
that's another problem).

Do you have to create new files?  I find that constantly a fraught problem,
For usage it's way easier to wrap stuff up with R code and work around all
these problems.

Cheers, Mike.

On Wed, 13 Feb 2019 at 12:47 Ben Tupper  wrote:

> Hi Ahmed,
>
> I can confirm that on my system I get the same result as you.
>
> Does your workflow confine you to using NetCDF?  I can successfully
> write/read to the native raster format.
>
> r2 <- writeRaster(r, filename = "r001.grd", overwrite=TRUE)
> r2
> #class   : RasterLayer
> #dimensions  : 1328, 1264, 1678592  (nrow, ncol, ncell)
> #resolution  : 6250, 6250  (x, y)
> #extent  : -395, 395, -395, 435  (xmin, xmax, ymin,
> ymax)
> #coord. ref. : +proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0
> +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs
> #data source : /home/btupper/r001.grd
> #names   : layer
> #values  : 0, 100  (min, max)
>
>
> r3 <- raster("r001.grd")
> r3
> #class   : RasterLayer
> #dimensions  : 1328, 1264, 1678592  (nrow, ncol, ncell)
> #resolution  : 6250, 6250  (x, y)
> #extent  : -395, 395, -395, 435  (xmin, xmax, ymin,
> ymax)
> #coord. ref. : +proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0
> +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs
> #data source : /home/btupper/r001.grd
> #names   : layer
> #values  : 0, 100  (min, max)
>
> Ben
>
> > sessionInfo()
> R version 3.5.1 (2018-07-02)
> Platform: x86_64-redhat-linux-gnu (64-bit)
> Running under: CentOS Linux 7 (Core)
>
> Matrix products: default
> BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
>
> locale:
>  [1] LC_CTYPE=en_US.UTF-8   LC_NUMERIC=C
>  [3] LC_TIME=en_US.UTF-8LC_COLLATE=en_US.UTF-8
>  [5] LC_MONETARY=en_US.UTF-8LC_MESSAGES=en_US.UTF-8
>  [7] LC_PAPER=en_US.UTF-8   LC_NAME=C
>  [9] LC_ADDRESS=C   LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats graphics  grDevices utils datasets  methods   base
>
> other attached packages:
> [1] raster_2.8-19 sp_1.3-1  ncdf4_1.16
>
> loaded via a namespace (and not attached):
> [1] compiler_3.5.1   rgdal_1.3-6  Rcpp_1.0.0   codetools_0.2-15
> [5] grid_3.5.1   lattice_0.20-35
>
>
>
>
>
> > On Feb 12, 2019, at 9:19 AM, Ahmed El-Gabbas  wrote:
> >
> > Thanks Ben for your suggestion.
> >
> > Adding *prj = TRUE *argument caused this error:
> >   Error in ncdf4::ncvar_def(varname, varunit, list(xdim, ydim),
> > NAvalue(x),  :
> >   unused argument (prj = T)
> >
> > Ahmed
> >
> > On Tue, Feb 12, 2019 at 3:13 PM Ben Tupper  wrote:
> >
> >> Hi,
> >>
> >> Have you tried using the 'prj' argument to writeRaster?  I don't know
> that
> >> it will be the solution, but according to ...
> >>
> >>
> >>
> https://www.rdocumentation.org/packages/raster/versions/2.8-19/topics/writeRaster
> >>
> >> It is documented among the '...' arguments.  Setting prj = TRUE will
> cause
> >> the crs to be written to the file.
> >>
> >> Cheers,
> >> Ben
> >>
> >> On Feb 12, 2019, at 8:34 AM, Ahmed El-Gabbas 
> wrote:
> >>
> >> Hello all,
> >>
> >> I am having a problem while converting a raster object as NetCDF (.nc)
> >> file, with keeping the CRS information in the output file.
> >> Here is a reproducible code:
> >>
> >> require(raster)
> >> require(ncdf4)
> >> CurrTemp <- tempfile()
> >> download.file(url = "
> https://s

Re: [R-sig-Geo] Sampling random directional lines within a polygon

[Hannah Justen]

Hi everyone,

I wanted to generate random lines between two spatial points (random points
in polygons). The lines should consist of segments (9 segments), after
following suggestions I received to my previous post, I ended up using
trajectories. My aim now is to restrict these random trajectories (see
points in the attachment as an example) to the extent of a polygon (shown
in the attachment in yellow). So the trajectories should lie within the
polygon.
Please find the code I used below:
I first generated a data frame with my start and end point and used it to
generate a trajectory (straight line in attachment). I then generated a
random trajectory and shifted it so that the start and endpoints were equal
to my first trajectory. I transformed the matrices to a Spatial Data frame
to be able to plot it. Does anyone know how I can restrict the extent of
the trajectory that I randomly generate?

#coordinates of start and end point and time steps
x<-c(1649786,-1636500)
y<-c(-2902593,738500)
times<-c(0,9)
start_end_points<-as.data.frame(cbind(x,y,times))

# generate a trajectory between start and endpoint
trj<-TrajFromCoords(point_trj,xCol=1,yCol=2,timeCol=3, spatialUnits="m",
timeUnits="d")

# generate a random trajectory with 9 segments
track<-TrajGenerate(n=9, random=T, stepLength=TrajLength(trj)/9)

# To be able to use StartEndTranslate transform trajectories into matrices
ma1<-as.matrix(trj)
ma2<-as.matrix(track)

# Translate, rotate and scale the points of traj2 using traj1. The new traj
will have the same start and end points as traj1.
p<-StartEndTranslate(ma1,ma2)

# to check if this worked transform into SpatialDataFrame
# the start and end point of the new matrix are the same as the traj1 which
is good, the rest of the points changed as well
p_df<-as.data.frame(p)

# change lat / longs to nummeric
p_df$x<-as.numeric(p_df$x)
p_df$y<-as.numeric(p_df$y)

# generate a spatial dataframe to check if the StartEndTranslate function
worked by plotting it
xy<-p_df[,c(1,2)]

spdf<-SpatialPointsDataFrame(coords=xy, data=p_df, proj4string =
CRS(projection))

Thank you very much!

On Wed, Feb 6, 2019 at 5:42 PM Barry Rowlingson 
wrote:

> Had another idea which is now implemented...
>
> Consider any segmented path of segments of lengths L_i at angles A_i. Its
> endpoint will be the vector sum of those segments, ie at (x,y) = (sum(L_i
> cos(A_i)), sum(L_i sin(A_i)).
>
> To create a segmented path to a given (x,y), solve that expression for the
> angles A_i. In R you can treat this as an optimisation problem - find a set
> of angles A_i that minimise the distance of the end of the segmented path
> from the target end point.
>
> Here's some code that does that for a path from 0,0 to 0,1:
>
> pathit <- function(segments){
> obj = function(angles){
> dxdy = dxdy(segments, angles)
> xerr = dxdy$dx-1
> yerr = dxdy$dy
> err = sqrt(xerr^2 + yerr^2)
> err
> }
> angles = runif(length(segments), 0 , 2*pi)
> optim(angles, obj)
> }
>
> dxdy = function(segments, angles){
> dx = sum(segments * cos(angles))
> dy = sum(segments * sin(angles))
> list(dx=dx, dy=dy)
> }
>
> plotsegs <- function(segments, angles){
> x = rep(NA, length(segments) +1)
> y = x
> x[1] = 0
> y[1] = 0
> for(i in 2:(length(segments)+1)){
> x[i] = x[i-1] + segments[i-1]*cos(angles[i-1])
> y[i] = y[i-1] + segments[i-1]*sin(angles[i-1])
> }
> cbind(x,y)
> }
>
> This is deliberately written naively for clarity.
>
> To use, set the segment sizes, optimise, and then plot:
>
>  s1 = c(.1,.3,.2,.1,.3,.3,.1)
>  a1 = pathit(s1)
>  plot(plotsegs(s1,a1$par),type="l")
>
> which should show a path of seven segments from 0,0 to 0,1 - since the
> initial starting values are random the model can find different solutions.
> Run again for a different path.
>
> To see what the space of paths looks like, do lots and overplot them:
>
>   lots = lapply(1:1000, function(i)plotsegs(s1,pathit(s1)$par))
>   plot(c(-.1,1.1),c(-1,1))
>   p = lapply(lots, function(xy){lines(xy)})
>
> this should show 1000 paths, and illustrates the "ellipse" of path
> possibles that I mentioned in the previous email.
>
> Sometimes the optimiser struggles to find a solution and so you should
> probably test the output from optim for convergence and to make sure the
> target function is close enough to zero for your purposes.  For the example
> above most of the time the end point is about 1e-5  from (1,0) but for
> harder problems such as s = rep(.1, 11) which only has 0.1 of extra "slack"
> length, the error can be 0.02 and failed convergence. Possibly longer optim
> runs would help or constraining the angles.
>
> Anyway, interesting problem
>
> Barry
>
>
>
>
>
>
>
> On Wed, Feb 6, 2019 at 8:23 PM Barry Rowlingson 
> wrote:
>
>> Do you want to generate these for input into some statistical process, or
>> to generate some test data that looks a bit like real data? I think
>> generating test data