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: clip a raster layer to a mask layer using raster algebra in R?

Mon, 12/10/2018 - 08:46
Thank you.

But the most expensive step is constructing the mask with NA values.

 >gebco_rsl[gebco < rsl] <- NA

The rest of the computations is quite fast indeed... Do you know a
faster alternative, for creating the mask, to the line above?

Best,
Christian

Am 09.12.18 um 09:59 schrieb obrl soil:
> Hello,
>
> in R you're working without the benefit of the GRASS 'region' concept
> for raster work (and all the hidden helper code that goes with it), so
> its never going to be a one liner - but the following may help.
>
> You can use raster algebra in R to do masking. Construct your mask so
> that the cells you want to keep are 0 and the cells you want to lose
> are NA.
>
> Then you can just add the raster and its mask together:
>
> masked_raster <- raster + mask
>
> because data + 0 = data, and data + NA = NA. This will still work even
> if the two rasters don't have the same extent, so long as they have
> the same cell size and alignment. The operation will throw a warning
> but will complete successfully, returning the intersection of the two
> objects. This should complete fairly quickly if the mask extent is
> small relative to the target dataset. Try using raster::trim() on the
> mask beforehand if it has a lot of whitespace.
>
> HTH
> @obrl_soil
>
> On Thu, Dec 6, 2018 at 6:38 PM Christian Willmes <[hidden email]> wrote:
>> Hello,
>>
>> this should be normally a very basic simple one-liner, but I can't get
>> my head around it. I only find solutions for cliping a raster against a
>> vector, but none for clipping a raster against a raster mask.
>>
>> I want to clip a layer to a mask. In GRASS GIS r.mapcalc the syntax is:
>>
>> result = if(RSLmask, gebco, null())
>>
>> This clips the GEBCO dataset to the extent of the given raster mask
>> 'rslmask'. The result raster has heigt values within the raster mask
>> extent, the rest of the raster layer has NULL or NA values.
>>
>> In R I did the following, but it does not work as expected:
>>
>> gebco <- raster(gebco2014) #load raster
>>
>> RSLmask <- gebco >= 0  #create mask layer
>>
>> RSL <- ifelse(RSLmask, 1, NA)
>>
>> The last command gives an error about vector type 'S4'?!
>>
>>   >In is.na(test) :
>>     >is.na() auf nicht-(Liste oder Vektor) des Typs 'S4' angewendet
>>
>> So, in short: what is the R way to simply clip a raster to a raster mask?
>>
>> Thanks and best regards,
>> Christian
>>
>> PS: If found many solutions for clipping a raster to a vector polygon,
>> but this is not what I want.
>>
>> --
>> Dr. Christian Willmes
>> AG GIS & Fernerkundung      | GIS & RS Group
>> Geographisches Institut     | Institute of Geography
>> Universität zu Köln         | University of Cologne
>> Tel.: +49 (0)221 470 6234
>> http://www.geographie.uni-koeln.de/14126.html
>> http://www.sfb806.de
>> http://crc806db.uni-koeln.de
>> http://publons.com/a/1316706/
>> http://orcid.org/0000-0002-5566-6542
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Dr. Christian Willmes
AG GIS & Fernerkundung      | GIS & RS Group
Geographisches Institut     | Institute of Geography
Universität zu Köln         | University of Cologne
Tel.: +49 (0)221 470 6234
http://www.geographie.uni-koeln.de/14126.html
http://www.sfb806.de
http://crc806db.uni-koeln.de
http://publons.com/a/1316706/
http://orcid.org/0000-0002-5566-6542

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

Re: clip a raster layer to a mask layer using raster algebra in R?

Sun, 12/09/2018 - 02:59
Hello,

in R you're working without the benefit of the GRASS 'region' concept
for raster work (and all the hidden helper code that goes with it), so
its never going to be a one liner - but the following may help.

You can use raster algebra in R to do masking. Construct your mask so
that the cells you want to keep are 0 and the cells you want to lose
are NA.

Then you can just add the raster and its mask together:

masked_raster <- raster + mask

because data + 0 = data, and data + NA = NA. This will still work even
if the two rasters don't have the same extent, so long as they have
the same cell size and alignment. The operation will throw a warning
but will complete successfully, returning the intersection of the two
objects. This should complete fairly quickly if the mask extent is
small relative to the target dataset. Try using raster::trim() on the
mask beforehand if it has a lot of whitespace.

HTH
@obrl_soil

On Thu, Dec 6, 2018 at 6:38 PM Christian Willmes <[hidden email]> wrote:
>
> Hello,
>
> this should be normally a very basic simple one-liner, but I can't get
> my head around it. I only find solutions for cliping a raster against a
> vector, but none for clipping a raster against a raster mask.
>
> I want to clip a layer to a mask. In GRASS GIS r.mapcalc the syntax is:
>
> result = if(RSLmask, gebco, null())
>
> This clips the GEBCO dataset to the extent of the given raster mask
> 'rslmask'. The result raster has heigt values within the raster mask
> extent, the rest of the raster layer has NULL or NA values.
>
> In R I did the following, but it does not work as expected:
>
> gebco <- raster(gebco2014) #load raster
>
> RSLmask <- gebco >= 0  #create mask layer
>
> RSL <- ifelse(RSLmask, 1, NA)
>
> The last command gives an error about vector type 'S4'?!
>
>  >In is.na(test) :
>    >is.na() auf nicht-(Liste oder Vektor) des Typs 'S4' angewendet
>
> So, in short: what is the R way to simply clip a raster to a raster mask?
>
> Thanks and best regards,
> Christian
>
> PS: If found many solutions for clipping a raster to a vector polygon,
> but this is not what I want.
>
> --
> Dr. Christian Willmes
> AG GIS & Fernerkundung      | GIS & RS Group
> Geographisches Institut     | Institute of Geography
> Universität zu Köln         | University of Cologne
> Tel.: +49 (0)221 470 6234
> http://www.geographie.uni-koeln.de/14126.html
> http://www.sfb806.de
> http://crc806db.uni-koeln.de
> http://publons.com/a/1316706/
> http://orcid.org/0000-0002-5566-6542
>
> _______________________________________________
> 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

Re: clip a raster layer to a mask layer using raster algebra in R?

Thu, 12/06/2018 - 07:22
Possibly there are alternatives, like exploring raster::overlay together
with raster:clusterR
Hugo

Christian Willmes <[hidden email]> escreveu no dia quinta,
6/12/2018 à(s) 13:34:

> Thank you.
>
> It seems that both approaches need an additional, very expensive
> computation (takes several minutes on a i7-4790 with SSD) to create the
> mask or the extent, and in my case the mask computed by this command
> below is already what I want. Only values above a given RSL (Sealevel)
> height. The rest of the raster needs to be NA or NULL.
>
>  >gebco_rsl <- gebco
>  >gebco_rsl[gebco < rsl] <- NA
>
> ....
>
> The easy (and faster) approach for creating a mask:
>
>  >gebco_mask <- gebco >= rsl
>
> writes 0's into the remaining cells, which does not work for me (and the
> mask() function).
>
> On the other hand, the GRASS mapcalc computation (given below) takes
> less than a minute on that same raster.
>
> I wanted to re-implement the whole analysis in R for better/easier
> reproducibility, but now I think I will stick to R + GRASS GIS.
>
> If anyone knows a faster approach for creating a mask in R, please let
> me know.
>
> Thanks and best regards,
> Christian
>
> Am 06.12.18 um 10:59 schrieb Bede-Fazekas Ákos:
> > Dear Christian,
> > raster::mask() and raster::crop() may solve your problem.
> >
> https://www.rdocumentation.org/packages/raster/versions/2.8-4/topics/mask
> >
> https://www.rdocumentation.org/packages/raster/versions/2.8-4/topics/crop
> >
> > HTH,
> > Ákos Bede-Fazekas
> > Hungarian Academy of Sciences
> >
> > 2018.12.06. 9:40 keltezéssel, Christian Willmes írta:
> >> Sorry, I had an error in the main command.
> >>
> >> It should be:
> >>
> >> gebco_clipped <- ifelse(RSLmask, gebco, NA)
> >>
> >> But it still gives the same error.
> >>
> >> Best,
> >> Christian
> >>
> >> Am 06.12.18 um 09:37 schrieb Christian Willmes:
> >>> Hello,
> >>>
> >>> this should be normally a very basic simple one-liner, but I can't
> >>> get my head around it. I only find solutions for cliping a raster
> >>> against a vector, but none for clipping a raster against a raster mask.
> >>>
> >>> I want to clip a layer to a mask. In GRASS GIS r.mapcalc the syntax is:
> >>>
> >>> result = if(RSLmask, gebco, null())
> >>>
> >>> This clips the GEBCO dataset to the extent of the given raster mask
> >>> 'rslmask'. The result raster has heigt values within the raster mask
> >>> extent, the rest of the raster layer has NULL or NA values.
> >>>
> >>> In R I did the following, but it does not work as expected:
> >>>
> >>> gebco <- raster(gebco2014) #load raster
> >>>
> >>> RSLmask <- gebco >= 0  #create mask layer
> >>>
> >>> RSL <- ifelse(RSLmask, 1, NA)
> >>>
> >>> The last command gives an error about vector type 'S4'?!
> >>>
> >>> >In is.na(test) :
> >>>   >is.na() auf nicht-(Liste oder Vektor) des Typs 'S4' angewendet
> >>>
> >>> So, in short: what is the R way to simply clip a raster to a raster
> >>> mask?
> >>>
> >>> Thanks and best regards,
> >>> Christian
> >>>
> >>> PS: If found many solutions for clipping a raster to a vector
> >>> polygon, but this is not what I want.
> >>>
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> --
> Dr. Christian Willmes
> AG GIS & Fernerkundung      | GIS & RS Group
> Geographisches Institut     | Institute of Geography
> Universität zu Köln         | University of Cologne
> Tel.: +49 (0)221 470 6234
> http://www.geographie.uni-koeln.de/14126.html
> http://www.sfb806.de
> http://crc806db.uni-koeln.de
> http://publons.com/a/1316706/
> http://orcid.org/0000-0002-5566-6542
>
> _______________________________________________
> 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: clip a raster layer to a mask layer using raster algebra in R?

