accessing spatial imagery via rgdal's web map tile service driver

classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|

accessing spatial imagery via rgdal's web map tile service driver

Anthony Fischbach
This post was updated on .
I see that WMS ("OGC Web Map Service") and WMTS ("OGC Web Mab Tile
Service") are supported within my rgdal installation , however, it is
unclear how to call for WMTS imagery.
Please advise on how to access geo-referenced imagery through a web map
tile service, such as described here
(
https://wiki.earthdata.nasa.gov/display/GIBS/Map+Library+Usage#expand-GDALBasics
)
for gdal using the current rgdal version (under R version 3.4.3 rgdal
version 1.2-16)

Please note that my motivation for accessing the WMTS geo-referenced imagery is so that I may directly manipulate the imagery, projecting and clipping it to my study area so that I may build small (~< 300Kb) map images for distribution to users working in low-bandwidth situations.  It would replace the MODIS subset download functions in this script: https://doi.org/10.5066/F7ZS2TRS
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Tony Fischbach, Wildlife Biologist
Walrus Research Program
Alaska Science Center
U.S. Geological Survey
4210 University Drive
Anchorage, AK 99508-4650

AFischbach@usgs.gov
http://alaska.usgs.gov/science/biology/walrus
Reply | Threaded
Open this post in threaded view
|

Re: accessing spatial imagery via rgdal's web map tile service driver

Roger Bivand
Administrator
On Fri, 26 Jan 2018, Fischbach, Anthony wrote:

> I see that WMS ("OGC Web Map Service") and WMTS ("OGC Web Mab Tile
> Service") are supported within my rgdal installation , however, it is
> unclear how to call for WMTS imagery.
> Please advise on how to access geo-referenced imagery through a web map
> tile service, such as described here
> (
> https://wiki.earthdata.nasa.gov/display/GIBS/Map+Library+Usage#expand-GDALBasics
> )
> for gdal using the current rgdal version (under R version 3.4.3 rgdal
> version 1.2-16)

Do what anyone using GDAL drivers would do, refer to the documentation, in
this case on http://www.gdal.org/frmt_wms.html.

Roger

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

--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
http://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
Reply | Threaded
Open this post in threaded view
|

Re: accessing spatial imagery via rgdal's web map tile service driver

Anthony Fischbach
The GDAL tile web map service appears to require an xml specification for
the TMS and a "window" specification indicating the bounds of the region to
be accessed, specified in the projected coordinate system using the
-projwin flag. I am uncertain how to pass this -projwin and its bounds to
the rgdal::readGDAL function.  The code snippet below stages the problem.
Please advise on how to acquire just the "window" of interest.

#Tile Web Mapping Service to earthdata imagery via rgdal
library(rgdal)  # r implementation of GDAL
require(raster) # r package for raster image manipulation, required for
raster capture of the readGDAL function result.
require(sp)     # r package for spatial class objects, including CRS
#
prj.geo<-CRS('+proj=longlat +datum=WGS84') # Geographic projection
prj.Polar_NSIDC <- CRS("+init=epsg:3413")
# (+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0
+datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0)
prj.StudyArea <- CRS("+proj=laea +lat_0=70 +lon_0=-170 +x_0=0 +y_0=0
+datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0") # study area
projection
##
xLimits <- sort(c(170, -135)) # define the study area polygon using
geographic coordinates (eastings WGS84)
yLimits <- sort(c(52, 75)) # (northings WGS84)
pts<-expand.grid(xLimits, yLimits)
names(pts)<-c('x','y')
coordinates(pts) <- ~x+y
proj4string(pts) <- prj.geo
pts.Polar_NSIDC<-spTransform(pts, prj.Polar_NSIDC)
bb <- round(bbox(pts.Polar_NSIDC)) #
#
yesterday <- Sys.Date() - 1
TileLevel <- 3
xml_specification <-
paste0('<GDAL_WMS>
    <Service name="TMS">
        <ServerUrl>
https://gibs.earthdata.nasa.gov/wmts/epsg3413/best/MODIS_Aqua_CorrectedReflectance_Bands721/default/
',
format(yesterday, "%Y-%m-%d"),
'/250m/${z}/${y}/${x}.jpg</ServerUrl>
    </Service>
    <DataWindow>
        <UpperLeftX>-4194304</UpperLeftX>
        <UpperLeftY>4194304</UpperLeftY>
        <LowerRightX>4194304</LowerRightX>
        <LowerRightY>-4194304</LowerRightY>
        <TileLevel>', TileLevel, '</TileLevel>
        <TileCountX>2</TileCountX>
        <TileCountY>2</TileCountY>
        <YOrigin>top</YOrigin>
    </DataWindow>
    <Projection>EPSG:3413</Projection>
    <BlockSizeX>512</BlockSizeX>
    <BlockSizeY>512</BlockSizeY>
    <BandsCount>3</BandsCount>
</GDAL_WMS>
')
#
cat(xml_specification)
GDALinfo(fname = xml_specification)
 (projWinParameters = paste('-projwin', bb[1, 1], bb[2,2], bb[1,2],
bb[2,1], sep = " "))
aqua <- readGDAL(fname = xml_specification) # , projWinParameters) ?? how
to implement the projwin constraints??
(aqua_stack <- stack(aqua))
#
plotRGB(aqua_stack) # plots the acquired image in RGB !! problem is that
the entire WMS domain is acquired.

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Tony Fischbach, Wildlife Biologist
Walrus Research Program
Alaska Science Center
U.S. Geological Survey
4210 University Drive
Anchorage, AK 99508-4650

AFischbach@usgs.gov
http://alaska.usgs.gov/science/biology/walrus
Reply | Threaded
Open this post in threaded view
|

Re: accessing spatial imagery via rgdal's web map tile service driver

Michael Sumner-2
Hi Anthony, what OS are you on? Can you share sessionInfo() or
devtools::session_info()?

But, it could be simpler, but  it doesn't work on my system, I'll find out
more:

https://gist.github.com/mdsumner/ca332f9c5112e00edcc40fdcf2a4b968

Cheers, Mike.

On Thu, 8 Feb 2018 at 06:28 Fischbach, Anthony <[hidden email]> wrote:

> The GDAL tile web map service appears to require an xml specification for
> the TMS and a "window" specification indicating the bounds of the region to
> be accessed, specified in the projected coordinate system using the
> -projwin flag. I am uncertain how to pass this -projwin and its bounds to
> the rgdal::readGDAL function.  The code snippet below stages the problem.
> Please advise on how to acquire just the "window" of interest.
>
> #Tile Web Mapping Service to earthdata imagery via rgdal
> library(rgdal)  # r implementation of GDAL
> require(raster) # r package for raster image manipulation, required for
> raster capture of the readGDAL function result.
> require(sp)     # r package for spatial class objects, including CRS
> #
> prj.geo<-CRS('+proj=longlat +datum=WGS84') # Geographic projection
> prj.Polar_NSIDC <- CRS("+init=epsg:3413")
> # (+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0
> +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0)
> prj.StudyArea <- CRS("+proj=laea +lat_0=70 +lon_0=-170 +x_0=0 +y_0=0
> +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0") # study area
> projection
> ##
> xLimits <- sort(c(170, -135)) # define the study area polygon using
> geographic coordinates (eastings WGS84)
> yLimits <- sort(c(52, 75)) # (northings WGS84)
> pts<-expand.grid(xLimits, yLimits)
> names(pts)<-c('x','y')
> coordinates(pts) <- ~x+y
> proj4string(pts) <- prj.geo
> pts.Polar_NSIDC<-spTransform(pts, prj.Polar_NSIDC)
> bb <- round(bbox(pts.Polar_NSIDC)) #
> #
> yesterday <- Sys.Date() - 1
> TileLevel <- 3
> xml_specification <-
> paste0('<GDAL_WMS>
>     <Service name="TMS">
>         <ServerUrl>
>
> https://gibs.earthdata.nasa.gov/wmts/epsg3413/best/MODIS_Aqua_CorrectedReflectance_Bands721/default/
> ',
> format(yesterday, "%Y-%m-%d"),
> '/250m/${z}/${y}/${x}.jpg</ServerUrl>
>     </Service>
>     <DataWindow>
>         <UpperLeftX>-4194304</UpperLeftX>
>         <UpperLeftY>4194304</UpperLeftY>
>         <LowerRightX>4194304</LowerRightX>
>         <LowerRightY>-4194304</LowerRightY>
>         <TileLevel>', TileLevel, '</TileLevel>
>         <TileCountX>2</TileCountX>
>         <TileCountY>2</TileCountY>
>         <YOrigin>top</YOrigin>
>     </DataWindow>
>     <Projection>EPSG:3413</Projection>
>     <BlockSizeX>512</BlockSizeX>
>     <BlockSizeY>512</BlockSizeY>
>     <BandsCount>3</BandsCount>
> </GDAL_WMS>
> ')
> #
> cat(xml_specification)
> GDALinfo(fname = xml_specification)
>  (projWinParameters = paste('-projwin', bb[1, 1], bb[2,2], bb[1,2],
> bb[2,1], sep = " "))
> aqua <- readGDAL(fname = xml_specification) # , projWinParameters) ?? how
> to implement the projwin constraints??
> (aqua_stack <- stack(aqua))
> #
> plotRGB(aqua_stack) # plots the acquired image in RGB !! problem is that
> the entire WMS domain is acquired.
>
>         [[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
Reply | Threaded
Open this post in threaded view
|

Re: accessing spatial imagery via rgdal's web map tile service driver

Anthony Fischbach
Mike, thank you for digging in on this.  I have submitted my sessionInfo()
and comments to your git script regarding an error on the raster function
call. - Tony

Anthony Fischbach, Wildlife Biologist
USGS Alaska Science Center Walrus Research Program
4210 University Drive
Anchorage, AK 99508

voice: (907) 786-7145

http://alaska.usgs.gov/science/biology/walrus/

On Thu, Feb 8, 2018 at 12:26 PM, Michael Sumner <[hidden email]> wrote:

> Hi Anthony, what OS are you on? Can you share sessionInfo() or
> devtools::session_info()?
>
> But, it could be simpler, but  it doesn't work on my system, I'll find out
> more:
>
> https://gist.github.com/mdsumner/ca332f9c5112e00edcc40fdcf2a4b968
>
> Cheers, Mike.
>
> On Thu, 8 Feb 2018 at 06:28 Fischbach, Anthony <[hidden email]>
> wrote:
>
>> The GDAL tile web map service appears to require an xml specification for
>> the TMS and a "window" specification indicating the bounds of the region
>> to
>> be accessed, specified in the projected coordinate system using the
>> -projwin flag. I am uncertain how to pass this -projwin and its bounds to
>> the rgdal::readGDAL function.  The code snippet below stages the problem.
>> Please advise on how to acquire just the "window" of interest.
>>
>> #Tile Web Mapping Service to earthdata imagery via rgdal
>> library(rgdal)  # r implementation of GDAL
>> require(raster) # r package for raster image manipulation, required for
>> raster capture of the readGDAL function result.
>> require(sp)     # r package for spatial class objects, including CRS
>> #
>> prj.geo<-CRS('+proj=longlat +datum=WGS84') # Geographic projection
>> prj.Polar_NSIDC <- CRS("+init=epsg:3413")
>> # (+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0
>> +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0)
>> prj.StudyArea <- CRS("+proj=laea +lat_0=70 +lon_0=-170 +x_0=0 +y_0=0
>> +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0") # study area
>> projection
>> ##
>> xLimits <- sort(c(170, -135)) # define the study area polygon using
>> geographic coordinates (eastings WGS84)
>> yLimits <- sort(c(52, 75)) # (northings WGS84)
>> pts<-expand.grid(xLimits, yLimits)
>> names(pts)<-c('x','y')
>> coordinates(pts) <- ~x+y
>> proj4string(pts) <- prj.geo
>> pts.Polar_NSIDC<-spTransform(pts, prj.Polar_NSIDC)
>> bb <- round(bbox(pts.Polar_NSIDC)) #
>> #
>> yesterday <- Sys.Date() - 1
>> TileLevel <- 3
>> xml_specification <-
>> paste0('<GDAL_WMS>
>>     <Service name="TMS">
>>         <ServerUrl>
>> https://gibs.earthdata.nasa.gov/wmts/epsg3413/best/MODIS_
>> Aqua_CorrectedReflectance_Bands721/default/
>> ',
>> format(yesterday, "%Y-%m-%d"),
>> '/250m/${z}/${y}/${x}.jpg</ServerUrl>
>>     </Service>
>>     <DataWindow>
>>         <UpperLeftX>-4194304</UpperLeftX>
>>         <UpperLeftY>4194304</UpperLeftY>
>>         <LowerRightX>4194304</LowerRightX>
>>         <LowerRightY>-4194304</LowerRightY>
>>         <TileLevel>', TileLevel, '</TileLevel>
>>         <TileCountX>2</TileCountX>
>>         <TileCountY>2</TileCountY>
>>         <YOrigin>top</YOrigin>
>>     </DataWindow>
>>     <Projection>EPSG:3413</Projection>
>>     <BlockSizeX>512</BlockSizeX>
>>     <BlockSizeY>512</BlockSizeY>
>>     <BandsCount>3</BandsCount>
>> </GDAL_WMS>
>> ')
>> #
>> cat(xml_specification)
>> GDALinfo(fname = xml_specification)
>>  (projWinParameters = paste('-projwin', bb[1, 1], bb[2,2], bb[1,2],
>> bb[2,1], sep = " "))
>> aqua <- readGDAL(fname = xml_specification) # , projWinParameters) ?? how
>> to implement the projwin constraints??
>> (aqua_stack <- stack(aqua))
>> #
>> plotRGB(aqua_stack) # plots the acquired image in RGB !! problem is that
>> the entire WMS domain is acquired.
>>
>>         [[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
> <https://maps.google.com/?q=203+Channel+Highway+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
> Kingston Tasmania 7050 Australia
> <https://maps.google.com/?q=203+Channel+Highway+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Tony Fischbach, Wildlife Biologist
Walrus Research Program
Alaska Science Center
U.S. Geological Survey
4210 University Drive
Anchorage, AK 99508-4650

AFischbach@usgs.gov
http://alaska.usgs.gov/science/biology/walrus
Reply | Threaded
Open this post in threaded view
|

Re: accessing spatial imagery via rgdal's web map tile service driver

Michael Sumner-2
Oh, I had a typo in the xml_specification due to blithe copy/paste from
email. It works fine as intended.

I'll update and report back, this is a good example for fleshing out some
issues.

Cheers!


On Fri, 9 Feb 2018 at 08:46 Fischbach, Anthony <[hidden email]> wrote:

> Mike, thank you for digging in on this.  I have submitted my sessionInfo()
> and comments to your git script regarding an error on the raster function
> call. - Tony
>
> Anthony Fischbach, Wildlife Biologist
> USGS Alaska Science Center Walrus Research Program
> 4210 University Drive
> <https://maps.google.com/?q=4210+University+DriveAnchorage,+AK+99508&entry=gmail&source=g>
> Anchorage, AK 99508
> <https://maps.google.com/?q=4210+University+DriveAnchorage,+AK+99508&entry=gmail&source=g>
>
> voice: (907) 786-7145
>
> http://alaska.usgs.gov/science/biology/walrus/
>
> On Thu, Feb 8, 2018 at 12:26 PM, Michael Sumner <[hidden email]>
> wrote:
>
>> Hi Anthony, what OS are you on? Can you share sessionInfo() or
>> devtools::session_info()?
>>
>> But, it could be simpler, but  it doesn't work on my system, I'll find
>> out more:
>>
>> https://gist.github.com/mdsumner/ca332f9c5112e00edcc40fdcf2a4b968
>>
>> Cheers, Mike.
>>
>> On Thu, 8 Feb 2018 at 06:28 Fischbach, Anthony <[hidden email]>
>> wrote:
>>
>>> The GDAL tile web map service appears to require an xml specification for
>>> the TMS and a "window" specification indicating the bounds of the region
>>> to
>>> be accessed, specified in the projected coordinate system using the
>>> -projwin flag. I am uncertain how to pass this -projwin and its bounds to
>>> the rgdal::readGDAL function.  The code snippet below stages the problem.
>>> Please advise on how to acquire just the "window" of interest.
>>>
>>> #Tile Web Mapping Service to earthdata imagery via rgdal
>>> library(rgdal)  # r implementation of GDAL
>>> require(raster) # r package for raster image manipulation, required for
>>> raster capture of the readGDAL function result.
>>> require(sp)     # r package for spatial class objects, including CRS
>>> #
>>> prj.geo<-CRS('+proj=longlat +datum=WGS84') # Geographic projection
>>> prj.Polar_NSIDC <- CRS("+init=epsg:3413")
>>> # (+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0
>>> +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0)
>>> prj.StudyArea <- CRS("+proj=laea +lat_0=70 +lon_0=-170 +x_0=0 +y_0=0
>>> +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0") # study area
>>> projection
>>> ##
>>> xLimits <- sort(c(170, -135)) # define the study area polygon using
>>> geographic coordinates (eastings WGS84)
>>> yLimits <- sort(c(52, 75)) # (northings WGS84)
>>> pts<-expand.grid(xLimits, yLimits)
>>> names(pts)<-c('x','y')
>>> coordinates(pts) <- ~x+y
>>> proj4string(pts) <- prj.geo
>>> pts.Polar_NSIDC<-spTransform(pts, prj.Polar_NSIDC)
>>> bb <- round(bbox(pts.Polar_NSIDC)) #
>>> #
>>> yesterday <- Sys.Date() - 1
>>> TileLevel <- 3
>>> xml_specification <-
>>> paste0('<GDAL_WMS>
>>>     <Service name="TMS">
>>>         <ServerUrl>
>>>
>>> https://gibs.earthdata.nasa.gov/wmts/epsg3413/best/MODIS_Aqua_CorrectedReflectance_Bands721/default/
>>> ',
>>> format(yesterday, "%Y-%m-%d"),
>>> '/250m/${z}/${y}/${x}.jpg</ServerUrl>
>>>     </Service>
>>>     <DataWindow>
>>>         <UpperLeftX>-4194304</UpperLeftX>
>>>         <UpperLeftY>4194304</UpperLeftY>
>>>         <LowerRightX>4194304</LowerRightX>
>>>         <LowerRightY>-4194304</LowerRightY>
>>>         <TileLevel>', TileLevel, '</TileLevel>
>>>         <TileCountX>2</TileCountX>
>>>         <TileCountY>2</TileCountY>
>>>         <YOrigin>top</YOrigin>
>>>     </DataWindow>
>>>     <Projection>EPSG:3413</Projection>
>>>     <BlockSizeX>512</BlockSizeX>
>>>     <BlockSizeY>512</BlockSizeY>
>>>     <BandsCount>3</BandsCount>
>>> </GDAL_WMS>
>>> ')
>>> #
>>> cat(xml_specification)
>>> GDALinfo(fname = xml_specification)
>>>  (projWinParameters = paste('-projwin', bb[1, 1], bb[2,2], bb[1,2],
>>> bb[2,1], sep = " "))
>>> aqua <- readGDAL(fname = xml_specification) # , projWinParameters) ?? how
>>> to implement the projwin constraints??
>>> (aqua_stack <- stack(aqua))
>>> #
>>> plotRGB(aqua_stack) # plots the acquired image in RGB !! problem is that
>>> the entire WMS domain is acquired.
>>>
>>>         [[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
>> <https://maps.google.com/?q=203+Channel+Highway+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
>> Kingston Tasmania 7050 Australia
>> <https://maps.google.com/?q=203+Channel+Highway+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
>>
>>
> --
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
Reply | Threaded
Open this post in threaded view
|

Re: accessing spatial imagery via rgdal's web map tile service driver

Michael Sumner-2
I've put a fuller example together, thanks for the prompt, it's nice to
know how to do this:

http://rpubs.com/cyclemumner/358029

Note that raster is using underlying rgdal offset/output.dim functionality
to drive GDAL itself here, something that's not exposed well - raster makes
it seamless, but the underlying tools are a bit obscure and could be used
more widely.

The example uses row/col indexing via the extent(RasterLayer, numeric, ...)
method, you can otherwise use any extent in the native (north Polar
Stereographic) coordinate system. I don't know if raster interprets
aggregate(, fact = ) calls into translations for rgdal's output.dim, but it
would be powerful to be easy to do that.

Cheers, Mike.

On Fri, 9 Feb 2018 at 10:32 Michael Sumner <[hidden email]> wrote:

> Oh, I had a typo in the xml_specification due to blithe copy/paste from
> email. It works fine as intended.
>
> I'll update and report back, this is a good example for fleshing out some
> issues.
>
> Cheers!
>
>
> On Fri, 9 Feb 2018 at 08:46 Fischbach, Anthony <[hidden email]>
> wrote:
>
>> Mike, thank you for digging in on this.  I have submitted my
>> sessionInfo() and comments to your git script regarding an error on the
>> raster function call. - Tony
>>
>> Anthony Fischbach, Wildlife Biologist
>> USGS Alaska Science Center Walrus Research Program
>> 4210 University Drive
>> <https://maps.google.com/?q=4210+University+DriveAnchorage,+AK+99508&entry=gmail&source=g>
>> Anchorage, AK 99508
>> <https://maps.google.com/?q=4210+University+DriveAnchorage,+AK+99508&entry=gmail&source=g>
>>
>> voice: (907) 786-7145
>>
>> http://alaska.usgs.gov/science/biology/walrus/
>>
>> On Thu, Feb 8, 2018 at 12:26 PM, Michael Sumner <[hidden email]>
>> wrote:
>>
>>> Hi Anthony, what OS are you on? Can you share sessionInfo() or
>>> devtools::session_info()?
>>>
>>> But, it could be simpler, but  it doesn't work on my system, I'll find
>>> out more:
>>>
>>> https://gist.github.com/mdsumner/ca332f9c5112e00edcc40fdcf2a4b968
>>>
>>> Cheers, Mike.
>>>
>>> On Thu, 8 Feb 2018 at 06:28 Fischbach, Anthony <[hidden email]>
>>> wrote:
>>>
>>>> The GDAL tile web map service appears to require an xml specification
>>>> for
>>>> the TMS and a "window" specification indicating the bounds of the
>>>> region to
>>>> be accessed, specified in the projected coordinate system using the
>>>> -projwin flag. I am uncertain how to pass this -projwin and its bounds
>>>> to
>>>> the rgdal::readGDAL function.  The code snippet below stages the
>>>> problem.
>>>> Please advise on how to acquire just the "window" of interest.
>>>>
>>>> #Tile Web Mapping Service to earthdata imagery via rgdal
>>>> library(rgdal)  # r implementation of GDAL
>>>> require(raster) # r package for raster image manipulation, required for
>>>> raster capture of the readGDAL function result.
>>>> require(sp)     # r package for spatial class objects, including CRS
>>>> #
>>>> prj.geo<-CRS('+proj=longlat +datum=WGS84') # Geographic projection
>>>> prj.Polar_NSIDC <- CRS("+init=epsg:3413")
>>>> # (+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0
>>>> +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0)
>>>> prj.StudyArea <- CRS("+proj=laea +lat_0=70 +lon_0=-170 +x_0=0 +y_0=0
>>>> +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0") # study
>>>> area
>>>> projection
>>>> ##
>>>> xLimits <- sort(c(170, -135)) # define the study area polygon using
>>>> geographic coordinates (eastings WGS84)
>>>> yLimits <- sort(c(52, 75)) # (northings WGS84)
>>>> pts<-expand.grid(xLimits, yLimits)
>>>> names(pts)<-c('x','y')
>>>> coordinates(pts) <- ~x+y
>>>> proj4string(pts) <- prj.geo
>>>> pts.Polar_NSIDC<-spTransform(pts, prj.Polar_NSIDC)
>>>> bb <- round(bbox(pts.Polar_NSIDC)) #
>>>> #
>>>> yesterday <- Sys.Date() - 1
>>>> TileLevel <- 3
>>>> xml_specification <-
>>>> paste0('<GDAL_WMS>
>>>>     <Service name="TMS">
>>>>         <ServerUrl>
>>>>
>>>> https://gibs.earthdata.nasa.gov/wmts/epsg3413/best/MODIS_Aqua_CorrectedReflectance_Bands721/default/
>>>> ',
>>>> format(yesterday, "%Y-%m-%d"),
>>>> '/250m/${z}/${y}/${x}.jpg</ServerUrl>
>>>>     </Service>
>>>>     <DataWindow>
>>>>         <UpperLeftX>-4194304</UpperLeftX>
>>>>         <UpperLeftY>4194304</UpperLeftY>
>>>>         <LowerRightX>4194304</LowerRightX>
>>>>         <LowerRightY>-4194304</LowerRightY>
>>>>         <TileLevel>', TileLevel, '</TileLevel>
>>>>         <TileCountX>2</TileCountX>
>>>>         <TileCountY>2</TileCountY>
>>>>         <YOrigin>top</YOrigin>
>>>>     </DataWindow>
>>>>     <Projection>EPSG:3413</Projection>
>>>>     <BlockSizeX>512</BlockSizeX>
>>>>     <BlockSizeY>512</BlockSizeY>
>>>>     <BandsCount>3</BandsCount>
>>>> </GDAL_WMS>
>>>> ')
>>>> #
>>>> cat(xml_specification)
>>>> GDALinfo(fname = xml_specification)
>>>>  (projWinParameters = paste('-projwin', bb[1, 1], bb[2,2], bb[1,2],
>>>> bb[2,1], sep = " "))
>>>> aqua <- readGDAL(fname = xml_specification) # , projWinParameters) ??
>>>> how
>>>> to implement the projwin constraints??
>>>> (aqua_stack <- stack(aqua))
>>>> #
>>>> plotRGB(aqua_stack) # plots the acquired image in RGB !! problem is that
>>>> the entire WMS domain is acquired.
>>>>
>>>>         [[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
>>> <https://maps.google.com/?q=203+Channel+Highway+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
>>> Kingston Tasmania 7050 Australia
>>> <https://maps.google.com/?q=203+Channel+Highway+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
>>>
>>>
>> --
> 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