[GRASS-user] WCS import into GRASS GIS - select region extent only

Sat Jan 4 14:32:08 PST 2014

Hi Stefan, hi Markus,
I just started to think over everything again after Stefans first
maill... thanks for the fast reaction. I just tried the proposed
HEIGHT and WIDTH parameter with gdal_translate without success.

I started to observe the HTTP-transfer with mitmproxy [1]. I can see,
that BBOX parameter is just ignored. Also i can see, that r.in.gdal by
itself seems to know about tile-size-limitations and starts to
download small 2MB-tiles, while gdal_translate just tries to get it
all in one file. If you are interested in more details, I can share
these observations on the new WCS-wikipage?

Concerning the "limit import to current region flag" I agree
absolutely with you two, I would just take it with a handkiss. Would
this -r flag be something easy and quick to implement? I see that
r.in.gdal is a compiled binary, so nothing I believe I can do myself.

Thanks to for the python-based-approach, which hopefully will respect
tile-size limits. I'll give it a try, until the -r flag will be
reality.

Cheers,
Martin

[1] http://mitmproxy.org/
--
Martin Zbinden
Riedacker 523
3154 Rüschegg Heubach
+41 78 628 28 82

2014/1/4 Blumentrath, Stefan <Stefan.Blumentrath at nina.no>:
> Hi Martin,
> Since you mentioned that you are familiar with python, here is a python-based approach:
> http://gis.stackexchange.com/questions/41123/how-to-extract-a-wcs-layer-with-gdal-translate
>
> Cheers
> Stefan
>
> -----Original Message-----
> From: grass-user-bounces at lists.osgeo.org [mailto:grass-user-bounces at lists.osgeo.org] On Behalf Of Blumentrath, Stefan
> Sent: 4. januar 2014 22:17
> To: Martin Zbinden; grass-user at lists.osgeo.org
> Subject: Re: [GRASS-user] WCS import into GRASS GIS - select region extent only
>
> Hi Martin,
>
> Just add a line like this (adjusted to your case) after the CoverageName line in the .vrt.
> <GetCoverageExtra>&SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&FORMAT=GeoTIFF&COVERAGE=layername&BBOX=381464.564874,7282848.63791,7310018.58218,422460.634623&CRS=EPSG:32633&WIDTH=800&HEIGHT=600</GetCoverageExtra>
>
> Likely you do not need all these GetCoverage specifications.
>
> However, with the BBOX set to your region you can import the .vrt directly to GRASS with r.in.gdal.
>
> (See http://www.gdal.org/frmt_wcs.html or http://wiki.esipfed.org/images/2/2c/BestPractices_WCS.pdf page 24)
>
> Cheers
> Stefan
>
>
> -----Original Message-----
> From: martin.zbinden at gmail.com [mailto:martin.zbinden at gmail.com] On Behalf Of Martin Zbinden
> Sent: 4. januar 2014 21:01
> To: Blumentrath, Stefan
> Subject: Re: [GRASS-user] WCS import into GRASS GIS - select region extent only
>
> Hey Stefan,
> Thanks, you helped me already a lot. I'm looking forward to build some python script which does these steps automatically. It seems to me that these .vrt-files can do kind of magic :-)
>
> However, one concern is still left. When I choose a too big projwin, then Mapserver gives an error.
>
> gdal_translate -co "COMPRESS=LZW" -projwin 2592000 1182500 2597500
> 1179000  dem_wcs.vrt dem_wcs.tif
> Computed -srcwin 53295 56717 2750 1750 from projected window.
>
> 0ERROR 1: msWCSGetCoverage(): WCS server error. Raster size out of range, width and height of resulting coverage must be no more than MAXSIZE=2048.
>
> It seems that r.in.gdal ist conscious about this and downloads data tile by tile. Maybe you can give me a hint me how to do the same with my own script. No way to tell the projwin through the GDAL-VRT directly, so I could load it directly with r.in.gdal?
>
> Cheers,
> Martin
> --
> Martin Zbinden
> Riedacker 523
> 3154 Rüschegg Heubach
> +41 78 628 28 82
>
>
> 2014/1/3 Blumentrath, Stefan <Stefan.Blumentrath at nina.no>:
>> Hei Martin,
>>
>> You could try making a (temporary) virtual GDAL-VRT dataset (it’s a simple textfile (XML) containing information like this, in the example below named "dem_wcs.vrt"):
>>
>> <WCS_GDAL>
>>  <ServiceURL>http://your_wcs_server/wcs.dtm?</ServiceURL>
>>  <CoverageName>layername</CoverageName>
>> </WCS_GDAL>
>>
>> For more info on vrt-format see:
>> http://www.gdal.org/gdal_vrttut.html
>>
>> After that you could use gdal_translate (http://www.gdal.org/gdal_translate.html), which has a projwin parameter where you can feed you GRASS region boundaries from g.region -up:
>> gdal_translate -co "COMPRESS=LZW" -projwin ulx uly lrx lry dem_wcs.vrt
>> dem_wcs.tif
>>
>> Hope that helps,
>>
>> Cheers
>> Stefan
>>
>> -----Original Message-----
>> From: grass-user-bounces at lists.osgeo.org
>> [mailto:grass-user-bounces at lists.osgeo.org] On Behalf Of Martin
>> Zbinden
>> Sent: 3. januar 2014 21:46
>> To: grass-user at lists.osgeo.org
>> Subject: [GRASS-user] WCS import into GRASS GIS - select region extent
>> only
>>
>> Hi,
>>
>> For my Bachelor-Thesis I want to optimize the erosion risk map of Switzerland for use in agricultural extension. Although I concentrate on 3 cantons only (BE, FR, SO),  I have to manage huge amounts of data. The fact that the raster data I got from the federal office are overlapping each other (following the defined fieldblocks) doesn't make things easier, but I've at least found r.patch which can combine them one by one.
>>
>> I thought, with Mapserver  (V 6.4.0) I'd found the perfect solution to all my problems. Now I'm stuck after 5+ hours of trying to import the data from my local WCS-Service. I have seen the hints on GRASS Wiki [1] and r.in.gdal starts to download geotiff-tiles from Mapserver. The problem is, that I don't want all the 12 GB of data, but only the grass region extent.
>>
>> The following http-query does what I expect, it makes me a GeoTIFF-File from my region defined by BBOX=... .
>> http://localhost/mywcs?SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&c
>> overage=blw.eros_z&FORMAT=image/tiff&CRS=EPSG:2056&BBOX=2596000,118100
>> 0,2597500,1182500&RESX=2&RESY=2
>>
>>
>> I tried to restrict the download to the  region by defining a BBOX in this WCS_GDAL xml file. However, this didn't work and isn't supposed to work [2]. Is there a way to say r.in.gdal to only get the data from the region extent?
>>
>> I know there is this new, very clever r.in.wms2 script, which does exactly what I want -- for WMS. Is there a way to get WMS-Mapserver to serve RAW-data (erosion in t per annum, elevation data)?  Or maybe better create a python script to download map by a custom http-request?
>>
>> Thanks in advance for any advice! I hope I can get it to work somehow.
>> Liebe Grüsse
>> Martin Zbinden
>>
>> PS: the WCS server given in GRASS wiki seems to be down this time.
>>
>> [1]
>> http://grasswiki.osgeo.org/wiki/Global_datasets#OGC_WCS_-_Albedo_examp
>> le [2]
>> http://lists.osgeo.org/pipermail/gdal-dev/2008-March/016552.html
>>
>>
>> --
>> Martin Zbinden
>> Riedacker 523
>> 3154 Rüschegg Heubach
>> +41 78 628 28 82
>> _______________________________________________
>> grass-user mailing list
>> grass-user at lists.osgeo.org
>> http://lists.osgeo.org/mailman/listinfo/grass-user
>> _______________________________________________
>> grass-user mailing list
>> grass-user at lists.osgeo.org
>> http://lists.osgeo.org/mailman/listinfo/grass-user
> _______________________________________________
> grass-user mailing list
> grass-user at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-user
> _______________________________________________
> grass-user mailing list
> grass-user at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-user