Thu, 12/06/2018 - 06:33
Thank you.

It seems that both approaches need an additional, very expensive
computation (takes several minutes on a i7-4790 with SSD) to create the
mask or the extent, and in my case the mask computed by this command
below is already what I want. Only values above a given RSL (Sealevel)
height. The rest of the raster needs to be NA or NULL.

 >gebco_rsl <- gebco
 >gebco_rsl[gebco < rsl] <- NA

....

The easy (and faster) approach for creating a mask:

 >gebco_mask <- gebco >= rsl

writes 0's into the remaining cells, which does not work for me (and the
mask() function).

On the other hand, the GRASS mapcalc computation (given below) takes
less than a minute on that same raster.

I wanted to re-implement the whole analysis in R for better/easier
reproducibility, but now I think I will stick to R + GRASS GIS.

If anyone knows a faster approach for creating a mask in R, please let
me know.

Thanks and best regards,
Christian

Am 06.12.18 um 10:59 schrieb Bede-Fazekas Ákos:
> Dear Christian,
> raster::mask() and raster::crop() may solve your problem.
> https://www.rdocumentation.org/packages/raster/versions/2.8-4/topics/mask
> https://www.rdocumentation.org/packages/raster/versions/2.8-4/topics/crop
>
> HTH,
> Ákos Bede-Fazekas
> Hungarian Academy of Sciences
>
> 2018.12.06. 9:40 keltezéssel, Christian Willmes írta:
>> Sorry, I had an error in the main command.
>>
>> It should be:
>>
>> gebco_clipped <- ifelse(RSLmask, gebco, NA)
>>
>> But it still gives the same error.
>>
>> Best,
>> Christian
>>
>> Am 06.12.18 um 09:37 schrieb Christian Willmes:
>>> Hello,
>>>
>>> this should be normally a very basic simple one-liner, but I can't
>>> get my head around it. I only find solutions for cliping a raster
>>> against a vector, but none for clipping a raster against a raster mask.
>>>
>>> I want to clip a layer to a mask. In GRASS GIS r.mapcalc the syntax is:
>>>
>>> result = if(RSLmask, gebco, null())
>>>
>>> This clips the GEBCO dataset to the extent of the given raster mask
>>> 'rslmask'. The result raster has heigt values within the raster mask
>>> extent, the rest of the raster layer has NULL or NA values.
>>>
>>> In R I did the following, but it does not work as expected:
>>>
>>> gebco <- raster(gebco2014) #load raster
>>>
>>> RSLmask <- gebco >= 0  #create mask layer
>>>
>>> RSL <- ifelse(RSLmask, 1, NA)
>>>
>>> The last command gives an error about vector type 'S4'?!
>>>
>>> >In is.na(test) :
>>>   >is.na() auf nicht-(Liste oder Vektor) des Typs 'S4' angewendet
>>>
>>> So, in short: what is the R way to simply clip a raster to a raster
>>> mask?
>>>
>>> Thanks and best regards,
>>> Christian
>>>
>>> PS: If found many solutions for clipping a raster to a vector
>>> polygon, but this is not what I want.
>>>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Dr. Christian Willmes
AG GIS & Fernerkundung      | GIS & RS Group
Geographisches Institut     | Institute of Geography
Universität zu Köln         | University of Cologne
Tel.: +49 (0)221 470 6234
http://www.geographie.uni-koeln.de/14126.html
http://www.sfb806.de
http://crc806db.uni-koeln.de
http://publons.com/a/1316706/
http://orcid.org/0000-0002-5566-6542

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

Re: clip a raster layer to a mask layer using raster algebra in R?

Thu, 12/06/2018 - 03:59
Dear Christian,
raster::mask() and raster::crop() may solve your problem.
https://www.rdocumentation.org/packages/raster/versions/2.8-4/topics/mask
https://www.rdocumentation.org/packages/raster/versions/2.8-4/topics/crop

HTH,
Ákos Bede-Fazekas
Hungarian Academy of Sciences

2018.12.06. 9:40 keltezéssel, Christian Willmes írta:
> Sorry, I had an error in the main command.
>
> It should be:
>
> gebco_clipped <- ifelse(RSLmask, gebco, NA)
>
> But it still gives the same error.
>
> Best,
> Christian
>
> Am 06.12.18 um 09:37 schrieb Christian Willmes:
>> Hello,
>>
>> this should be normally a very basic simple one-liner, but I can't
>> get my head around it. I only find solutions for cliping a raster
>> against a vector, but none for clipping a raster against a raster mask.
>>
>> I want to clip a layer to a mask. In GRASS GIS r.mapcalc the syntax is:
>>
>> result = if(RSLmask, gebco, null())
>>
>> This clips the GEBCO dataset to the extent of the given raster mask
>> 'rslmask'. The result raster has heigt values within the raster mask
>> extent, the rest of the raster layer has NULL or NA values.
>>
>> In R I did the following, but it does not work as expected:
>>
>> gebco <- raster(gebco2014) #load raster
>>
>> RSLmask <- gebco >= 0  #create mask layer
>>
>> RSL <- ifelse(RSLmask, 1, NA)
>>
>> The last command gives an error about vector type 'S4'?!
>>
>> >In is.na(test) :
>>   >is.na() auf nicht-(Liste oder Vektor) des Typs 'S4' angewendet
>>
>> So, in short: what is the R way to simply clip a raster to a raster
>> mask?
>>
>> Thanks and best regards,
>> Christian
>>
>> PS: If found many solutions for clipping a raster to a vector
>> polygon, but this is not what I want.
>>
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: clip a raster layer to a mask layer using raster algebra in R?

Thu, 12/06/2018 - 02:40
Sorry, I had an error in the main command.

It should be:

gebco_clipped <- ifelse(RSLmask, gebco, NA)

But it still gives the same error.

Best,
Christian

Am 06.12.18 um 09:37 schrieb Christian Willmes:
> Hello,
>
> this should be normally a very basic simple one-liner, but I can't get
> my head around it. I only find solutions for cliping a raster against
> a vector, but none for clipping a raster against a raster mask.
>
> I want to clip a layer to a mask. In GRASS GIS r.mapcalc the syntax is:
>
> result = if(RSLmask, gebco, null())
>
> This clips the GEBCO dataset to the extent of the given raster mask
> 'rslmask'. The result raster has heigt values within the raster mask
> extent, the rest of the raster layer has NULL or NA values.
>
> In R I did the following, but it does not work as expected:
>
> gebco <- raster(gebco2014) #load raster
>
> RSLmask <- gebco >= 0  #create mask layer
>
> RSL <- ifelse(RSLmask, 1, NA)
>
> The last command gives an error about vector type 'S4'?!
>
> >In is.na(test) :
>   >is.na() auf nicht-(Liste oder Vektor) des Typs 'S4' angewendet
>
> So, in short: what is the R way to simply clip a raster to a raster mask?
>
> Thanks and best regards,
> Christian
>
> PS: If found many solutions for clipping a raster to a vector polygon,
> but this is not what I want.
> --
Dr. Christian Willmes
AG GIS & Fernerkundung      | GIS & RS Group
Geographisches Institut     | Institute of Geography
Universität zu Köln         | University of Cologne
Tel.: +49 (0)221 470 6234
http://www.geographie.uni-koeln.de/14126.html
http://www.sfb806.de
http://crc806db.uni-koeln.de
http://publons.com/a/1316706/
http://orcid.org/0000-0002-5566-6542

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

clip a raster layer to a mask layer using raster algebra in R?

Thu, 12/06/2018 - 02:37
Hello,

this should be normally a very basic simple one-liner, but I can't get
my head around it. I only find solutions for cliping a raster against a
vector, but none for clipping a raster against a raster mask.

I want to clip a layer to a mask. In GRASS GIS r.mapcalc the syntax is:

result = if(RSLmask, gebco, null())

This clips the GEBCO dataset to the extent of the given raster mask
'rslmask'. The result raster has heigt values within the raster mask
extent, the rest of the raster layer has NULL or NA values.

In R I did the following, but it does not work as expected:

gebco <- raster(gebco2014) #load raster

RSLmask <- gebco >= 0  #create mask layer

RSL <- ifelse(RSLmask, 1, NA)

The last command gives an error about vector type 'S4'?!

 >In is.na(test) :
   >is.na() auf nicht-(Liste oder Vektor) des Typs 'S4' angewendet

So, in short: what is the R way to simply clip a raster to a raster mask?

Thanks and best regards,
Christian

PS: If found many solutions for clipping a raster to a vector polygon,
but this is not what I want.

--
Dr. Christian Willmes
AG GIS & Fernerkundung      | GIS & RS Group
Geographisches Institut     | Institute of Geography
Universität zu Köln         | University of Cologne
Tel.: +49 (0)221 470 6234
http://www.geographie.uni-koeln.de/14126.html
http://www.sfb806.de
http://crc806db.uni-koeln.de
http://publons.com/a/1316706/
http://orcid.org/0000-0002-5566-6542

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

Re: How to test for a layer in GRASS GIS mapset using rgrass7

Wed, 12/05/2018 - 02:24

Am 04.12.18 um 20:55 schrieb Roger Bivand:
> You need to wrap the test in try() if it may not succeed (in nc):
>
> oo <- try(execGRASS("g.findfile", element="cell", file="slope",
>       mapset="'.'", intern=TRUE), silent=TRUE)
> ifelse(class(oo) == "try-error", FALSE, TRUE)
>
> and use intern to capture the output if the file exists.

