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: 2 hours 49 min ago

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

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

Thu, 11/29/2018 - 11:53
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

How to manipulatechange the space between the numbers of axes (equal interval)? How to add a legend to a map (plot raster)?

Thu, 11/29/2018 - 10:42
Hello,
 I´m trying to plot a map in R. Is it possible to change the values of
x.axis and y.axis and have the same interval (space between the numbers)?
 How to add a legend for the map?

https://www.dropbox.com/s/60tavy6vez33mw3/Plot_map.tiff?dl=0

> class(proj22PA$P_EMmeanByROC_mergedAlgo_mergedRun_mergedData)
[1] "RasterLayer"
attr(,"package")
[1] "raster"
>
> proj22PA$P_EMmeanByROC_mergedAlgo_mergedRun_mergedData
class       : RasterLayer
band        : 2  (of  2  bands)
dimensions  : 277, 412, 114124  (nrow, ncol, ncell)
resolution  : 100, 100  (x, y)
extent      : 461214, 502414, 4272892, 4300592  (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : grid
names       : P_EMmeanByROC_mergedAlgo_mergedRun_mergedData
values      : 42, 924  (min, max)

>

Code:

plot(proj22PA$P_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
axes=FALSE,legend=FALSE,box=FALSE,par(mar=c(3,3,0.5,0.5),mgp=c(2,1,0)))
axis(side = 1, at =c(470000,480000,500000),labels= c("470000", "480000",
"500000"))
axis(side = 2, at =c(4272892,4300592),labels= c("4272892", "4300592"))
box()

Regards,

Silva

        [[alternative HTML version deleted]]

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

Re: How to change axis (plot map)?

Sat, 11/24/2018 - 11:50
Hi Hugo,
thank you for your help and advices.

Regards,

Silva

Hugo Costa <[hidden email]> escreveu no dia sexta, 23/11/2018 à(s)
20:27:

> Probably you're looking for arguments xaxp and yaxp, used like this:
> plot(1:100,xaxp=c(0,100,4), yaxp=c(0,100,2))
>
> Please look for more info in the help. Also, if you intend to make maps,
> possibly you should spend some time getting familiar with packages, such as
> ggplot2 or tmap.
>
> Cheers
> Hugo
>
>
> Lara Silva <[hidden email]> escreveu no dia sexta, 23/11/2018
> à(s) 14:31:
>
>> Hello,
>>
>>  I am trying to change the scale/number of x-axis and y-axis (map plot -
>> Figure 1).
>>
>> https://www.dropbox.com/sh/34nddj3405q3a6p/AACe-_23Wx7bZRETYIOYjsiQa?dl=0
>> For example:
>>
>> x-axis: 465000,480000, 500000 (min. med. max values)
>> The same idea for y-axis ....
>>
>> Code:
>>
>> plot(proj22PA$OB_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
>> main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis)
>>
>>
>> plot(proj22PA$PittOB_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
>> main= "", xlab ="x.coords", ylab="y.coords",
>> cex.axis=caxis,xlim=c(465000,500000,480000), ylim=c(4270000,4300000))
>>
>> Do not change ... The map has shifted. Probably, I need a new function or
>> modified the argumente for link the information.
>>
>>
>> Another hypothesis
>>
>> plot(proj22PA$Pittosporum_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
>> main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis,
>> xaxt="n",yaxt="n")
>>
>> axis(1, seq(450000,480000,500000)) - x.axis
>> axis(2, seq(4270000,425000,4300000)) - y.axis
>>
>> Output - Figure 2. I can not add the "numbers".
>>
>> https://www.dropbox.com/sh/34nddj3405q3a6p/AACe-_23Wx7bZRETYIOYjsiQa?dl=0
>>
>> Is it possible change the scale of legend (200, 400, 600, 800)?
>>
>> Any suggestion for my doubts/problems?
>>
>> Thanks,
>>
>> Silva
>>
>>         [[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: How to sum area of adjacent patches?

Fri, 11/23/2018 - 17:48
Hi Ben,
thank you for your help. Clump is a possibility, however I was looking for
something different. In my real data, 9 rasters have to be merged (1
central raster and the 8 surrounding rasters), and running clump after
merging 9 rasters is slow and memory consuming. And this has to be repeated
thousands of times (hopefully in parallel with the rasters in memory). I
didn't explain these details to make the message short.
I'll continue searching for more alternatives, and all ideas are welcome.
Thanks
Hugo

Ben Tupper <[hidden email]> escreveu no dia sexta, 23/11/2018 à(s)
23:58:

> Hi,
>
> I think you may need to use raster::clump() and table() to label your
> regions.  Check out an old message about clump...
>
> https://stat.ethz.ch/pipermail/r-sig-geo/2017-January/025346.html
>
> It's sort of cumbersome, but I think the following works...
>
>
> library(raster)
> r1<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=0, ymx=1)
> r2<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=1, ymx=2)
> r1[]<-c(rep(24,24), rep(0,76))
> r2[]<-c(rep(0,86), rep(14,14))
> par(mfrow=c(3,2)); plot(r1,main="r1"); plot(r2,main="r2")
>
>
> r3<-merge(r1,r2) ; plot(r3, main="merged rasters")
>
> r3 <- r3
> r4[r4>0]<-38     ; plot(r4, main="desired result")
>
> r5 <- clump(r3) ; plot(r5, main = 'clumped')
> r5[is.na(r5)] <- 0
>
> r5t <- table(getValues(r5))
> for (i in names(r5t)[-1]){
>         id <- as.numeric(i)
>         print(id)
>         r5[r5 == id] <- r5t[i]
> }
>
> plot(r5, main = 'relabeled with area')
>
>
> Cheers,
> Ben
>
> P.S. Years ago I used a language called IDL.  It's histogram function has
> an output keyword called reverse_indices (more here
> https://www.harrisgeospatial.com/docs/HISTOGRAM.html#H_835179117_677188
> and here http://www.idlcoyote.com/tips/histogram_tutorial.html )  - super
> useful!  Your situation just begs for that same solution in R.
>
>
>
>
>
> > On Nov 23, 2018, at 5:13 PM, Hugo Costa <[hidden email]> wrote:
> >
> > Deal all,
> >
> > I'm merging several rasters representing patches of forest, and the pixel
> > values represent the area of the patches. The patches' area should be
> > summed along the edges of the merged rasters. I provide a reproducible
> > example. Any help is welcome.
> > Hugo
> >
> > library(raster)
> > r1<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=0, ymx=1)
> > r2<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=1, ymx=2)
> > r1[]<-c(rep(24,24), rep(0,76))
> > r2[]<-c(rep(0,86), rep(14,14))
> > par(mfrow=c(2,2)); plot(r1,main="r1"); plot(r2,main="r2")
> >
> >
> > r3<-merge(r1,r2) ; plot(r3, main="merged rasters")
> >
> > # do something
> >
> > r3[r3>0]<-38     ; plot(r3, main="desired result")
> >
> >       [[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/
>
>
>
>
>
>
        [[alternative HTML version deleted]]

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

Re: How to sum area of adjacent patches?

Fri, 11/23/2018 - 16:58
Hi,

I think you may need to use raster::clump() and table() to label your regions.  Check out an old message about clump...

https://stat.ethz.ch/pipermail/r-sig-geo/2017-January/025346.html

It's sort of cumbersome, but I think the following works...


library(raster)
r1<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=0, ymx=1)
r2<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=1, ymx=2)
r1[]<-c(rep(24,24), rep(0,76))
r2[]<-c(rep(0,86), rep(14,14))
par(mfrow=c(3,2)); plot(r1,main="r1"); plot(r2,main="r2")


