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: 1 hour 41 min ago

Re: Use rbind for more than two data frame

Thu, 02/21/2019 - 06:37
Yes,

alter columname   "bio15_FP" to "bio15" in  table 'sdmdata_F'

:-)

greetings  Peter


Op 21-2-2019 om 13:29 schreef Lara Silva:
> Hello,
>
> I am trying to join different data frames using rbind but I have the
> following mensage:
>
> sdmata_global <-
> rbind(sdmdata_SM,sdmdata_P,sdmdata_T,sdmdata_F,sdmdata_SJ)
>
> Error in match.names(clabs, names(xi)) :
> names do not match previous names
>
>> class(sdmdata_SM)
> [1] "data.frame"
>> class(sdmdata_P)
> [1] "data.frame"
>> class(sdmdata_T)
> [1] "data.frame"
>> class(sdmdata_F)
> [1] "data.frame"
>> class(sdmdata_SJ)
> [1] "data.frame"
>> dim(sdmdata_SM)
> [1] 10003    22
>> dim(sdmdata_P)
> [1] 10005    22
>> dim(sdmdata_T)
> [1] 10004    22
>> dim(sdmdata_F)
> [1] 10001    22
>> dim(sdmdata_SJ)
> [1] 10001    22
>
>
>> dim(sdmdata_SM)
> [1] 10003    22
>> dim(sdmdata_SM)
> [1] 10003    22
>> dim(sdmdata_P)
> [1] 10005    22
>> dim(sdmdata_T)
> [1] 10004    22
>> dim(sdmdata_F)
> [1] 10001    22
>> dim(sdmdata_SJ)
> [1] 10001    22
>> names(sdmdata_SM)
>   [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
> [10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"  "bio6"
> [19] "bio7"  "bio8"  "bio9"  "code"
>> names(sdmdata_P)
>   [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
> [10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"  "bio6"
> [19] "bio7"  "bio8"  "bio9"  "code"
>> names(sdmdata_T)
>   [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
> [10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"  "bio6"
> [19] "bio7"  "bio8"  "bio9"  "code"
>> names(sdmdata_F)
>   [1] "pb"       "bio1"     "bio10"    "bio11"    "bio12"    "bio13"
>   [7] "bio14"    "bio15_FP" "bio16"    "bio17"    "bio18"    "bio19"
> [13] "bio2"     "bio3"     "bio4"     "bio4a"    "bio5"     "bio6"
> [19] "bio7"     "bio8"     "bio9"     "code"
>> names(sdmdata_SJ)
>   [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
> [10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"  "bio6"
> [19] "bio7"  "bio8"  "bio9"  "code"
> Any suggestion?
>
> Thanks!
>
> Lara
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Use rbind for more than two data frame

Thu, 02/21/2019 - 06:29
Hello,

I am trying to join different data frames using rbind but I have the
following mensage:

sdmata_global <-
rbind(sdmdata_SM,sdmdata_P,sdmdata_T,sdmdata_F,sdmdata_SJ)

Error in match.names(clabs, names(xi)) :
names do not match previous names

> class(sdmdata_SM)
[1] "data.frame"
> class(sdmdata_P)
[1] "data.frame"
> class(sdmdata_T)
[1] "data.frame"
> class(sdmdata_F)
[1] "data.frame"
> class(sdmdata_SJ)
[1] "data.frame"
> dim(sdmdata_SM)
[1] 10003    22
> dim(sdmdata_P)
[1] 10005    22
> dim(sdmdata_T)
[1] 10004    22
> dim(sdmdata_F)
[1] 10001    22
> dim(sdmdata_SJ)
[1] 10001    22


> dim(sdmdata_SM)
[1] 10003    22
> dim(sdmdata_SM)
[1] 10003    22
> dim(sdmdata_P)
[1] 10005    22
> dim(sdmdata_T)
[1] 10004    22
> dim(sdmdata_F)
[1] 10001    22
> dim(sdmdata_SJ)
[1] 10001    22
>

> names(sdmdata_SM)
 [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
[10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"  "bio6"
[19] "bio7"  "bio8"  "bio9"  "code"
> names(sdmdata_P)
 [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
[10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"  "bio6"
[19] "bio7"  "bio8"  "bio9"  "code"
> names(sdmdata_T)
 [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
[10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"  "bio6"
[19] "bio7"  "bio8"  "bio9"  "code"
> names(sdmdata_F)
 [1] "pb"       "bio1"     "bio10"    "bio11"    "bio12"    "bio13"
 [7] "bio14"    "bio15_FP" "bio16"    "bio17"    "bio18"    "bio19"
[13] "bio2"     "bio3"     "bio4"     "bio4a"    "bio5"     "bio6"
[19] "bio7"     "bio8"     "bio9"     "code"
> names(sdmdata_SJ)
 [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
[10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"  "bio6"
[19] "bio7"  "bio8"  "bio9"  "code"
>

Any suggestion?

Thanks!

Lara

        [[alternative HTML version deleted]]

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

Re: Reverse color in stplot

Thu, 02/21/2019 - 01:45
I think that it accepts an argument col.regions where you specify a
color palette; if you have a particular one, called pal, and want to
revert it, use rev(pal).

On 2/21/19 7:02 AM, Mahsa Sameti wrote:
> Dear all
>
> I want to reverse color in my spatio-temporal maps that are created with
> stplot function. How can i solve this problem?
>
> Thanks alot
> Mahsa Sameti
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

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

pEpkey.asc (2K) Download Attachment

Reverse color in stplot

Wed, 02/20/2019 - 23:54
Dear all

I want to reverse color in my spatio-temporal maps that are created with
stplot function. How can i solve this problem?

Thanks alot
Mahsa Sameti

        [[alternative HTML version deleted]]

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

Re: [R-sig-eco] Problem in generate random samples in r

Wed, 02/20/2019 - 05:29
Thanks!

Lara

Henrik Eckermann <[hidden email]> escreveu no dia quarta,
20/02/2019 à(s) 10:25:

> Hi Lara,
>
> the error is quite informative. You use the sample function twice. At
> least for one the length of the vector where you sample from is shorter
> (contains less values) than than the number of samples you wanna draw.
>
> best,
>
> Henrik
>
> > On 20. Feb 2019, at 12:20, Lara Silva <[hidden email]> wrote:
> >
> > Hello,
> >
> > I am trying to generate random samples from the following code
> >
> > ### Setting random seed to always create the same
> > ### Random set of points
> > set.seed(0)
> >
> > absences_15000<-absences[sample(nrow(absences), 15000),]
> > points(absences_15000, cex=0.1)
> >
> > ## Subsample_10000
> > set.seed(0)
> >
> > absences_10000<-absences_15000[sample(nrow(absences_15000), 10000),]
> > dim(absences_10000)
> >
> > I get the following error:
> >
> > Error in sample.int(length(x), size, replace, prob) :
> >  cannot take a sample larger than the population when 'replace = FALSE'
> >
> > Any advice?
> >
> > Regards,
> >
> > Lara
> >
> >       [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-ecology mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>
>
        [[alternative HTML version deleted]]

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

Re: [R-sig-eco] Problem in generate random samples in r

Wed, 02/20/2019 - 05:28
My guess would be that absences_150000 doesn't have 10000 rows. Can you
confirm or refute this?

Cheers,
Roman

On Wed, Feb 20, 2019 at 12:20 PM Lara Silva <[hidden email]>
wrote:

> Hello,
>
> I am trying to generate random samples from the following code
>
> ### Setting random seed to always create the same
> ### Random set of points
> set.seed(0)
>
> absences_15000<-absences[sample(nrow(absences), 15000),]
> points(absences_15000, cex=0.1)
>
> ## Subsample_10000
> set.seed(0)
>
> absences_10000<-absences_15000[sample(nrow(absences_15000), 10000),]
> dim(absences_10000)
>
> I get the following error:
>
> Error in sample.int(length(x), size, replace, prob) :
>   cannot take a sample larger than the population when 'replace = FALSE'
>
> Any advice?
>
> Regards,
>
> Lara
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-ecology mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>

--
In God we trust, all others bring data.

        [[alternative HTML version deleted]]

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

Re: Plot coordinates on world map with Robinson CRS - ggplot2

Mon, 02/18/2019 - 05:43
Thank you very much. This works perfectly. Appreciate it.

Sincerely,

Milu

On Mon, Feb 18, 2019 at 11:01 AM Michael Sumner <[hidden email]> wrote:

> You are very close, just needed the expert-knowledge for the incantation
> required to transform the longlat points (as was already done for the
> prj.coord, btw).
>
> HTH: https://gist.github.com/mdsumner/afd9228aa1dc4652f8bf193cad4245aa
>
> It does come as a surprise to many that R's plotting is not unit or
> projection aware (well, it is sometimes, but not usually ... a modern
> example is ggplot2::geom_sf  and the sf package but you need quite a lot of
> prior-won expertise there too, good luck).
>
> Cheers, Mike.
>
> On Mon, 18 Feb 2019 at 20:01 Miluji Sb <[hidden email]> wrote:
>
>> Dear all,
>>
>> I am trying to plot coordinates on a world map with Robinson CRS. While
>> the
>> world map is generated without any issues, when I try to plot the points -
>> I only get a single point.
>>
>> The code I am using and the coordinates data is below. What am I doing
>> wrong? Any help/suggestions will be highly appreciated.
>>
>> library(data.table)
>> library(foreign)
>> library(readstata13)
>> library(rgdal)
>> library(maptools)
>> library(ggplot2)
>> library(dplyr)
>>
>> load(url("
>>
>> https://github.com/valentinitnelav/RandomScripts/blob/master/NaturalEarth.RData?raw=true
>> "))
>>
>> PROJ <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84
>> +units=m +no_defs"
>>
>> NE_countries_rob  <- spTransform(NE_countries, CRSobj = PROJ)
>> NE_graticules_rob <- spTransform(NE_graticules, CRSobj = PROJ)
>> NE_box_rob        <- spTransform(NE_box, CRSobj = PROJ)
>>
>> # project long-lat coordinates for graticule label data frames
>> # (two extra columns with projected XY are created)
>> prj.coord <- project(cbind(lbl.Y$lon, lbl.Y$lat), proj=PROJ)
>> lbl.Y.prj <- cbind(prj.coord, lbl.Y)
>> names(lbl.Y.prj)[1:2] <- c("X.prj","Y.prj")
>>
>> prj.coord <- project(cbind(lbl.X$lon, lbl.X$lat), proj=PROJ)
>> lbl.X.prj <- cbind(prj.coord, lbl.X)
>> names(lbl.X.prj)[1:2] <- c("X.prj","Y.prj")
>>
>> m <- ggplot() +
>>   # add Natural Earth countries projected to Robinson, give black border
>> and fill with gray
>>   geom_polygon(data=NE_countries_rob, aes(long,lat, group=group),
>> colour="black", fill="gray80", size = 0.25) +
>>   # Note: "Regions defined for each Polygons" warning has to do with
>> fortify transformation. Might get deprecated in future!
>>   # alternatively, use use map_data(NE_countries) to transform to data
>> frame and then use project() to change to desired projection.
>>   # add Natural Earth box projected to Robinson
>>   geom_polygon(data=NE_box_rob, aes(x=long, y=lat), colour="black",
>> fill="transparent", size = 0.25) +
>>   # add graticules projected to Robinson
>>   geom_path(data=NE_graticules_rob, aes(long, lat, group=group),
>> linetype="dotted", color="grey50", size = 0.25) +
>>   # add graticule labels - latitude and longitude
>>   geom_text(data = lbl.Y.prj, aes(x = X.prj, y = Y.prj, label = lbl),
>> color="grey50", size=2) +
>>   geom_text(data = lbl.X.prj, aes(x = X.prj, y = Y.prj, label = lbl),
>> color="grey50", size=2) +
>>   # the default, ratio = 1 in coord_fixed ensures that one unit on the
>> x-axis is the same length as one unit on the y-axis
>>   coord_fixed(ratio = 1) +
>>   geom_point(data=df,
>>              aes(x=lon, y=lat), colour="Deep Pink",
>>              fill="Pink",pch=21, size=2, alpha=I(0.7))
>>   # remove the background and default gridlines
>>   theme_void()
>>
>> ## coordinates dataframe
>> dput(df)
>> structure(list(lon = c(2.67569724621467, 17.5766416259819,
>> 28.4126232192772,
>> 23.8147674538232, 29.8917589327105), lat = c(28.1503115976162,
>> -12.3388787380201, 9.78891068739477, -22.1873831176644, -3.36546931479253
>> )), class = "data.frame", row.names = c(NA, -5L))
>>
>> Sincerely,
>>
>> Milu
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
> --
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> Kingston Tasmania 7050 Australia
>
>
        [[alternative HTML version deleted]]

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

Re: Plot coordinates on world map with Robinson CRS - ggplot2

Mon, 02/18/2019 - 04:00
You are very close, just needed the expert-knowledge for the incantation
required to transform the longlat points (as was already done for the
prj.coord, btw).

HTH: https://gist.github.com/mdsumner/afd9228aa1dc4652f8bf193cad4245aa

It does come as a surprise to many that R's plotting is not unit or
projection aware (well, it is sometimes, but not usually ... a modern
example is ggplot2::geom_sf  and the sf package but you need quite a lot of
prior-won expertise there too, good luck).

Cheers, Mike.

On Mon, 18 Feb 2019 at 20:01 Miluji Sb <[hidden email]> wrote:

> Dear all,
>
> I am trying to plot coordinates on a world map with Robinson CRS. While the
> world map is generated without any issues, when I try to plot the points -
> I only get a single point.
>
> The code I am using and the coordinates data is below. What am I doing
> wrong? Any help/suggestions will be highly appreciated.
>
> library(data.table)
> library(foreign)
> library(readstata13)
> library(rgdal)
> library(maptools)
> library(ggplot2)
> library(dplyr)
>
> load(url("
>
> https://github.com/valentinitnelav/RandomScripts/blob/master/NaturalEarth.RData?raw=true
> "))
>
> PROJ <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84
> +units=m +no_defs"
>
> NE_countries_rob  <- spTransform(NE_countries, CRSobj = PROJ)
> NE_graticules_rob <- spTransform(NE_graticules, CRSobj = PROJ)
> NE_box_rob        <- spTransform(NE_box, CRSobj = PROJ)
>
> # project long-lat coordinates for graticule label data frames
> # (two extra columns with projected XY are created)
> prj.coord <- project(cbind(lbl.Y$lon, lbl.Y$lat), proj=PROJ)
> lbl.Y.prj <- cbind(prj.coord, lbl.Y)
> names(lbl.Y.prj)[1:2] <- c("X.prj","Y.prj")
>
> prj.coord <- project(cbind(lbl.X$lon, lbl.X$lat), proj=PROJ)
> lbl.X.prj <- cbind(prj.coord, lbl.X)
> names(lbl.X.prj)[1:2] <- c("X.prj","Y.prj")
>
> m <- ggplot() +
>   # add Natural Earth countries projected to Robinson, give black border
> and fill with gray
>   geom_polygon(data=NE_countries_rob, aes(long,lat, group=group),
> colour="black", fill="gray80", size = 0.25) +
>   # Note: "Regions defined for each Polygons" warning has to do with
> fortify transformation. Might get deprecated in future!
>   # alternatively, use use map_data(NE_countries) to transform to data
> frame and then use project() to change to desired projection.
>   # add Natural Earth box projected to Robinson
>   geom_polygon(data=NE_box_rob, aes(x=long, y=lat), colour="black",
> fill="transparent", size = 0.25) +
>   # add graticules projected to Robinson
>   geom_path(data=NE_graticules_rob, aes(long, lat, group=group),
> linetype="dotted", color="grey50", size = 0.25) +
>   # add graticule labels - latitude and longitude
>   geom_text(data = lbl.Y.prj, aes(x = X.prj, y = Y.prj, label = lbl),
> color="grey50", size=2) +
>   geom_text(data = lbl.X.prj, aes(x = X.prj, y = Y.prj, label = lbl),
> color="grey50", size=2) +
>   # the default, ratio = 1 in coord_fixed ensures that one unit on the
> x-axis is the same length as one unit on the y-axis
>   coord_fixed(ratio = 1) +
>   geom_point(data=df,
>              aes(x=lon, y=lat), colour="Deep Pink",
>              fill="Pink",pch=21, size=2, alpha=I(0.7))
>   # remove the background and default gridlines
>   theme_void()
>
> ## coordinates dataframe
> dput(df)
> structure(list(lon = c(2.67569724621467, 17.5766416259819,
> 28.4126232192772,
> 23.8147674538232, 29.8917589327105), lat = c(28.1503115976162,
> -12.3388787380201, 9.78891068739477, -22.1873831176644, -3.36546931479253
> )), class = "data.frame", row.names = c(NA, -5L))
>
> Sincerely,
>
> Milu
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

        [[alternative HTML version deleted]]

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

Re: Plot coordinates on world map with Robinson CRS - ggplot2

Mon, 02/18/2019 - 03:17
On Mon, 18 Feb 2019, Miluji Sb wrote:

> Dear all,
>
> I am trying to plot coordinates on a world map with Robinson CRS. While the
> world map is generated without any issues, when I try to plot the points -
> I only get a single point.

Well, the point is correct, since you didn't project df, did you? This
seems to be typical of the arcane nature of ggplot invocations. Did you
try tmap? Almost all the loaded&attached packages are superfluous too.

Roger

>
> The code I am using and the coordinates data is below. What am I doing
> wrong? Any help/suggestions will be highly appreciated.
>
> library(data.table)
> library(foreign)
> library(readstata13)
> library(rgdal)
> library(maptools)
> library(ggplot2)
> library(dplyr)
>
> load(url("
> https://github.com/valentinitnelav/RandomScripts/blob/master/NaturalEarth.RData?raw=true
> "))
>
> PROJ <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84
> +units=m +no_defs"
>
> NE_countries_rob  <- spTransform(NE_countries, CRSobj = PROJ)
> NE_graticules_rob <- spTransform(NE_graticules, CRSobj = PROJ)
> NE_box_rob        <- spTransform(NE_box, CRSobj = PROJ)
>
> # project long-lat coordinates for graticule label data frames
> # (two extra columns with projected XY are created)
> prj.coord <- project(cbind(lbl.Y$lon, lbl.Y$lat), proj=PROJ)
> lbl.Y.prj <- cbind(prj.coord, lbl.Y)
> names(lbl.Y.prj)[1:2] <- c("X.prj","Y.prj")
>
> prj.coord <- project(cbind(lbl.X$lon, lbl.X$lat), proj=PROJ)
> lbl.X.prj <- cbind(prj.coord, lbl.X)
> names(lbl.X.prj)[1:2] <- c("X.prj","Y.prj")
>
> m <- ggplot() +
>  # add Natural Earth countries projected to Robinson, give black border
> and fill with gray
>  geom_polygon(data=NE_countries_rob, aes(long,lat, group=group),
> colour="black", fill="gray80", size = 0.25) +
>  # Note: "Regions defined for each Polygons" warning has to do with
> fortify transformation. Might get deprecated in future!
>  # alternatively, use use map_data(NE_countries) to transform to data
> frame and then use project() to change to desired projection.
>  # add Natural Earth box projected to Robinson
>  geom_polygon(data=NE_box_rob, aes(x=long, y=lat), colour="black",
> fill="transparent", size = 0.25) +
>  # add graticules projected to Robinson
>  geom_path(data=NE_graticules_rob, aes(long, lat, group=group),
> linetype="dotted", color="grey50", size = 0.25) +
>  # add graticule labels - latitude and longitude
>  geom_text(data = lbl.Y.prj, aes(x = X.prj, y = Y.prj, label = lbl),
> color="grey50", size=2) +
>  geom_text(data = lbl.X.prj, aes(x = X.prj, y = Y.prj, label = lbl),
> color="grey50", size=2) +
>  # the default, ratio = 1 in coord_fixed ensures that one unit on the
> x-axis is the same length as one unit on the y-axis
>  coord_fixed(ratio = 1) +
>  geom_point(data=df,
>             aes(x=lon, y=lat), colour="Deep Pink",
>             fill="Pink",pch=21, size=2, alpha=I(0.7))
>  # remove the background and default gridlines
>  theme_void()
>
> ## coordinates dataframe
> dput(df)
> structure(list(lon = c(2.67569724621467, 17.5766416259819,
> 28.4126232192772,
> 23.8147674538232, 29.8917589327105), lat = c(28.1503115976162,
> -12.3388787380201, 9.78891068739477, -22.1873831176644, -3.36546931479253
> )), class = "data.frame", row.names = c(NA, -5L))
>
> Sincerely,
>
> Milu
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> 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

Plot coordinates on world map with Robinson CRS - ggplot2

Mon, 02/18/2019 - 02:57
Dear all,

I am trying to plot coordinates on a world map with Robinson CRS. While the
world map is generated without any issues, when I try to plot the points -
I only get a single point.

The code I am using and the coordinates data is below. What am I doing
wrong? Any help/suggestions will be highly appreciated.

library(data.table)
library(foreign)
library(readstata13)
library(rgdal)
library(maptools)
library(ggplot2)
library(dplyr)

load(url("
https://github.com/valentinitnelav/RandomScripts/blob/master/NaturalEarth.RData?raw=true
"))

PROJ <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84
+units=m +no_defs"

NE_countries_rob  <- spTransform(NE_countries, CRSobj = PROJ)
NE_graticules_rob <- spTransform(NE_graticules, CRSobj = PROJ)
NE_box_rob        <- spTransform(NE_box, CRSobj = PROJ)

# project long-lat coordinates for graticule label data frames
# (two extra columns with projected XY are created)
prj.coord <- project(cbind(lbl.Y$lon, lbl.Y$lat), proj=PROJ)
lbl.Y.prj <- cbind(prj.coord, lbl.Y)
names(lbl.Y.prj)[1:2] <- c("X.prj","Y.prj")

prj.coord <- project(cbind(lbl.X$lon, lbl.X$lat), proj=PROJ)
lbl.X.prj <- cbind(prj.coord, lbl.X)
names(lbl.X.prj)[1:2] <- c("X.prj","Y.prj")

m <- ggplot() +
  # add Natural Earth countries projected to Robinson, give black border
and fill with gray
  geom_polygon(data=NE_countries_rob, aes(long,lat, group=group),
colour="black", fill="gray80", size = 0.25) +
  # Note: "Regions defined for each Polygons" warning has to do with
fortify transformation. Might get deprecated in future!
  # alternatively, use use map_data(NE_countries) to transform to data
frame and then use project() to change to desired projection.
  # add Natural Earth box projected to Robinson
  geom_polygon(data=NE_box_rob, aes(x=long, y=lat), colour="black",
fill="transparent", size = 0.25) +
  # add graticules projected to Robinson
  geom_path(data=NE_graticules_rob, aes(long, lat, group=group),
linetype="dotted", color="grey50", size = 0.25) +
  # add graticule labels - latitude and longitude
  geom_text(data = lbl.Y.prj, aes(x = X.prj, y = Y.prj, label = lbl),
color="grey50", size=2) +
  geom_text(data = lbl.X.prj, aes(x = X.prj, y = Y.prj, label = lbl),
color="grey50", size=2) +
  # the default, ratio = 1 in coord_fixed ensures that one unit on the
x-axis is the same length as one unit on the y-axis
  coord_fixed(ratio = 1) +
  geom_point(data=df,
             aes(x=lon, y=lat), colour="Deep Pink",
             fill="Pink",pch=21, size=2, alpha=I(0.7))
  # remove the background and default gridlines
  theme_void()

## coordinates dataframe
dput(df)
structure(list(lon = c(2.67569724621467, 17.5766416259819,
28.4126232192772,
23.8147674538232, 29.8917589327105), lat = c(28.1503115976162,
-12.3388787380201, 9.78891068739477, -22.1873831176644, -3.36546931479253
)), class = "data.frame", row.names = c(NA, -5L))

Sincerely,

Milu

        [[alternative HTML version deleted]]

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

Euclidean distance raster

Sun, 02/17/2019 - 20:46
Hi all,

I want to create a nearest distance raster in R. As far as I know, the
common practice is to compute the distance matrix of all the points in the
mask raster to all the features in the vector layer, then calculate the min
distance for each point. But this makes a huge (m by n) matrix and it
involves a lot of computation which is too slow.

To go around the computation, I wrote the following function that first
gets the index of the nearest feature to each point in the mask raster,
then only calculates the distance between each point and its nearest
feature. This solution is 3 times faster than the previous one, but still
not fast enough.

I wonder if anyone knows a better solution to this problem.

Thank you in advance,
Roozbeh Valavi

rasterDistance <- function(feature, rastermask){
  require(raster)
  require(sf)
  require(progress)
  p <- st_as_sf(rasterToPoints(rastermask, spatial = TRUE))
  p$indx <- st_nearest_feature(st_geometry(p), st_geometry(feature))
  pb <- progress::progress_bar$new(format = " Progress [:bar] :percent in
:elapsed",
                                   total=nrow(p), clear=FALSE, width=75) #
add progress bar
  for(i in 1:nrow(p)){
    p$dist[i] <- st_distance(st_geometry(p[i,]),
st_geometry(feature[p$indx[i],]))
    pb$tick() # update progress bar
  }
  output <- raster::rasterize(p, rastermask, field = "dist")
  return(output)
}

--
*Roozbeh Valavi*
PhD Candidate
The Quantitative & Applied Ecology Group <http://qaeco.com/>
School of BioSciences | Faculty of Science
The University of Melbourne, VIC 3010, Australia
Mobile: +61 423 283 238

        [[alternative HTML version deleted]]

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

Error with running bfastSpatial on MODIS data

Sun, 02/17/2019 - 06:35
Dear all,
I got error with running bfastSpatial on MODIS data for period of 2001-2019
. What could be my problem here?

###########

library(bfastSpatial)



##########

#MODIS data for period of 2001-2019

  n <- timeStackMODIS(list.files(path = dir, pattern = ".tif"))

   t <- system.time(

    bfm <- bfmSpatial(n, start=c(2009, 1), order=1)

 )



 # extract change raster

 change <- raster(bfm, 1)

months <- changeMonth(change)



###### Error

Error in setValues(r, values = methods::callGeneric(getValues(e1), e2)) :

  length(values) is not equal to ncell(x), or to 1

        [[alternative HTML version deleted]]

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

Re: Is there a way to merge different stacks?

Wed, 02/13/2019 - 09:28
Thank you for your answer.

Lara

Hugo Costa <[hidden email]> escreveu no dia quarta, 13/02/2019 à(s)
14:21:

> You are possibly looking for function mosaic (raster package).
> Hugo
>
> Lara Silva <[hidden email]> escreveu no dia quarta, 13/02/2019
> à(s) 14:47:
>
>> Hello everybody,
>> I am having a problem to merge several stacks. My stacks have different
>> extensions ... It is a problem.
>>
>> If possible, I want to merge several (5 stacks) asc grids into in a
>> "global
>> grid/stack"
>> How can I do it? "Merging", "joining" ...?
>> Any sugestion?
>>
>> Thanks,
>>
>> Lara Silva
>>
>> My outputs:
>>
>> 1st stack
>>
>> > grids_climatic_SMP <- list.files(pattern='asc', full.names=T)# All "asc"
>> files in folder
>> > climatic_SMP<- stack(grids_climatic_SMP)# Make a "stack" of EGV #Install
>> "stack"
>> > climatic_SMP
>> class       : RasterStack
>> dimensions  : 901, 901, 811801, 20  (nrow, ncol, ncell, nlayers)
>> resolution  : 100, 100  (x, y)
>> extent      : 586974, 677074, 4140847, 4230947  (xmin, xmax, ymin, ymax)
>> coord. ref. : NA
>> names       : bio1_SMP, bio10_SMP, bio11_SMP, bio12_SMP, bio13_SMP,
>> bio14_SMP, bio15_SMP, bio16_SMP, bio17_SMP, bio18_SMP, bio19_SMP,
>> bio2_SMP,
>> bio3_SMP, bio4_SMP, bio4a_SMP, ...
>>
>> 2nd stack
>>
>> > grids_climatic_PP <- list.files(pattern='asc', full.names=T)# All "asc"
>> files in folder
>> > climatic_PP<- stack(grids_climatic_PP)# Make a "stack" of EGV #Install
>> "stack"
>> > climatic_PP
>> class       : RasterStack
>> dimensions  : 601, 601, 361201, 20  (nrow, ncol, ncell, nlayers)
>> resolution  : 100, 100  (x, y)
>> extent      : 358603, 418703, 4229376, 4289476  (xmin, xmax, ymin, ymax)
>> coord. ref. : NA
>> names       : bio1_PP, bio10_PP, bio11_PP, bio12_PP, bio13_PP, bio14_PP,
>> bio15_PP, bio16_PP, bio17_PP, bio18_PP, bio19_PP, bio2_PP, bio3_PP,
>> bio4_PP, bio4a_PP, ...
>>
>> 3rd stack
>>
>> > grids_climatic_TP <- list.files(pattern='asc', full.names=T)# All "asc"
>> files in folder
>> > climatic_TP<- stack(grids_climatic_TP)# Make a "stack" of EGV #Install
>> "stack"
>> > climatic_TP
>> class       : RasterStack
>> dimensions  : 277, 412, 114124, 20  (nrow, ncol, ncell, nlayers)
>> resolution  : 100, 100  (x, y)
>> extent      : 461214, 502414, 4272892, 4300592  (xmin, xmax, ymin, ymax)
>> coord. ref. : NA
>> names       : bio1_TP, bio10_TP, bio11_TP, bio12_TP, bio13_TP, bio14_TP,
>> bio15_TP, bio16_TP, bio17_TP, bio18_TP, bio19_TP, bio2_TP, bio3_TP,
>> bio4_TP, bio4a_TP, ...
>>
>>
>> 4th stack
>>
>> > grids_climatic_FP <- list.files(pattern='asc', full.names=T)# All "asc"
>> files in folder
>> > climatic_FP<- stack(grids_climatic_FP)# Make a "stack" of EGV #Install
>> "stack"
>> > climatic_FP
>> class       : RasterStack
>> dimensions  : 301, 301, 90601, 20  (nrow, ncol, ncell, nlayers)
>> resolution  : 100, 100  (x, y)
>> extent      : 336088, 366188, 4256440, 4286540  (xmin, xmax, ymin, ymax)
>> coord. ref. : NA
>> names       : bio1_FP, bio10_FP, bio11_FP, bio12_FP, bio13_FP, bio14_FP,
>> bio15_FP, bio16_FP, bio17_FP, bio18_FP, bio19_FP, bio2_FP, bio3_FP,
>> bio4_FP, bio4a_FP, ...
>>
>>
>> 5th stack
>>
>> > grids_climatic_SJP <- list.files(pattern='asc', full.names=T)# All "asc"
>> files in folder
>> > climatic_SJP<- stack(grids_climatic_SJP)# Make a "stack" of EGV #Install
>> "stack"
>> > climatic_SJP
>> class       : RasterStack
>> dimensions  : 602, 602, 362404, 20  (nrow, ncol, ncell, nlayers)
>> resolution  : 100, 100  (x, y)
>> extent      : 381100, 441300, 4249487, 4309687  (xmin, xmax, ymin, ymax)
>> coord. ref. : NA
>> names       : bio1_SJP, bio10_SJP, bio11_SJP, bio12_SJP, bio13_SJP,
>> bio14_SJP, bio15_SJP, bio16_SJP, bio17_SJP, bio18_SJP, bio19_SJP,
>> bio2_SJP,
>> bio3_SJP, bio4_SJP, bio4a_SJP, ...
>>
>> >
>>
>>         [[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: Is there a way to merge different stacks?

Wed, 02/13/2019 - 09:21
You are possibly looking for function mosaic (raster package).
Hugo

Lara Silva <[hidden email]> escreveu no dia quarta, 13/02/2019
à(s) 14:47:

> Hello everybody,
> I am having a problem to merge several stacks. My stacks have different
> extensions ... It is a problem.
>
> If possible, I want to merge several (5 stacks) asc grids into in a "global
> grid/stack"
> How can I do it? "Merging", "joining" ...?
> Any sugestion?
>
> Thanks,
>
> Lara Silva
>
> My outputs:
>
> 1st stack
>
> > grids_climatic_SMP <- list.files(pattern='asc', full.names=T)# All "asc"
> files in folder
> > climatic_SMP<- stack(grids_climatic_SMP)# Make a "stack" of EGV #Install
> "stack"
> > climatic_SMP
> class       : RasterStack
> dimensions  : 901, 901, 811801, 20  (nrow, ncol, ncell, nlayers)
> resolution  : 100, 100  (x, y)
> extent      : 586974, 677074, 4140847, 4230947  (xmin, xmax, ymin, ymax)
> coord. ref. : NA
> names       : bio1_SMP, bio10_SMP, bio11_SMP, bio12_SMP, bio13_SMP,
> bio14_SMP, bio15_SMP, bio16_SMP, bio17_SMP, bio18_SMP, bio19_SMP, bio2_SMP,
> bio3_SMP, bio4_SMP, bio4a_SMP, ...
>
> 2nd stack
>
> > grids_climatic_PP <- list.files(pattern='asc', full.names=T)# All "asc"
> files in folder
> > climatic_PP<- stack(grids_climatic_PP)# Make a "stack" of EGV #Install
> "stack"
> > climatic_PP
> class       : RasterStack
> dimensions  : 601, 601, 361201, 20  (nrow, ncol, ncell, nlayers)
> resolution  : 100, 100  (x, y)
> extent      : 358603, 418703, 4229376, 4289476  (xmin, xmax, ymin, ymax)
> coord. ref. : NA
> names       : bio1_PP, bio10_PP, bio11_PP, bio12_PP, bio13_PP, bio14_PP,
> bio15_PP, bio16_PP, bio17_PP, bio18_PP, bio19_PP, bio2_PP, bio3_PP,
> bio4_PP, bio4a_PP, ...
>
> 3rd stack
>
> > grids_climatic_TP <- list.files(pattern='asc', full.names=T)# All "asc"
> files in folder
> > climatic_TP<- stack(grids_climatic_TP)# Make a "stack" of EGV #Install
> "stack"
> > climatic_TP
> class       : RasterStack
> dimensions  : 277, 412, 114124, 20  (nrow, ncol, ncell, nlayers)
> resolution  : 100, 100  (x, y)
> extent      : 461214, 502414, 4272892, 4300592  (xmin, xmax, ymin, ymax)
> coord. ref. : NA
> names       : bio1_TP, bio10_TP, bio11_TP, bio12_TP, bio13_TP, bio14_TP,
> bio15_TP, bio16_TP, bio17_TP, bio18_TP, bio19_TP, bio2_TP, bio3_TP,
> bio4_TP, bio4a_TP, ...
>
>
> 4th stack
>
> > grids_climatic_FP <- list.files(pattern='asc', full.names=T)# All "asc"
> files in folder
> > climatic_FP<- stack(grids_climatic_FP)# Make a "stack" of EGV #Install
> "stack"
> > climatic_FP
> class       : RasterStack
> dimensions  : 301, 301, 90601, 20  (nrow, ncol, ncell, nlayers)
> resolution  : 100, 100  (x, y)
> extent      : 336088, 366188, 4256440, 4286540  (xmin, xmax, ymin, ymax)
> coord. ref. : NA
> names       : bio1_FP, bio10_FP, bio11_FP, bio12_FP, bio13_FP, bio14_FP,
> bio15_FP, bio16_FP, bio17_FP, bio18_FP, bio19_FP, bio2_FP, bio3_FP,
> bio4_FP, bio4a_FP, ...
>
>
> 5th stack
>
> > grids_climatic_SJP <- list.files(pattern='asc', full.names=T)# All "asc"
> files in folder
> > climatic_SJP<- stack(grids_climatic_SJP)# Make a "stack" of EGV #Install
> "stack"
> > climatic_SJP
> class       : RasterStack
> dimensions  : 602, 602, 362404, 20  (nrow, ncol, ncell, nlayers)
> resolution  : 100, 100  (x, y)
> extent      : 381100, 441300, 4249487, 4309687  (xmin, xmax, ymin, ymax)
> coord. ref. : NA
> names       : bio1_SJP, bio10_SJP, bio11_SJP, bio12_SJP, bio13_SJP,
> bio14_SJP, bio15_SJP, bio16_SJP, bio17_SJP, bio18_SJP, bio19_SJP, bio2_SJP,
> bio3_SJP, bio4_SJP, bio4a_SJP, ...
>
> >
>
>         [[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

Is there a way to merge different stacks?

Wed, 02/13/2019 - 08:47
Hello everybody,
I am having a problem to merge several stacks. My stacks have different
extensions ... It is a problem.

If possible, I want to merge several (5 stacks) asc grids into in a "global
grid/stack"
How can I do it? "Merging", "joining" ...?
Any sugestion?

Thanks,

Lara Silva

My outputs:

1st stack

> grids_climatic_SMP <- list.files(pattern='asc', full.names=T)# All "asc"
files in folder
> climatic_SMP<- stack(grids_climatic_SMP)# Make a "stack" of EGV #Install
"stack"
> climatic_SMP
class       : RasterStack
dimensions  : 901, 901, 811801, 20  (nrow, ncol, ncell, nlayers)
resolution  : 100, 100  (x, y)
extent      : 586974, 677074, 4140847, 4230947  (xmin, xmax, ymin, ymax)
coord. ref. : NA
names       : bio1_SMP, bio10_SMP, bio11_SMP, bio12_SMP, bio13_SMP,
bio14_SMP, bio15_SMP, bio16_SMP, bio17_SMP, bio18_SMP, bio19_SMP, bio2_SMP,
bio3_SMP, bio4_SMP, bio4a_SMP, ...

2nd stack

> grids_climatic_PP <- list.files(pattern='asc', full.names=T)# All "asc"
files in folder
> climatic_PP<- stack(grids_climatic_PP)# Make a "stack" of EGV #Install
"stack"
> climatic_PP
class       : RasterStack
dimensions  : 601, 601, 361201, 20  (nrow, ncol, ncell, nlayers)
resolution  : 100, 100  (x, y)
extent      : 358603, 418703, 4229376, 4289476  (xmin, xmax, ymin, ymax)
coord. ref. : NA
names       : bio1_PP, bio10_PP, bio11_PP, bio12_PP, bio13_PP, bio14_PP,
bio15_PP, bio16_PP, bio17_PP, bio18_PP, bio19_PP, bio2_PP, bio3_PP,
bio4_PP, bio4a_PP, ...

3rd stack

> grids_climatic_TP <- list.files(pattern='asc', full.names=T)# All "asc"
files in folder
> climatic_TP<- stack(grids_climatic_TP)# Make a "stack" of EGV #Install
"stack"
> climatic_TP
class       : RasterStack
dimensions  : 277, 412, 114124, 20  (nrow, ncol, ncell, nlayers)
resolution  : 100, 100  (x, y)
extent      : 461214, 502414, 4272892, 4300592  (xmin, xmax, ymin, ymax)
coord. ref. : NA
names       : bio1_TP, bio10_TP, bio11_TP, bio12_TP, bio13_TP, bio14_TP,
bio15_TP, bio16_TP, bio17_TP, bio18_TP, bio19_TP, bio2_TP, bio3_TP,
bio4_TP, bio4a_TP, ...


4th stack

> grids_climatic_FP <- list.files(pattern='asc', full.names=T)# All "asc"
files in folder
> climatic_FP<- stack(grids_climatic_FP)# Make a "stack" of EGV #Install
"stack"
> climatic_FP
class       : RasterStack
dimensions  : 301, 301, 90601, 20  (nrow, ncol, ncell, nlayers)
resolution  : 100, 100  (x, y)
extent      : 336088, 366188, 4256440, 4286540  (xmin, xmax, ymin, ymax)
coord. ref. : NA
names       : bio1_FP, bio10_FP, bio11_FP, bio12_FP, bio13_FP, bio14_FP,
bio15_FP, bio16_FP, bio17_FP, bio18_FP, bio19_FP, bio2_FP, bio3_FP,
bio4_FP, bio4a_FP, ...


5th stack

> grids_climatic_SJP <- list.files(pattern='asc', full.names=T)# All "asc"
files in folder
> climatic_SJP<- stack(grids_climatic_SJP)# Make a "stack" of EGV #Install
"stack"
> climatic_SJP
class       : RasterStack
dimensions  : 602, 602, 362404, 20  (nrow, ncol, ncell, nlayers)
resolution  : 100, 100  (x, y)
extent      : 381100, 441300, 4249487, 4309687  (xmin, xmax, ymin, ymax)
coord. ref. : NA
names       : bio1_SJP, bio10_SJP, bio11_SJP, bio12_SJP, bio13_SJP,
bio14_SJP, bio15_SJP, bio16_SJP, bio17_SJP, bio18_SJP, bio19_SJP, bio2_SJP,
bio3_SJP, bio4_SJP, bio4a_SJP, ...

>

        [[alternative HTML version deleted]]

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

Sirad for hourly extraterrestrial solar radiation

Wed, 02/13/2019 - 05:06
Hello everybody,

I use SIrad for estimating daily extraterrestrial solar radiation for
Sardinia using dem, in this way

 

DAY    = "2019-02-13"

dem$Ra <- extrat(dayOfYear(DAY),
radians(dem$lat_250))$ExtraTerrestrialSolarRadiationDaily

 

now I need to map hourly values for specific hours, but when I calculate
$ExtraTerrestrialSolarRadiationHorly

I obtain the whole set of 24 values for each grid point, ma question is ..
how can I extract and plot only specific hourly values, for example maximum
values (12 o' clock)?  

Thank you

Regards, Michele

 

 

Michele Fiori

ARPAS - Agenzia Regionale per la Protezione dell'Ambiente della Sardegna
Dipartimento Meteoclimatico - Servizio Meteorologico Agrometeorologico ed
Ecosistemi

Viale Porto Torres 119 - 07100 Sassari
Tel + 39 079 258617 Fax +39 262681

 

 <http://www.sardegnaambiente.it/arpas> www.sardegnaambiente.it/arpas

 


        [[alternative HTML version deleted]]

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

Re: Sampling random directional lines within a polygon

Tue, 02/12/2019 - 23:48
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 endpointtrj<-TrajFromCoords(point_trj,xCol=1,yCol=2,timeCol=3, spatialUnits="m", timeUnits="d")
# generate a random trajectory with 9 segmentstrack<-TrajGenerate(n=9, random=T, stepLength=TrajLength(trj)/9)
# To be able to use StartEndTranslate transform trajectories into matricesma1<-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 wellp_df<-as.data.frame(p)
# change lat / longs to nummericp_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 itxy<-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 <[hidden email]> 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 <[hidden email]> 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 isn't too difficult, but anything that you might want to put into a statistical test (eg testing some hypothesis about the birds maximum deviation from the straight line A-B) needs a lot more care in formulating the path generating process.
Here's some thoughts - if you consider one of the red segments in your map as a piece of string (rather than 10 segments) anchored at the points, then you can stretch it taut with a pencil and draw an ellipse with A and B as the foci. Any path created with that string - taut or slack as in your map - has to strictly lie within the ellipse. Now you can wiggle that string inside that ellipse and create an infinity of paths from A to B of the same length. I'm not sure how you can sample uniformly from that infinity such that any path has an equal sampling probability. Your problem is similar but has the additional rigid segment constraint.

Any two adjacent segments of a chain, eg 1----2-------3, as long as it isn't taut (ie straight) can be perturbed by holding 1 and 3 still and moving 2 to the "mirror image" point over the straight line from 1 to 3. You can also take three segments 1--2--3--4 and hold 1 and 4 still and perturb 2 and 3 fairly easily. In this way you could set up an initial chain and then run multiple perturbations on the chain to get a "random" chain, but quite what set of all chains it would be a sample from is not clear. It could at least generate reasonable looking paths, but I wouldn't want to test a hypothesis against it.
I'm going to generate a path from my office to my home now.
Barry





On Wed, Feb 6, 2019 at 7:50 PM Hannah Justen <[hidden email]> wrote:
Hello everyone,
Thank you very much for the suggestions.

Regarding more details (please also see attached figure):  I would like to simulate a bird's migration between breeding (starting polygon - blue in the figure) and wintering grounds (end polygon - green in the figure). The lines can start from anywhere within the starting polygon and end anywhere in the end polygon. The lines shall consist of 10 connected segments with variable length between 300 and 1000 km (see red lines in figure; two examples of possible lines with 10 segments between the polygons). 
Thank you very much for you help, Hannah
--- PhD Student |Ecology and Evolutionary Biology Texas A&M University
On Wed, Feb 6, 2019 at 5:12 AM Barry Rowlingson <[hidden email]> wrote:
Interesting, but I think we need more details...
Do the lines have to start and finish at specific locations in the polygons - like the centroid, or anywhere?

So one line might be 3 segments of 10km each connecting two polygon centroids that are 15km apart? Imagining three rigid rods of length 10 connected at their ends and with the first and last also connected to two fixed points tells me there's an infinite number of possible solutions. There's probably also a number of ways of sampling from those solutions. Its going to get very complicated with a larger number of segments.

Hmmmmm.....
Barry



On Tue, Feb 5, 2019 at 11:30 PM Hannah Justen <[hidden email]> wrote:
Hi everyone,

I am studying migratory tracks of birds for my dissertation and I would
like to model possible pathways between two polygons. Therefore, I would
like to sample random lines between the polygons. These lines can differ in
total length but should consist of x - number of fragments of equal length.
Each fragment can have slightly different orientation but overall the lines
should connect the two polygons.

I fail to find the appropriate R package that will allow me to do this type
of analysis. Does anyone have a suggestion how to approach analysis?

Thank you,
Hannah

---
PhD Student |Ecology and Evolutionary Biology
Texas A&M University

        [[alternative HTML version deleted]]

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

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

restrict_trajectory_to_extent_of_polygon2.jpg (19K) Download Attachment

Re: converting raster to netcdf .nc with keeping CRS information

Tue, 02/12/2019 - 23:11
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      : -3950000, 3950000, -3950000, 4350000  (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 <[hidden email]> 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      : -3950000, 3950000, -3950000, 4350000  (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      : -3950000, 3950000, -3950000, 4350000  (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-8        LC_COLLATE=en_US.UTF-8
>  [5] LC_MONETARY=en_US.UTF-8    LC_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 <[hidden email]> 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 <[hidden email]> 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 <[hidden email]>
> 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(-3950000, 3950000, -3950000, 4350000)
> >> 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      :
> -3950000, 3950000, -3950000, 4350000  (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      : -3950000, 3950000,
> -3950000, 4350000  (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      : -3950000, 3950000, -3950000, 4350000  (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
> >> <
> https://stackoverflow.com/questions/54593552/saving-r-raster-to-netcdf-nc-file-with-keeping-crs-information
> >,
> >> without a solution so far.
> >>
> >> Thanks
> >>
> >> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> >> *Dr. Ahmed El-Gabbas,*
> >> Ocean Acoustics Lab,
> >> Alfred-Wegener-Institut
> >> Helmholtz-Zentrum für Polar und Meeresforschung
> >> <image.png>
> >> My Website:  https://elgabbas.github.io
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> [hidden email]
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>
> >>
> >> Ben Tupper
> >> Bigelow Laboratory for Ocean Sciences
> >> 60 Bigelow Drive, P.O. Box 380
> >> East Boothbay, Maine 04544
> >> http://www.bigelow.org
> >>
> >> Ecological Forecasting: https://eco.bigelow.org/
> >>
> >>
> >>
> >>
> >>
> >>
> >
> >       [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> Ben Tupper
> Bigelow Laboratory for Ocean Sciences
> 60 Bigelow Drive, P.O. Box 380
> East Boothbay, Maine 04544
> http://www.bigelow.org
>
> Ecological Forecasting: https://eco.bigelow.org/
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

        [[alternative HTML version deleted]]

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

Re: converting raster to netcdf .nc with keeping CRS information

Tue, 02/12/2019 - 19:46
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      : -3950000, 3950000, -3950000, 4350000  (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      : -3950000, 3950000, -3950000, 4350000  (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-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_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 <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> 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(-3950000, 3950000, -3950000, 4350000)
>> 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      : -3950000, 3950000, -3950000, 4350000  (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      : -3950000, 3950000, -3950000, 4350000  (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      : -3950000, 3950000, -3950000, 4350000  (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
>> <https://stackoverflow.com/questions/54593552/saving-r-raster-to-netcdf-nc-file-with-keeping-crs-information>,
>> without a solution so far.
>>
>> Thanks
>>
>> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>> *Dr. Ahmed El-Gabbas,*
>> Ocean Acoustics Lab,
>> Alfred-Wegener-Institut
>> Helmholtz-Zentrum für Polar und Meeresforschung
>> <image.png>
>> My Website:  https://elgabbas.github.io
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>> Ben Tupper
>> Bigelow Laboratory for Ocean Sciences
>> 60 Bigelow Drive, P.O. Box 380
>> East Boothbay, Maine 04544
>> http://www.bigelow.org
>>
>> Ecological Forecasting: https://eco.bigelow.org/
>>
>>
>>
>>
>>
>>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org

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

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

Re: converting raster to netcdf .nc with keeping CRS information

Tue, 02/12/2019 - 08:19
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 <[hidden email]> 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 <[hidden email]> 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(-3950000, 3950000, -3950000, 4350000)
> 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      : -3950000, 3950000, -3950000, 4350000  (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      : -3950000, 3950000, -3950000, 4350000  (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      : -3950000, 3950000, -3950000, 4350000  (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
> <https://stackoverflow.com/questions/54593552/saving-r-raster-to-netcdf-nc-file-with-keeping-crs-information>,
> without a solution so far.
>
> Thanks
>
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> *Dr. Ahmed El-Gabbas,*
> Ocean Acoustics Lab,
> Alfred-Wegener-Institut
> Helmholtz-Zentrum für Polar und Meeresforschung
> <image.png>
> My Website:  https://elgabbas.github.io
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
> Ben Tupper
> Bigelow Laboratory for Ocean Sciences
> 60 Bigelow Drive, P.O. Box 380
> East Boothbay, Maine 04544
> http://www.bigelow.org
>
> Ecological Forecasting: https://eco.bigelow.org/
>
>
>
>
>
>
        [[alternative HTML version deleted]]

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

Pages