Thank you very much, this was exactly what solved the issue.

Best regards,
Christian



--
Dr. Christian Willmes
AG GIS & Fernerkundung      | GIS & RS Group
Geographisches Institut     | Institute of Geography
Universität zu Köln         | University of Cologne
Tel.: +49 (0)221 470 6234
http://www.geographie.uni-koeln.de/14126.html
http://www.sfb806.de
http://crc806db.uni-koeln.de
http://publons.com/a/1316706/
http://orcid.org/0000-0002-5566-6542

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

Re: How to test for a layer in GRASS GIS mapset using rgrass7

Tue, 12/04/2018 - 15:42

On 12/4/18 5:42 PM, Christian Willmes wrote:
Hello,

I am not sure if this is the correct list to ask this specific question about executing GRASS GIS from R Scripts.

My problem is how to determine if a layer exists in the current GRASS GIS mapset or not.

During my script I use r.mapcalc from GRASS GIS to do some expensive computation. So I want to store and reuse the resulting layer, wich is produced according to a variable value, in case I need this layer for the specific variable value again.

If I just omit the overwrite tag it does not work, becasue GRASS GIS stops the execution on this event.

Error in execGRASS("r.mapcalc", expression = expr) : The command:
r.mapcalc expression="rsl40 = GEBCO_2014_2D_4326 >= -40"
produced an error (1) during execution:
FEHLER: output map <rsl40> exists. To overwrite, use the --overwrite flag
Error in execGRASS("r.mapcalc", expression = expr) : The command:
r.mapcalc expression="land_NA_1_rsl40  =  if( rsl40 , 1 ,null())"
produced an error (1) during execution:
FEHLER: output map <land_NA_1_rsl40> exists. To overwrite, use the
        --overwrite flag

So, I tryed to check if a certain layer already exists in the mapset to test if the computation needs/can be executed or not.

Using the following:

if(execGRASS("g.findfile", element="cell", file=lyrname, mapset='"."')){
    return(FALSE)
  }else{
    return(TRUE)
  }

Here it stops, if the layer does not exist:

Fehler in execGRASS("g.findfile", element = "cell", file = lyrname, mapset = "\".\"") :
  The command:
g.findfile element=cell file=rsl40 mapset="."
produced an error (1) during execution:


I think that the mapset parameter requires the actual mapset name, *not* the current working directory.




Does anyone know a solution to this?

Thank you very much!

Best,
Christian


-- Micha Silver Ben Gurion Univ. Sde Boker, Remote Sensing Lab cell: +972-523-665918
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: How to test for a layer in GRASS GIS mapset using rgrass7

Tue, 12/04/2018 - 13:55
On Tue, 4 Dec 2018, Christian Willmes wrote:

> Hello,
>
> I am not sure if this is the correct list to ask this specific question about
> executing GRASS GIS from R Scripts.
>
> My problem is how to determine if a layer exists in the current GRASS GIS
> mapset or not.
>
> During my script I use r.mapcalc from GRASS GIS to do some expensive
> computation. So I want to store and reuse the resulting layer, wich is
> produced according to a variable value, in case I need this layer for the
> specific variable value again.
>
> If I just omit the overwrite tag it does not work, becasue GRASS GIS stops
> the execution on this event.
>
> Error in execGRASS("r.mapcalc", expression = expr) : The command:
> r.mapcalc expression="rsl40 = GEBCO_2014_2D_4326 >= -40"
> produced an error (1) during execution:
> FEHLER: output map <rsl40> exists. To overwrite, use the --overwrite flag
> Error in execGRASS("r.mapcalc", expression = expr) : The command:
> r.mapcalc expression="land_NA_1_rsl40  =  if( rsl40 , 1 ,null())"
> produced an error (1) during execution:
> FEHLER: output map <land_NA_1_rsl40> exists. To overwrite, use the
>         --overwrite flag
>
> So, I tryed to check if a certain layer already exists in the mapset to test
> if the computation needs/can be executed or not.
>
> Using the following:
>
> if(execGRASS("g.findfile", element="cell", file=lyrname, mapset='"."')){
>     return(FALSE)
>   }else{
>     return(TRUE)
>   } You need to wrap the test in try() if it may not succeed (in nc):

oo <- try(execGRASS("g.findfile", element="cell", file="slope",
       mapset="'.'", intern=TRUE), silent=TRUE)
ifelse(class(oo) == "try-error", FALSE, TRUE)

and use intern to capture the output if the file exists. However, you can,
as Rich says, simply overwrite existing files if you wish.

>
> Here it stops, if the layer does not exist:
>

Note that your R script does not stop, simply g.findfile has returned 1
rather than 0 on exit.

Hope this clarifies,

Roger


> Fehler in execGRASS("g.findfile", element = "cell", file = lyrname, mapset =
> "\".\"") :
>   The command:
> g.findfile element=cell file=rsl40 mapset="."
> produced an error (1) during execution:
>
>
> Does anyone know a solution to this?
>
> Thank you very much!
>
> Best,
> Christian
>
>
> --
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Re: How to test for a layer in GRASS GIS mapset using rgrass7

Tue, 12/04/2018 - 12:25
On Tue, 4 Dec 2018, Christian Willmes wrote:

> During my script I use r.mapcalc from GRASS GIS to do some expensive
> computation. So I want to store and reuse the resulting layer, wich is
> produced according to a variable value, in case I need this layer for the
> specific variable value again.
>
> If I just omit the overwrite tag it does not work, becasue GRASS GIS stops
> the execution on this event.

Christian,

   Yes, grass wants the --o flag to overwrite an existing file. You can learn
what raster maps exist using 'g.list type=rast | grep <mapname>'.

HTH,

Rich

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

How to test for a layer in GRASS GIS mapset using rgrass7

Tue, 12/04/2018 - 09:42
Hello,

I am not sure if this is the correct list to ask this specific question
about executing GRASS GIS from R Scripts.

My problem is how to determine if a layer exists in the current GRASS
GIS mapset or not.

During my script I use r.mapcalc from GRASS GIS to do some expensive
computation. So I want to store and reuse the resulting layer, wich is
produced according to a variable value, in case I need this layer for
the specific variable value again.

If I just omit the overwrite tag it does not work, becasue GRASS GIS
stops the execution on this event.

Error in execGRASS("r.mapcalc", expression = expr) : The command:
r.mapcalc expression="rsl40 = GEBCO_2014_2D_4326 >= -40"
produced an error (1) during execution:
FEHLER: output map <rsl40> exists. To overwrite, use the --overwrite flag
Error in execGRASS("r.mapcalc", expression = expr) : The command:
r.mapcalc expression="land_NA_1_rsl40  =  if( rsl40 , 1 ,null())"
produced an error (1) during execution:
FEHLER: output map <land_NA_1_rsl40> exists. To overwrite, use the
         --overwrite flag

So, I tryed to check if a certain layer already exists in the mapset to
test if the computation needs/can be executed or not.

Using the following:

if(execGRASS("g.findfile", element="cell", file=lyrname, mapset='"."')){
     return(FALSE)
   }else{
     return(TRUE)
   }

Here it stops, if the layer does not exist:

Fehler in execGRASS("g.findfile", element = "cell", file = lyrname,
mapset = "\".\"") :
   The command:
g.findfile element=cell file=rsl40 mapset="."
produced an error (1) during execution:


Does anyone know a solution to this?

Thank you very much!

Best,
Christian


--
Dr. Christian Willmes
AG GIS & Fernerkundung      | GIS & RS Group
Geographisches Institut     | Institute of Geography
Universität zu Köln         | University of Cologne
Tel.: +49 (0)221 470 6234
http://www.geographie.uni-koeln.de/14126.html
http://www.sfb806.de
http://crc806db.uni-koeln.de
http://publons.com/a/1316706/
http://orcid.org/0000-0002-5566-6542

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

JSS Special Issue on Bayesian Statistics

Fri, 11/30/2018 - 13:22
The Journal of Statistical Software (JSS) now reached an impact factor of over 22, and recently ended highest in the category "statistics and probability" on the ISI web of knowledge ranking.

Special issues of JSS are very popular and get a lot of citations.  We now invite papers for a special issue of the Journal of Statistical Software on "Bayesian Statistics". The conditions are the regular ones of JSS submissions [1], but the special issue aims at bringing together articles presenting open-source software developed by the authors focusing on (but not limited to) the following topics:

* Bayesian modelling and inference.
* Software for Bayesian Data Analysis.
* Software for Bayesian Statistics.
* Model selection and assessment.
* Priors in Bayesian analysis.
* Visualization of results.
* Bayesian data analysis of highly structured data.
* High-performance computing for Bayesian data analysis.
* Bayesian methods for the analysis of ‘Big Data’.

Authors who intend to submit an article for this special issue are strongly encouraged to submit paper title, authors, and a preliminary abstract to [hidden email] or [hidden email] before December 31st, 2018.

Guest editors for this special issue are:

 Dr. Michela Cameletti (University of Bergamo, Italy)
 Dr. Virgilio Gomez-Rubio (University of Castilla-La Mancha, Spain)
 Prof. Martyn Plummer (Warwick University, UK)

Deadline for papers is Jun 30th, 2019. Submissions must clearly state that they are to be considered for this Special Issue, by mentioning it in the cover letter or by leaving a comment in the upload form.

[1] http://www.jstatsoft.org/

With best regards,

---
Virgilio Gómez Rubio
Escuela de Ingenieros Industriales
Universidad de Castilla-La Mancha
Avda España s/n
02071 Albacete (Spain)

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

Summer School on spatial and spatiotemporal computing: processing large-scale Earth observation (IfGI, University of Münster, Sept 1–7, 2019)

Fri, 11/30/2018 - 03:36