r3<-merge(r1,r2) ; plot(r3, main="merged rasters")

r3 <- r3
r4[r4>0]<-38     ; plot(r4, main="desired result")

r5 <- clump(r3) ; plot(r5, main = 'clumped')
r5[is.na(r5)] <- 0

r5t <- table(getValues(r5))
for (i in names(r5t)[-1]){
        id <- as.numeric(i)
        print(id)
        r5[r5 == id] <- r5t[i]
}

plot(r5, main = 'relabeled with area')


Cheers,
Ben

P.S. Years ago I used a language called IDL.  It's histogram function has an output keyword called reverse_indices (more here https://www.harrisgeospatial.com/docs/HISTOGRAM.html#H_835179117_677188 and here http://www.idlcoyote.com/tips/histogram_tutorial.html )  - super useful!  Your situation just begs for that same solution in R.





> On Nov 23, 2018, at 5:13 PM, Hugo Costa <[hidden email]> wrote:
>
> Deal all,
>
> I'm merging several rasters representing patches of forest, and the pixel
> values represent the area of the patches. The patches' area should be
> summed along the edges of the merged rasters. I provide a reproducible
> example. Any help is welcome.
> Hugo
>
> library(raster)
> r1<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=0, ymx=1)
> r2<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=1, ymx=2)
> r1[]<-c(rep(24,24), rep(0,76))
> r2[]<-c(rep(0,86), rep(14,14))
> par(mfrow=c(2,2)); plot(r1,main="r1"); plot(r2,main="r2")
>
>
> r3<-merge(r1,r2) ; plot(r3, main="merged rasters")
>
> # do something
>
> r3[r3>0]<-38     ; plot(r3, main="desired result")
>
> [[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

How to sum area of adjacent patches?

Fri, 11/23/2018 - 16:13
Deal all,

I'm merging several rasters representing patches of forest, and the pixel
values represent the area of the patches. The patches' area should be
summed along the edges of the merged rasters. I provide a reproducible
example. Any help is welcome.
Hugo

library(raster)
r1<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=0, ymx=1)
r2<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=1, ymx=2)
r1[]<-c(rep(24,24), rep(0,76))
r2[]<-c(rep(0,86), rep(14,14))
par(mfrow=c(2,2)); plot(r1,main="r1"); plot(r2,main="r2")


r3<-merge(r1,r2) ; plot(r3, main="merged rasters")

# do something

r3[r3>0]<-38     ; plot(r3, main="desired result")

        [[alternative HTML version deleted]]

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

Re: How to change axis (plot map)?

Fri, 11/23/2018 - 15:27
Probably you're looking for arguments xaxp and yaxp, used like this:
plot(1:100,xaxp=c(0,100,4), yaxp=c(0,100,2))

Please look for more info in the help. Also, if you intend to make maps,
possibly you should spend some time getting familiar with packages, such as
ggplot2 or tmap.

Cheers
Hugo


Lara Silva <[hidden email]> escreveu no dia sexta, 23/11/2018
à(s) 14:31:

> Hello,
>
>  I am trying to change the scale/number of x-axis and y-axis (map plot -
> Figure 1).
>
> https://www.dropbox.com/sh/34nddj3405q3a6p/AACe-_23Wx7bZRETYIOYjsiQa?dl=0
> For example:
>
> x-axis: 465000,480000, 500000 (min. med. max values)
> The same idea for y-axis ....
>
> Code:
>
> plot(proj22PA$OB_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
> main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis)
>
>
> plot(proj22PA$PittOB_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
> main= "", xlab ="x.coords", ylab="y.coords",
> cex.axis=caxis,xlim=c(465000,500000,480000), ylim=c(4270000,4300000))
>
> Do not change ... The map has shifted. Probably, I need a new function or
> modified the argumente for link the information.
>
>
> Another hypothesis
>
> plot(proj22PA$Pittosporum_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
> main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis,
> xaxt="n",yaxt="n")
>
> axis(1, seq(450000,480000,500000)) - x.axis
> axis(2, seq(4270000,425000,4300000)) - y.axis
>
> Output - Figure 2. I can not add the "numbers".
>
> https://www.dropbox.com/sh/34nddj3405q3a6p/AACe-_23Wx7bZRETYIOYjsiQa?dl=0
>
> Is it possible change the scale of legend (200, 400, 600, 800)?
>
> Any suggestion for my doubts/problems?
>
> Thanks,
>
> Silva
>
>         [[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

How to change axis (plot map)?

Fri, 11/23/2018 - 07:31
Hello,

 I am trying to change the scale/number of x-axis and y-axis (map plot -
Figure 1).

https://www.dropbox.com/sh/34nddj3405q3a6p/AACe-_23Wx7bZRETYIOYjsiQa?dl=0
For example:

x-axis: 465000,480000, 500000 (min. med. max values)
The same idea for y-axis ....

Code:

plot(proj22PA$OB_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis)


plot(proj22PA$PittOB_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
main= "", xlab ="x.coords", ylab="y.coords",
cex.axis=caxis,xlim=c(465000,500000,480000), ylim=c(4270000,4300000))

Do not change ... The map has shifted. Probably, I need a new function or
modified the argumente for link the information.


Another hypothesis

plot(proj22PA$Pittosporum_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis,
xaxt="n",yaxt="n")

axis(1, seq(450000,480000,500000)) - x.axis
axis(2, seq(4270000,425000,4300000)) - y.axis

Output - Figure 2. I can not add the "numbers".

https://www.dropbox.com/sh/34nddj3405q3a6p/AACe-_23Wx7bZRETYIOYjsiQa?dl=0

Is it possible change the scale of legend (200, 400, 600, 800)?

Any suggestion for my doubts/problems?

Thanks,

Silva

        [[alternative HTML version deleted]]

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

Re: Problem with initGRASS

Fri, 11/23/2018 - 02:52
Jaime,

if the location of your GRASS installation is not known you may also
give the link2GI package a try.

try link2GI::linkGRASS7() for a "brute force" searching and linking.
This should setup your enviroment for using the newest GRASS version.

cheers Chris

On 23/11/2018 09:44, Roger Bivand wrote:
> On Fri, 23 Nov 2018, Jaime Burbano Girón wrote:
>
>> Hi everyone,
>>
>> I'm trying to use GRASS in R, but I can't get to initialize it. I'm
>> using
>> Ubuntu 18.04 and GRASS 7.4.2. This is the system information:
>>
>> Versión de GRASS: 7.4.2
>>
>> Revisión SVN de GRASS: exported
>>
>> Fecha de compilación: 2018-10-28
>>
>> Construir plataforma: x86_64-pc-linux-gnu
>>
>> GDAL: 2.2.3
>>
>> PROJ.4: 4.9.3
>>
>> GEOS: 3.6.2
>>
>> SQLite: 3.22.0
>>
>> Python: 2.7.15rc1
>>
>> wxPython: 3.0.2.0
>>
>> Plataforma: Linux-4.15.0-38-generic-x86_64-with-Ubuntu-18.04-bionic
>>
>>
>> I'm running: > initGRASS(gisBase="/usr/lib/grass74/bin",override=T)
>> And this is the error:
>
> No, ?initGRASS says:
>
> gisBase: The directory path to GRASS binaries and libraries
>
> and grass74/bin only contains (for me):
>
> $ ls bin
> grass74
>
> which is a Python script, ASCII text executable. Trying
>
> grass74/grass-7.4.2
>
> (no idea where your installer puts this):
>
> $ ls grass-7.4.2
> AUTHORS           contributors_extra.csv  fonts REQUIREMENTS.html
> bin               COPYING                 GPL.TXT  scripts
> CHANGES           demolocation            gui      share
> CITING            docs                    include  tools
> config.status     driver                  INSTALL  translators.csv
> contributors.csv  etc                     lib
>
> which contains the binaries in bin and the libraries in lib. The
> required string is what is set in the shell variable GISBASE when
> running GRASS interactively:
>
> echo $GISBASE
>
> hence the name of the argument. initGRASS() cannot discover this as
> GRASS is not running, and has not set its environment variables. I
> have added a test looking for a directory called bin in the directory
> given as the gisBase argument, check on
> https://r-forge.r-project.org/R/?group_id=2020 to see when rgrass7
> 0.1-13 (Rev: 68) has build status current, try that with your current
> setting, the correct setting and report back.
>
> Hope this clarifies,
>
> Roger
>
>
>>
>> sh: 1: g.gisenv: not foundsh: 1: g.gisenv: not foundsh: 1: g.gisenv:
>> not foundsh: 1: g.gisenv: not foundsh: 1: g.gisenv: not foundsh: 1:
>> g.version: not foundError in system(paste("g.version", get("addEXE",
>> envir = .GRASS_CACHE),  :
>>  error in running commandAdemás: Warning messages:1: In
>> system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
>> error in running command2: In system(paste(paste("g.gisenv",
>> get("addEXE", envir = .GRASS_CACHE),  :  error in running command3: In
>> system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
>> error in running command4: In system(paste(paste("g.gisenv",
>> get("addEXE", envir = .GRASS_CACHE),  :  error in running command5: In
>> system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
>> error in running command
>>
>>
>> I verified that g.gisenv and g.version are located in
>> /usr/lib/grass74/bin.
>> Any suggestion?
>>
>> Thanks in advance.
>>
>> Best,
>>
>> Jaime.
>>
>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Dr Christoph Reudenbach, Philipps-University of Marburg, Faculty of Geography, GIS and Environmental Modeling, Deutschhausstr. 10, D-35032 Marburg, fon: ++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web: gis-ma.org, giswerk.org, moc.environmentalinformatics-marburg.de


        [[alternative HTML version deleted]]

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

Pages