Registrations for the 2019 Summer School on spatial and spatiotemporal
computing: processing large-scale Earth observation (University of
Münster, Sept 1–7, 2019) are now open. For more info see:

https://opengeohub.org/summer_school_2019

Registration deadline: 15th of February 2019 24:00 CET.

This summer school is limited to 65 participants. In the case of higher
number of applications invitations candidates will be selected based on
a ranking system, which is based on: solidarity, academic output and
contributions to the open source projects. To remove geographical bias,
participants coming from more distant areas have a priority on the
rankings list.

Lecturers:
- Edzer Pebesma: Analyzing large amounts of Earth Observation data with
R and openEO
- Roger Bivand: Not just R-spatial: sustaining open source geospatial
software stacks / Data, data everywhere, nor any drop to analyse
(without making brave assumptions about how the data represent
underlying processes)
- Michael Sumner: Computer graphics data structures for geo-spatial /
Building a data library and R toolkit for domain-specific research group
/ Challenges of working with data in polar regions
- Markus Neteler: Cloud based processing of geo and Earth observation
data / Introduction to GRASS GIS / Advanced data analysis in GRASS GIS
- Veronica Andreo: Analysis of space-time satellite data for disease
ecology applications with GRASS GIS and R stats / Analyzing space-time
satellite data with GRASS GIS for environmental monitoring
- Martijn Tennekes: Creating thematic maps in R (tmap package)
- Tomislav Hengl: Computing with large rasters in R: introduction to
tiling and parallelization / Spatial and Spatiotemporal prediction using
ensemble Machine Learning
- Hanna Meyer: Machine learning strategies for spatio-temporal data
- Anita Graser: Analyzing movement data
- Madlene Nussbaum: Mastering machine learning for spatial prediction -
overview and introduction in methods / model selection and
interpretation, uncertainty
- Meng Lu: Assessment of global air pollution exposure

Dates:
- February 15th 2019 — Registration deadline;
- Mid March 2019 — All invitation letters send to applicants;
- May 15th 2019 — Deadline for settling registration fees (working
programme confirmed);
- July 15th 2019 — Final programme, data sets and exercises published;
- Sun 1st September to Sun 8th September 2019 (arrival Sunday, departure
Sunday; 7 night accommodation) Summer School;

Note: OpenGeoHub Summer school will be held week after the FOSS4G
conference (Bucharest, 26–30 August 2019).

Registration fees:
The registrations fees for this Summer School will be in the range
400–500 EUR. Registration fees cover costs of using facilities, lunch
and coffee breaks, administration costs, local travel costs, and costs
of travel and accommodation for lecturers. Participants from ODA
countries (employed by an organization or company in ODA-listed country;
http://www.oecd.org/dac/stats/daclist.htm) typically receive a
subsidized price for the registration costs. Also, full-time students
(MSc or PhD level) are offered subsidized price, provided that they are,
at the moment of application, not employed at any University or research
institute.

Summer school hosts:
This Summer School is hosted by the Institute for Geoinformatics (ifgi),
Heisenbergstr. 2, 48149 Münster. Local organizing committee:
- prof. dr. Edzer Pebesma: Summer School programme, discussion sessions,
- dr. Christian Knoth: logistics, lecture rooms, accommodation, social
programme,

OpenGeoHub is a not-for-profit research foundation with headquarters in
Wageningen, the Netherlands (Stichting OpenGeoHub, KvK 71844570). The
main goal of the OpenGeoHub is to promote publishing and sharing of Open
Geographical and Geoscientific Data and using and developing of Open
Source Software.

Follow us:
https://twitter.com/opengeohub
https://www.youtube.com/c/OpenGeoHubFoundation
https://business.facebook.com/Opengeohub-foundation-274552396602535/
https://plus.google.com/communities/103560544243938791437

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

Re: Reproject MODIS data using R (results in NAs or no spatial extent)

Thu, 11/29/2018 - 22:09
Thank you Katherine and Mike.   You helped identify the major problem, which is that the metadata is incorrect.  The MODIS GLASS product is a climate modeling dataset and as such it is in cylindrical equidistant projection (the metadata incorrectly label it as the sinusoidal grid).  [I'm writing this just to clarify,  I certainly didn't understand this before your emails earlier today.]  The AVHRR data appear to be in the same projection (cylindrical equidistant).  The next step should just be to convert the MODIS units to lat/long.  However, the AVHRR data does not specify a datum per se (see below).


As a follow up, because I am new to this, when converting the MODIS units to lat/long, what should I use as a datum?  Or, given that both MODIS and AVHRR appear to be in the climate modeling grid format, can we assume their datums already match?


AVHRR Datum info:

GEOGCS[\"Unknown datum based upon the Clarke 1866 ellipsoid\","
 [6] "    DATUM[\"Not specified (based on Clarke 1866 spheroid)\","
 [7] "        SPHEROID[\"Clarke 1866\",6378206.4,294.9786982138982,"
 [8] "            AUTHORITY[\"EPSG\",\"7008\"]]],"
 [9] "    PRIMEM[\"Greenwich\",0],"
[10] "    UNIT[\"degree\",0.0174532925199433]]"
[11] "Origin = (-180.000000000000000,90.000000000000000)"
[12] "Pixel Size = (0.050000000000000,-0.050000000000000)"



*Note ESPG 7008<https://epsg.io/7008-ellipsoid> says the units are in meters.  Other options exist, such as ESPG 4008<https://epsg.io/4008>, unknown datum based on the Clarke 1866 elipsoid exist.


Thanks for the help so far and for any future help that might come my way.


Elizabeth

________________________________
From: Kilpatrick, Katherine A <[hidden email]>
Sent: Thursday, November 29, 2018 3:40 PM
To: Michael Sumner
Cc: Elizabeth Webb; [hidden email]
Subject: Re: [R-sig-Geo] Reproject MODIS data using R (results in NAs or no spatial extent)

FYI The link that Mike provided is for ocean color products the GLASS land products use a different grid�see this link.

https://modis-land.gsfc.nasa.gov/MODLAND_grid.html<https://urldefense.proofpoint.com/v2/url?u=https-3A__modis-2Dland.gsfc.nasa.gov_MODLAND-5Fgrid.html&d=DwMGaQ&c=pZJPUDQ3SB9JplYbifm4nt2lEVG5pWx2KikqINpWlZM&r=KpmI6Vdu2He-hM565ejK5Q&m=6y9mIOW6mOGqZ8k_68J89wbwswkZx5iN0QPOMcOQCIw&s=nZLE2iJ2s4P03oznAf9jhIYLo8_O4nbdqj3YRCf4irs&e=>

K




On Nov 29, 2018, at 3:16 PM, Michael Sumner <[hidden email]<mailto:[hidden email]>> wrote:

To answer the "Any ideas on why reprojecting this MODIS data is so
difficult?" - it's that the sinusoidal projection is actually pretty
exotic, it matches the daily aggregation of L2 geophysical variables from
individual swaths (multiple per day) into L3 bins (aggregated daily, then
into longer time periods with bin identity maintained). I personally think
it's a bad choice for a regular gridded product, but it's a way of
"regularizing" the L3 bins, which are quite abstract and hard to use. The
L3 bins are a ragged array of bins at regular latitude intervals, a
diminishing number of bins as absolute latitude increases towards each
pole.  Described here:

https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Foceancolor.gsfc.nasa.gov%2Fdocs%2Fformat%2Fl3bins%2F&amp;data=02%7C01%7Ckkilpatrick%40rsmas.miami.edu%7C08935d04ff104c1b75fd08d65637a4f5%7C2a144b72f23942d48c0e6f0f17c48e33%7C0%7C0%7C636791194328942432&amp;sdata=r4qRfvtrxz2uNp0QQtKoz0X3eAkOcPTrUqg4fUHv5GI%3D&amp;reserved=0<https://urldefense.proofpoint.com/v2/url?u=https-3A__na01.safelinks.protection.outlook.com_-3Furl-3Dhttps-253A-252F-252Foceancolor.gsfc.nasa.gov-252Fdocs-252Fformat-252Fl3bins-252F-26amp-3Bdata-3D02-257C01-257Ckkilpatrick-2540rsmas.miami.edu-257C08935d04ff104c1b75fd08d65637a4f5-257C2a144b72f23942d48c0e6f0f17c48e33-257C0-257C0-257C636791194328942432-26amp-3Bsdata-3Dr4qRfvtrxz2uNp0QQtKoz0X3eAkOcPTrUqg4fUHv5GI-253D-26amp-3Breserved-3D0&d=DwMGaQ&c=pZJPUDQ3SB9JplYbifm4nt2lEVG5pWx2KikqINpWlZM&r=KpmI6Vdu2He-hM565ejK5Q&m=6y9mIOW6mOGqZ8k_68J89wbwswkZx5iN0QPOMcOQCIw&s=j5HW0q-bCYjevZppbmv_FnjcMCcz3MDRIA3O7k3uGEM&e=>

My answer to the "how to reproject the sinusoidal grid to regular longlat"
is - you should not do that. The sinusoidal grid is faithful to the
statistical properties of the L3 bins, though not to the bin geometries
themselves. One is meant to query the sinusoidal grid for the data needed,
it's really inappropriate to stretch those data out by resampling/
reprojection.  (I don't know if the sinusoidal-regular-grid creators have
written down this rationale).

Generally one can reproject vector points, polygons, lines into the
sinusoidal projection and use extractions in that coordinate system, and
that's what I'd recommend. (Using the L3 bins is also possible in R but
it's a long story about how to do that).

Just a question: do you download the auxialiary .hdf.xml files that
accompany the .hdf files?   They perhaps contain extra information that
GDAL can interpret.

Using some inside knowledge, delve into subdatasets with stars/sf:

(sds <- sf::gdal_subdatasets("GLASS02B06.V04.A2013169.2017128.hdf"))

x <- stars::read_stars(sds[[1]])

x
stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
HDF4_EOS:EOS_GRID:"GLASS02B06.V04.A2013169.2017128.hdf":GLASS02B06:ABD_BSA_VIS
Min.   : NA

1st Qu.: NA

Median : NA

Mean   :NaN

3rd Qu.: NA

Max.   : NA

NA's   :1e+05

dimension(s):
 from   to    offset    delta                       refsys point values
x    1 7200 -20015109  154.438 +proj=sinu +lon_0=0 +x_0=...    NA   NULL [x]
y    1 3600  10007555 -308.875 +proj=sinu +lon_0=0 +x_0=...    NA   NULL [y]


plot(x)
sf::st_bbox(x)
    xmin      ymin      xmax      ymax
-20015109   8895604 -18903159  10007555

And OMG! To my disgust we see these data are not in sinusoidal projection
at all, but in Equidistant Cylindrical (the -20015109 is approximately the
distance from the prime meridian to the anti meridian along the equator
(i.e about pi * 6378137m).   So, I don't trust these in-situ metadata at
all.  I see this stuff gets messed up fairly frequently, perhaps it's a
GDAL problem interpreting the metadata from the file, or perhaps it's
actually mis-applied in the hdf - finding out requires careful examination
of the metadata, as does choice of datum (as follows).

Using stars to read and raster to "fix it" (raster just because I'm more
sure of my steps that way):

r <- setExtent(stars:::st_as_raster(x), extent(-180, 180, -90, 90))
projection(r) <- "+proj=longlat +datum=WGS84"
plot(r)
maps::map(add = T)
plot(extent(r), add = TRUE)

## convince yourself the extent is correct with things like:
abline(h = c(-90, 90), col = "red")

That use of "datum=WGS84" may well be incorrect, but compared to trying to
project data from a completely mis applied projection it's an improvement.

HTH

Cheers, Mike.




On Fri, 30 Nov 2018 at 06:39 Michael Sumner <[hidden email]> wrote:

Ah, never mind - it's the subdataset discovery that's probably not easy
with rgdal.
Sorry for the noise.

Mike.

On Fri, 30 Nov 2018 at 06:38 Michael Sumner <[hidden email]> wrote:

Fwiw there shouldn't be any need to convert from hdf to tif - could you
please try this?

x <- readGDAL( ".../GLASS02B05.V04.A1990161.2018062.hdf")

If that works it at least removes some steps from your process.

Cheers, Mike.

On Fri, Nov 30, 2018, 05:03 Elizabeth Webb <[hidden email]> wrote:

I am using GLASS albedo data stored here<
ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/AVHRR/> for pre-2000 (AVHRR) data
and here<ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/MODIS/0.05D/> for
post-2000 data (MODIS). My end goal is to create a raster stack of each
month that contains white sky albedo data from 1982-2015. The problem I
have run into is that the MODIS and AVHRR data are in different spatial
reference systems and I can't seem to reproject them to be in the same
system.

I convert from hdf to tif using R like this:

fileavhrr <- ".../GLASS02B05.V04.A1990161.2018062.hdf"
filemodis<-".../GLASS02B06.V04.A2013169.2017128.hdf"
gdal_translate(get_subdatasets(filemodis)[10], dst_dataset =
       ".../modis.tif")
gdal_translate(get_subdatasets(fileavhrr)[8], projwin =
c(-180,90,180,50), dst_dataset = ".../avhrr.tif") #ideally I'd only like
data north of 50 degrees

avhrr<- raster(".../avhrr.tif")

#class       : RasterLayer
#dimensions  : 800, 7200, 5760000  (nrow, ncol, ncell)
#resolution  : 0.05, 0.05  (x, y)
#extent      : -180, 180, 50, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
#values      : -32768, 32767  (min, max)

modis<- raster(".../modis.tif")

#class       : RasterLayer
#dimensions  : 3600, 7200, 25920000  (nrow, ncol, ncell)
#resolution  : 154.4376, 308.8751  (x, y)
#extent   : -20015109, -18903159, 8895604, 10007555  (xmin, xmax, ymin,
ymax)
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
   +b=6371007.181 +units=m +no_defs
#values      : -32768, 32767  (min, max)

Here are things I have tried:

1.) Use the MODIS Reprojection Tool<
https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Flpdaac.usgs.gov%2Ftools%2Fmodis_reprojection_tool&amp;data=02%7C01%7Ckkilpatrick%40rsmas.miami.edu%7C08935d04ff104c1b75fd08d65637a4f5%7C2a144b72f23942d48c0e6f0f17c48e33%7C0%7C0%7C636791194328942432&amp;sdata=WnSIC3kBe%2Fc5%2FeqiAlJxVjhm2BPs9h1XC%2FiPFfp4JCg%3D&amp;reserved=0>. For whatever
reason, this tool seems to think the subdatasets of the MODIS .hdf files
are only one tile (the upper left most tile, tile 0,0) and not the global
dataset. My understanding is that the MODIS data are global (not in
tiles?), so I do not know why the MRT is doing this.


2.) Use the raster package in R.

projectedMODIS <- projectRaster(modis,avhrr,method="bilinear")

This returns a raster with values that are all NA:

class       : RasterLayer
dimensions  : 800, 7200, 5760000  (nrow,> ncol, ncell)
resolution  : 0.05, 0.05  (x, y)
extent      : -180, 180,> 50, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
values      : NA, NA  (min, max)

3.) Use the gdalUtils package in R:

gdalwarp(srcfile=get_subdatasets(filemodis)[10], dstfile=
".../gdalMODIS_avhrr.tif", s_srs = crs(modis), t_srs =crs(avhrr) )

This returns a raster with essentially no spatial extent.

gdalMODISavhrr<-raster(".../gdalMODIS_avhrr.tif")
#class       : RasterLayer
#dimensions  : 357, 12850, 4587450  (nrow, ncol, ncell)
#resolution  : 0.02801551, 0.02801573  (x, y)
#extent      : -180, 179.9993, 79.99838, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
#values      : -32768, 32767  (min, max)

Any ideas on why reprojecting this MODIS data is so difficult?


       [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&amp;data=02%7C01%7Ckkilpatrick%40rsmas.miami.edu%7C08935d04ff104c1b75fd08d65637a4f5%7C2a144b72f23942d48c0e6f0f17c48e33%7C0%7C0%7C636791194328942432&amp;sdata=NLV%2FZLET%2FHxlxlsl93Laqra4z5A5bDdleHi198Uxh0I%3D&amp;reserved=0

--
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

--
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

--
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://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&amp;data=02%7C01%7Ckkilpatrick%40rsmas.miami.edu%7C08935d04ff104c1b75fd08d65637a4f5%7C2a144b72f23942d48c0e6f0f17c48e33%7C0%7C0%7C636791194328942432&amp;sdata=NLV%2FZLET%2FHxlxlsl93Laqra4z5A5bDdleHi198Uxh0I%3D&amp;reserved=0


        [[alternative HTML version deleted]]


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

Re: Reproject MODIS data using R (results in NAs or no spatial extent)

Thu, 11/29/2018 - 14:40
FYI The link that Mike provided is for ocean color products the GLASS land products use a different grid…see this link.

https://modis-land.gsfc.nasa.gov/MODLAND_grid.html

K




On Nov 29, 2018, at 3:16 PM, Michael Sumner <[hidden email]<mailto:[hidden email]>> wrote:

To answer the "Any ideas on why reprojecting this MODIS data is so
difficult?" - it's that the sinusoidal projection is actually pretty
exotic, it matches the daily aggregation of L2 geophysical variables from
individual swaths (multiple per day) into L3 bins (aggregated daily, then
into longer time periods with bin identity maintained). I personally think
it's a bad choice for a regular gridded product, but it's a way of
"regularizing" the L3 bins, which are quite abstract and hard to use. The
L3 bins are a ragged array of bins at regular latitude intervals, a
diminishing number of bins as absolute latitude increases towards each
pole.  Described here:

https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Foceancolor.gsfc.nasa.gov%2Fdocs%2Fformat%2Fl3bins%2F&amp;data=02%7C01%7Ckkilpatrick%40rsmas.miami.edu%7C08935d04ff104c1b75fd08d65637a4f5%7C2a144b72f23942d48c0e6f0f17c48e33%7C0%7C0%7C636791194328942432&amp;sdata=r4qRfvtrxz2uNp0QQtKoz0X3eAkOcPTrUqg4fUHv5GI%3D&amp;reserved=0

My answer to the "how to reproject the sinusoidal grid to regular longlat"
is - you should not do that. The sinusoidal grid is faithful to the
statistical properties of the L3 bins, though not to the bin geometries
themselves. One is meant to query the sinusoidal grid for the data needed,
it's really inappropriate to stretch those data out by resampling/
reprojection.  (I don't know if the sinusoidal-regular-grid creators have
written down this rationale).

Generally one can reproject vector points, polygons, lines into the
sinusoidal projection and use extractions in that coordinate system, and
that's what I'd recommend. (Using the L3 bins is also possible in R but
it's a long story about how to do that).

Just a question: do you download the auxialiary .hdf.xml files that
accompany the .hdf files?   They perhaps contain extra information that
GDAL can interpret.

Using some inside knowledge, delve into subdatasets with stars/sf:

(sds <- sf::gdal_subdatasets("GLASS02B06.V04.A2013169.2017128.hdf"))

x <- stars::read_stars(sds[[1]])

x
stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
HDF4_EOS:EOS_GRID:"GLASS02B06.V04.A2013169.2017128.hdf":GLASS02B06:ABD_BSA_VIS
Min.   : NA

1st Qu.: NA

Median : NA

Mean   :NaN

3rd Qu.: NA

Max.   : NA

NA's   :1e+05

dimension(s):
 from   to    offset    delta                       refsys point values
x    1 7200 -20015109  154.438 +proj=sinu +lon_0=0 +x_0=...    NA   NULL [x]
y    1 3600  10007555 -308.875 +proj=sinu +lon_0=0 +x_0=...    NA   NULL [y]


plot(x)
sf::st_bbox(x)
    xmin      ymin      xmax      ymax
-20015109   8895604 -18903159  10007555

And OMG! To my disgust we see these data are not in sinusoidal projection
at all, but in Equidistant Cylindrical (the -20015109 is approximately the
distance from the prime meridian to the anti meridian along the equator
(i.e about pi * 6378137m).   So, I don't trust these in-situ metadata at
all.  I see this stuff gets messed up fairly frequently, perhaps it's a
GDAL problem interpreting the metadata from the file, or perhaps it's
actually mis-applied in the hdf - finding out requires careful examination
of the metadata, as does choice of datum (as follows).

Using stars to read and raster to "fix it" (raster just because I'm more
sure of my steps that way):

r <- setExtent(stars:::st_as_raster(x), extent(-180, 180, -90, 90))
projection(r) <- "+proj=longlat +datum=WGS84"
plot(r)
maps::map(add = T)
plot(extent(r), add = TRUE)

## convince yourself the extent is correct with things like:
abline(h = c(-90, 90), col = "red")

That use of "datum=WGS84" may well be incorrect, but compared to trying to
project data from a completely mis applied projection it's an improvement.

HTH

Cheers, Mike.




On Fri, 30 Nov 2018 at 06:39 Michael Sumner <[hidden email]> wrote:

Ah, never mind - it's the subdataset discovery that's probably not easy
with rgdal.
Sorry for the noise.

Mike.

On Fri, 30 Nov 2018 at 06:38 Michael Sumner <[hidden email]> wrote:

Fwiw there shouldn't be any need to convert from hdf to tif - could you
please try this?

x <- readGDAL( ".../GLASS02B05.V04.A1990161.2018062.hdf")

If that works it at least removes some steps from your process.

Cheers, Mike.

On Fri, Nov 30, 2018, 05:03 Elizabeth Webb <[hidden email]> wrote:

I am using GLASS albedo data stored here<
ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/AVHRR/> for pre-2000 (AVHRR) data
and here<ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/MODIS/0.05D/> for
post-2000 data (MODIS). My end goal is to create a raster stack of each
month that contains white sky albedo data from 1982-2015. The problem I
have run into is that the MODIS and AVHRR data are in different spatial
reference systems and I can't seem to reproject them to be in the same
system.

I convert from hdf to tif using R like this:

fileavhrr <- ".../GLASS02B05.V04.A1990161.2018062.hdf"
filemodis<-".../GLASS02B06.V04.A2013169.2017128.hdf"
gdal_translate(get_subdatasets(filemodis)[10], dst_dataset =
       ".../modis.tif")
gdal_translate(get_subdatasets(fileavhrr)[8], projwin =
c(-180,90,180,50), dst_dataset = ".../avhrr.tif") #ideally I'd only like
data north of 50 degrees

avhrr<- raster(".../avhrr.tif")

#class       : RasterLayer
#dimensions  : 800, 7200, 5760000  (nrow, ncol, ncell)
#resolution  : 0.05, 0.05  (x, y)
#extent      : -180, 180, 50, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
#values      : -32768, 32767  (min, max)

modis<- raster(".../modis.tif")

#class       : RasterLayer
#dimensions  : 3600, 7200, 25920000  (nrow, ncol, ncell)
#resolution  : 154.4376, 308.8751  (x, y)
#extent   : -20015109, -18903159, 8895604, 10007555  (xmin, xmax, ymin,
ymax)
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
   +b=6371007.181 +units=m +no_defs
#values      : -32768, 32767  (min, max)

Here are things I have tried:

1.) Use the MODIS Reprojection Tool<
https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Flpdaac.usgs.gov%2Ftools%2Fmodis_reprojection_tool&amp;data=02%7C01%7Ckkilpatrick%40rsmas.miami.edu%7C08935d04ff104c1b75fd08d65637a4f5%7C2a144b72f23942d48c0e6f0f17c48e33%7C0%7C0%7C636791194328942432&amp;sdata=WnSIC3kBe%2Fc5%2FeqiAlJxVjhm2BPs9h1XC%2FiPFfp4JCg%3D&amp;reserved=0>. For whatever
reason, this tool seems to think the subdatasets of the MODIS .hdf files
are only one tile (the upper left most tile, tile 0,0) and not the global
dataset. My understanding is that the MODIS data are global (not in
tiles?), so I do not know why the MRT is doing this.


2.) Use the raster package in R.

projectedMODIS <- projectRaster(modis,avhrr,method="bilinear")

This returns a raster with values that are all NA:

class       : RasterLayer
dimensions  : 800, 7200, 5760000  (nrow,> ncol, ncell)
resolution  : 0.05, 0.05  (x, y)
extent      : -180, 180,> 50, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
values      : NA, NA  (min, max)

3.) Use the gdalUtils package in R:

gdalwarp(srcfile=get_subdatasets(filemodis)[10], dstfile=
".../gdalMODIS_avhrr.tif", s_srs = crs(modis), t_srs =crs(avhrr) )

This returns a raster with essentially no spatial extent.

gdalMODISavhrr<-raster(".../gdalMODIS_avhrr.tif")
#class       : RasterLayer
#dimensions  : 357, 12850, 4587450  (nrow, ncol, ncell)
#resolution  : 0.02801551, 0.02801573  (x, y)
#extent      : -180, 179.9993, 79.99838, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
#values      : -32768, 32767  (min, max)

Any ideas on why reprojecting this MODIS data is so difficult?


       [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&amp;data=02%7C01%7Ckkilpatrick%40rsmas.miami.edu%7C08935d04ff104c1b75fd08d65637a4f5%7C2a144b72f23942d48c0e6f0f17c48e33%7C0%7C0%7C636791194328942432&amp;sdata=NLV%2FZLET%2FHxlxlsl93Laqra4z5A5bDdleHi198Uxh0I%3D&amp;reserved=0

--
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

--
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

--
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://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&amp;data=02%7C01%7Ckkilpatrick%40rsmas.miami.edu%7C08935d04ff104c1b75fd08d65637a4f5%7C2a144b72f23942d48c0e6f0f17c48e33%7C0%7C0%7C636791194328942432&amp;sdata=NLV%2FZLET%2FHxlxlsl93Laqra4z5A5bDdleHi198Uxh0I%3D&amp;reserved=0


        [[alternative HTML version deleted]]

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

Re: Reproject MODIS data using R (results in NAs or no spatial extent)

Thu, 11/29/2018 - 14:16
To answer the "Any ideas on why reprojecting this MODIS data is so
difficult?" - it's that the sinusoidal projection is actually pretty
exotic, it matches the daily aggregation of L2 geophysical variables from
individual swaths (multiple per day) into L3 bins (aggregated daily, then
into longer time periods with bin identity maintained). I personally think
it's a bad choice for a regular gridded product, but it's a way of
"regularizing" the L3 bins, which are quite abstract and hard to use. The
L3 bins are a ragged array of bins at regular latitude intervals, a
diminishing number of bins as absolute latitude increases towards each
pole.  Described here:

https://oceancolor.gsfc.nasa.gov/docs/format/l3bins/

My answer to the "how to reproject the sinusoidal grid to regular longlat"
is - you should not do that. The sinusoidal grid is faithful to the
statistical properties of the L3 bins, though not to the bin geometries
themselves. One is meant to query the sinusoidal grid for the data needed,
it's really inappropriate to stretch those data out by resampling/
reprojection.  (I don't know if the sinusoidal-regular-grid creators have
written down this rationale).

Generally one can reproject vector points, polygons, lines into the
sinusoidal projection and use extractions in that coordinate system, and
that's what I'd recommend. (Using the L3 bins is also possible in R but
it's a long story about how to do that).

 Just a question: do you download the auxialiary .hdf.xml files that
accompany the .hdf files?   They perhaps contain extra information that
GDAL can interpret.

Using some inside knowledge, delve into subdatasets with stars/sf:

(sds <- sf::gdal_subdatasets("GLASS02B06.V04.A2013169.2017128.hdf"))

x <- stars::read_stars(sds[[1]])

x
stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
 HDF4_EOS:EOS_GRID:"GLASS02B06.V04.A2013169.2017128.hdf":GLASS02B06:ABD_BSA_VIS
 Min.   : NA

 1st Qu.: NA

 Median : NA

 Mean   :NaN

 3rd Qu.: NA

 Max.   : NA

 NA's   :1e+05

dimension(s):
  from   to    offset    delta                       refsys point values
x    1 7200 -20015109  154.438 +proj=sinu +lon_0=0 +x_0=...    NA   NULL [x]
y    1 3600  10007555 -308.875 +proj=sinu +lon_0=0 +x_0=...    NA   NULL [y]


plot(x)
sf::st_bbox(x)
     xmin      ymin      xmax      ymax
-20015109   8895604 -18903159  10007555

And OMG! To my disgust we see these data are not in sinusoidal projection
at all, but in Equidistant Cylindrical (the -20015109 is approximately the
distance from the prime meridian to the anti meridian along the equator
(i.e about pi * 6378137m).   So, I don't trust these in-situ metadata at
all.  I see this stuff gets messed up fairly frequently, perhaps it's a
GDAL problem interpreting the metadata from the file, or perhaps it's
actually mis-applied in the hdf - finding out requires careful examination
of the metadata, as does choice of datum (as follows).

Using stars to read and raster to "fix it" (raster just because I'm more
sure of my steps that way):

r <- setExtent(stars:::st_as_raster(x), extent(-180, 180, -90, 90))
projection(r) <- "+proj=longlat +datum=WGS84"
plot(r)
maps::map(add = T)
plot(extent(r), add = TRUE)

## convince yourself the extent is correct with things like:
abline(h = c(-90, 90), col = "red")

That use of "datum=WGS84" may well be incorrect, but compared to trying to
project data from a completely mis applied projection it's an improvement.

HTH

Cheers, Mike.




On Fri, 30 Nov 2018 at 06:39 Michael Sumner <[hidden email]> wrote:

> Ah, never mind - it's the subdataset discovery that's probably not easy
> with rgdal.
> Sorry for the noise.
>
> Mike.
>
> On Fri, 30 Nov 2018 at 06:38 Michael Sumner <[hidden email]> wrote:
>
>> Fwiw there shouldn't be any need to convert from hdf to tif - could you
>> please try this?
>>
>> x <- readGDAL( ".../GLASS02B05.V04.A1990161.2018062.hdf")
>>
>> If that works it at least removes some steps from your process.
>>
>> Cheers, Mike.
>>
>> On Fri, Nov 30, 2018, 05:03 Elizabeth Webb <[hidden email]> wrote:
>>
>>> I am using GLASS albedo data stored here<
>>> ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/AVHRR/> for pre-2000 (AVHRR) data
>>> and here<ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/MODIS/0.05D/> for
>>> post-2000 data (MODIS). My end goal is to create a raster stack of each
>>> month that contains white sky albedo data from 1982-2015. The problem I
>>> have run into is that the MODIS and AVHRR data are in different spatial
>>> reference systems and I can't seem to reproject them to be in the same
>>> system.
>>>
>>> I convert from hdf to tif using R like this:
>>>
>>> fileavhrr <- ".../GLASS02B05.V04.A1990161.2018062.hdf"
>>> filemodis<-".../GLASS02B06.V04.A2013169.2017128.hdf"
>>> gdal_translate(get_subdatasets(filemodis)[10], dst_dataset =
>>>         ".../modis.tif")
>>> gdal_translate(get_subdatasets(fileavhrr)[8], projwin =
>>> c(-180,90,180,50), dst_dataset = ".../avhrr.tif") #ideally I'd only like
>>> data north of 50 degrees
>>>
>>> avhrr<- raster(".../avhrr.tif")
>>>
>>> #class       : RasterLayer
>>> #dimensions  : 800, 7200, 5760000  (nrow, ncol, ncell)
>>> #resolution  : 0.05, 0.05  (x, y)
>>> #extent      : -180, 180, 50, 90  (xmin, xmax, ymin, ymax)
>>> #coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
>>> #values      : -32768, 32767  (min, max)
>>>
>>> modis<- raster(".../modis.tif")
>>>
>>> #class       : RasterLayer
>>> #dimensions  : 3600, 7200, 25920000  (nrow, ncol, ncell)
>>> #resolution  : 154.4376, 308.8751  (x, y)
>>> #extent   : -20015109, -18903159, 8895604, 10007555  (xmin, xmax, ymin,
>>> ymax)
>>> #coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
>>>     +b=6371007.181 +units=m +no_defs
>>> #values      : -32768, 32767  (min, max)
>>>
>>> Here are things I have tried:
>>>
>>> 1.) Use the MODIS Reprojection Tool<
>>> https://lpdaac.usgs.gov/tools/modis_reprojection_tool>. For whatever
>>> reason, this tool seems to think the subdatasets of the MODIS .hdf files
>>> are only one tile (the upper left most tile, tile 0,0) and not the global
>>> dataset. My understanding is that the MODIS data are global (not in
>>> tiles?), so I do not know why the MRT is doing this.
>>>
>>>
>>> 2.) Use the raster package in R.
>>>
>>> projectedMODIS <- projectRaster(modis,avhrr,method="bilinear")
>>>
>>> This returns a raster with values that are all NA:
>>>
>>> class       : RasterLayer
>>> dimensions  : 800, 7200, 5760000  (nrow,> ncol, ncell)
>>> resolution  : 0.05, 0.05  (x, y)
>>> extent      : -180, 180,> 50, 90  (xmin, xmax, ymin, ymax)
>>> coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
>>> values      : NA, NA  (min, max)
>>>
>>> 3.) Use the gdalUtils package in R:
>>>
>>> gdalwarp(srcfile=get_subdatasets(filemodis)[10], dstfile=
>>> ".../gdalMODIS_avhrr.tif", s_srs = crs(modis), t_srs =crs(avhrr) )
>>>
>>> This returns a raster with essentially no spatial extent.
>>>
>>> gdalMODISavhrr<-raster(".../gdalMODIS_avhrr.tif")
>>> #class       : RasterLayer
>>> #dimensions  : 357, 12850, 4587450  (nrow, ncol, ncell)
>>> #resolution  : 0.02801551, 0.02801573  (x, y)
>>> #extent      : -180, 179.9993, 79.99838, 90  (xmin, xmax, ymin, ymax)
>>> #coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
>>> #values      : -32768, 32767  (min, max)
>>>
>>> Any ideas on why reprojecting this MODIS data is so difficult?
>>>
>>>
>>>         [[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
>>
>> --
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> Kingston Tasmania 7050 Australia
>
> -- 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: Reproject MODIS data using R (results in NAs or no spatial extent)

Thu, 11/29/2018 - 14:08
gdalUtils has a get_subdatasets() which is really helpful here, it's
just a gdalinfo with grep.

Here's what my code looks like for reprojecting, what I found is that
even if the hdf has the proj definition, gdal tools don't seem to use it
right, so I specify the proj string as the s_srs or a_srs.

.modis_warp <- function(output, t_srs="EPSG:4326", ...){
  # Multiple sources say this is the proj definition of the MODIS Sinusoidal
  MODIS_SRS <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+b=6371007.181 +units=m +no_defs"

  CO <- c("COMPRESS=DEFLATE","PREDICTOR=2","ZLEVEL=9")

  catcher <- tryCatch(
    gdalwarp(output[1]
             , output[2]
             , s_srs=MODIS_SRS
             , t_srs=t_srs
             , of="GTiff"
             #, srcnodata
             #, dstnodata
             , "-multi"
             #, paste0("-wo NUM_THREADS=",ncpu)#Multithread computatation
             , CO
             , output)
  )

  return(catcher)
}

Sure you could transform the raster in "memory" but that will probably
right a temp grd file to the system anyways, might as well just convert
it to a multi-band tiff for ease of use.

Thanks,
Alex

On 11/29/18 11:39, Michael Sumner wrote:
> Ah, never mind - it's the subdataset discovery that's probably not easy
> with rgdal.
> Sorry for the noise.
>
> Mike.
>
> On Fri, 30 Nov 2018 at 06:38 Michael Sumner <[hidden email]> wrote:
>
>> Fwiw there shouldn't be any need to convert from hdf to tif - could you
>> please try this?
>>
>> x <- readGDAL( ".../GLASS02B05.V04.A1990161.2018062.hdf")
>>
>> If that works it at least removes some steps from your process.
>>
>> Cheers, Mike.
>>
>> On Fri, Nov 30, 2018, 05:03 Elizabeth Webb <[hidden email]> wrote:
>>
>>> I am using GLASS albedo data stored here<
>>> ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/AVHRR/> for pre-2000 (AVHRR) data
>>> and here<ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/MODIS/0.05D/> for
>>> post-2000 data (MODIS). My end goal is to create a raster stack of each
>>> month that contains white sky albedo data from 1982-2015. The problem I
>>> have run into is that the MODIS and AVHRR data are in different spatial
>>> reference systems and I can't seem to reproject them to be in the same
>>> system.
>>>
>>> I convert from hdf to tif using R like this:
>>>
>>> fileavhrr <- ".../GLASS02B05.V04.A1990161.2018062.hdf"
>>> filemodis<-".../GLASS02B06.V04.A2013169.2017128.hdf"
>>> gdal_translate(get_subdatasets(filemodis)[10], dst_dataset =
>>>         ".../modis.tif")
>>> gdal_translate(get_subdatasets(fileavhrr)[8], projwin =
>>> c(-180,90,180,50), dst_dataset = ".../avhrr.tif") #ideally I'd only like
>>> data north of 50 degrees
>>>
>>> avhrr<- raster(".../avhrr.tif")
>>>
>>> #class       : RasterLayer
>>> #dimensions  : 800, 7200, 5760000  (nrow, ncol, ncell)
>>> #resolution  : 0.05, 0.05  (x, y)
>>> #extent      : -180, 180, 50, 90  (xmin, xmax, ymin, ymax)
>>> #coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
>>> #values      : -32768, 32767  (min, max)
>>>
>>> modis<- raster(".../modis.tif")
>>>
>>> #class       : RasterLayer
>>> #dimensions  : 3600, 7200, 25920000  (nrow, ncol, ncell)
>>> #resolution  : 154.4376, 308.8751  (x, y)
>>> #extent   : -20015109, -18903159, 8895604, 10007555  (xmin, xmax, ymin,
>>> ymax)
>>> #coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
>>>     +b=6371007.181 +units=m +no_defs
>>> #values      : -32768, 32767  (min, max)
>>>
>>> Here are things I have tried:
>>>
>>> 1.) Use the MODIS Reprojection Tool<
>>> https://lpdaac.usgs.gov/tools/modis_reprojection_tool>. For whatever
>>> reason, this tool seems to think the subdatasets of the MODIS .hdf files
>>> are only one tile (the upper left most tile, tile 0,0) and not the global
>>> dataset. My understanding is that the MODIS data are global (not in
>>> tiles?), so I do not know why the MRT is doing this.
>>>
>>>
>>> 2.) Use the raster package in R.
>>>
>>> projectedMODIS <- projectRaster(modis,avhrr,method="bilinear")
>>>
>>> This returns a raster with values that are all NA:
>>>
>>> class       : RasterLayer
>>> dimensions  : 800, 7200, 5760000  (nrow,> ncol, ncell)
>>> resolution  : 0.05, 0.05  (x, y)
>>> extent      : -180, 180,> 50, 90  (xmin, xmax, ymin, ymax)
>>> coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
>>> values      : NA, NA  (min, max)
>>>
>>> 3.) Use the gdalUtils package in R:
>>>
>>> gdalwarp(srcfile=get_subdatasets(filemodis)[10], dstfile=
>>> ".../gdalMODIS_avhrr.tif", s_srs = crs(modis), t_srs =crs(avhrr) )
>>>
>>> This returns a raster with essentially no spatial extent.
>>>
>>> gdalMODISavhrr<-raster(".../gdalMODIS_avhrr.tif")
>>> #class       : RasterLayer
>>> #dimensions  : 357, 12850, 4587450  (nrow, ncol, ncell)
>>> #resolution  : 0.02801551, 0.02801573  (x, y)
>>> #extent      : -180, 179.9993, 79.99838, 90  (xmin, xmax, ymin, ymax)
>>> #coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
>>> #values      : -32768, 32767  (min, max)
>>>
>>> Any ideas on why reprojecting this MODIS data is so difficult?
>>>
>>>
>>>         [[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
>>
>> --
> 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
>
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Reproject MODIS data using R (results in NAs or no spatial extent)

Thu, 11/29/2018 - 13:39
Ah, never mind - it's the subdataset discovery that's probably not easy
with rgdal.
Sorry for the noise.

Mike.

On Fri, 30 Nov 2018 at 06:38 Michael Sumner <[hidden email]> wrote:

> Fwiw there shouldn't be any need to convert from hdf to tif - could you
> please try this?
>
> x <- readGDAL( ".../GLASS02B05.V04.A1990161.2018062.hdf")
>
> If that works it at least removes some steps from your process.
>
> Cheers, Mike.
>
> On Fri, Nov 30, 2018, 05:03 Elizabeth Webb <[hidden email]> wrote:
>
>> I am using GLASS albedo data stored here<
>> ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/AVHRR/> for pre-2000 (AVHRR) data
>> and here<ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/MODIS/0.05D/> for
>> post-2000 data (MODIS). My end goal is to create a raster stack of each
>> month that contains white sky albedo data from 1982-2015. The problem I
>> have run into is that the MODIS and AVHRR data are in different spatial
>> reference systems and I can't seem to reproject them to be in the same
>> system.
>>
>> I convert from hdf to tif using R like this:
>>
>> fileavhrr <- ".../GLASS02B05.V04.A1990161.2018062.hdf"
>> filemodis<-".../GLASS02B06.V04.A2013169.2017128.hdf"
>> gdal_translate(get_subdatasets(filemodis)[10], dst_dataset =
>>         ".../modis.tif")
>> gdal_translate(get_subdatasets(fileavhrr)[8], projwin =
>> c(-180,90,180,50), dst_dataset = ".../avhrr.tif") #ideally I'd only like
>> data north of 50 degrees
>>
>> avhrr<- raster(".../avhrr.tif")
>>
>> #class       : RasterLayer
>> #dimensions  : 800, 7200, 5760000  (nrow, ncol, ncell)
>> #resolution  : 0.05, 0.05  (x, y)
>> #extent      : -180, 180, 50, 90  (xmin, xmax, ymin, ymax)
>> #coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
>> #values      : -32768, 32767  (min, max)
>>
>> modis<- raster(".../modis.tif")
>>
>> #class       : RasterLayer
>> #dimensions  : 3600, 7200, 25920000  (nrow, ncol, ncell)
>> #resolution  : 154.4376, 308.8751  (x, y)
>> #extent   : -20015109, -18903159, 8895604, 10007555  (xmin, xmax, ymin,
>> ymax)
>> #coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
>>     +b=6371007.181 +units=m +no_defs
>> #values      : -32768, 32767  (min, max)
>>
>> Here are things I have tried:
>>
>> 1.) Use the MODIS Reprojection Tool<
>> https://lpdaac.usgs.gov/tools/modis_reprojection_tool>. For whatever
>> reason, this tool seems to think the subdatasets of the MODIS .hdf files
>> are only one tile (the upper left most tile, tile 0,0) and not the global
>> dataset. My understanding is that the MODIS data are global (not in
>> tiles?), so I do not know why the MRT is doing this.
>>
>>
>> 2.) Use the raster package in R.
>>
>> projectedMODIS <- projectRaster(modis,avhrr,method="bilinear")
>>
>> This returns a raster with values that are all NA:
>>
>> class       : RasterLayer
>> dimensions  : 800, 7200, 5760000  (nrow,> ncol, ncell)
>> resolution  : 0.05, 0.05  (x, y)
>> extent      : -180, 180,> 50, 90  (xmin, xmax, ymin, ymax)
>> coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
>> values      : NA, NA  (min, max)
>>
>> 3.) Use the gdalUtils package in R:
>>
>> gdalwarp(srcfile=get_subdatasets(filemodis)[10], dstfile=
>> ".../gdalMODIS_avhrr.tif", s_srs = crs(modis), t_srs =crs(avhrr) )
>>
>> This returns a raster with essentially no spatial extent.
>>
>> gdalMODISavhrr<-raster(".../gdalMODIS_avhrr.tif")
>> #class       : RasterLayer
>> #dimensions  : 357, 12850, 4587450  (nrow, ncol, ncell)
>> #resolution  : 0.02801551, 0.02801573  (x, y)
>> #extent      : -180, 179.9993, 79.99838, 90  (xmin, xmax, ymin, ymax)
>> #coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
>> #values      : -32768, 32767  (min, max)
>>
>> Any ideas on why reprojecting this MODIS data is so difficult?
>>
>>
>>         [[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
>
> -- 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: Reproject MODIS data using R (results in NAs or no spatial extent)

Thu, 11/29/2018 - 13:38
Fwiw there shouldn't be any need to convert from hdf to tif - could you
please try this?

x <- readGDAL( ".../GLASS02B05.V04.A1990161.2018062.hdf")

If that works it at least removes some steps from your process.

Cheers, Mike.

On Fri, Nov 30, 2018, 05:03 Elizabeth Webb <[hidden email]> wrote:

> I am using GLASS albedo data stored here<
> ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/AVHRR/> for pre-2000 (AVHRR) data
> and here<ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/MODIS/0.05D/> for
> post-2000 data (MODIS). My end goal is to create a raster stack of each
> month that contains white sky albedo data from 1982-2015. The problem I
> have run into is that the MODIS and AVHRR data are in different spatial
> reference systems and I can't seem to reproject them to be in the same
> system.
>
> I convert from hdf to tif using R like this:
>
> fileavhrr <- ".../GLASS02B05.V04.A1990161.2018062.hdf"
> filemodis<-".../GLASS02B06.V04.A2013169.2017128.hdf"
> gdal_translate(get_subdatasets(filemodis)[10], dst_dataset =
>         ".../modis.tif")
> gdal_translate(get_subdatasets(fileavhrr)[8], projwin = c(-180,90,180,50),
> dst_dataset = ".../avhrr.tif") #ideally I'd only like data north of 50
> degrees
>
> avhrr<- raster(".../avhrr.tif")
>
> #class       : RasterLayer
> #dimensions  : 800, 7200, 5760000  (nrow, ncol, ncell)
> #resolution  : 0.05, 0.05  (x, y)
> #extent      : -180, 180, 50, 90  (xmin, xmax, ymin, ymax)
> #coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
> #values      : -32768, 32767  (min, max)
>
> modis<- raster(".../modis.tif")
>
> #class       : RasterLayer
> #dimensions  : 3600, 7200, 25920000  (nrow, ncol, ncell)
> #resolution  : 154.4376, 308.8751  (x, y)
> #extent   : -20015109, -18903159, 8895604, 10007555  (xmin, xmax, ymin,
> ymax)
> #coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
>     +b=6371007.181 +units=m +no_defs
> #values      : -32768, 32767  (min, max)
>
> Here are things I have tried:
>
> 1.) Use the MODIS Reprojection Tool<
> https://lpdaac.usgs.gov/tools/modis_reprojection_tool>. For whatever
> reason, this tool seems to think the subdatasets of the MODIS .hdf files
> are only one tile (the upper left most tile, tile 0,0) and not the global
> dataset. My understanding is that the MODIS data are global (not in
> tiles?), so I do not know why the MRT is doing this.
>
>
> 2.) Use the raster package in R.
>
> projectedMODIS <- projectRaster(modis,avhrr,method="bilinear")
>
> This returns a raster with values that are all NA:
>
> class       : RasterLayer
> dimensions  : 800, 7200, 5760000  (nrow,> ncol, ncell)
> resolution  : 0.05, 0.05  (x, y)
> extent      : -180, 180,> 50, 90  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
> values      : NA, NA  (min, max)
>
> 3.) Use the gdalUtils package in R:
>
> gdalwarp(srcfile=get_subdatasets(filemodis)[10], dstfile=
> ".../gdalMODIS_avhrr.tif", s_srs = crs(modis), t_srs =crs(avhrr) )
>
> This returns a raster with essentially no spatial extent.
>
> gdalMODISavhrr<-raster(".../gdalMODIS_avhrr.tif")
> #class       : RasterLayer
> #dimensions  : 357, 12850, 4587450  (nrow, ncol, ncell)
> #resolution  : 0.02801551, 0.02801573  (x, y)
> #extent      : -180, 179.9993, 79.99838, 90  (xmin, xmax, ymin, ymax)
> #coord. ref. : +proj=longlat +ellps=clrk66 +no_defs
> #values      : -32768, 32767  (min, max)
>
> Any ideas on why reprojecting this MODIS data is so difficult?
>
>
>         [[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

Pages