Re: [GRASS-user] WCS import into GRASS GIS - select region extent only
On 05/01/14 22:48, Blumentrath, Stefan wrote: As for the r-flag Markus mentioned: personally I would not miss it in r.external because data is not really duplicated, or are there drawbacks when registering big external datasets? In r.in.gdal it might be useful, but does the module change the resolution or would it do so if a bounding box is applied to the input-raster? One option would be to have the -r flag only apply to the bounding box and not to the resolution of the current region, but adjusting that bounding box outwards as much as necessary to fit the resolution of the imported raster. Moritz ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] WCS import into GRASS GIS - select region extent only
Sorry, there were some line-breaks missing in the script. It should look like
this:
source_prj=`eval g.proj -fj`
r.tileset sourceproj="$source_prj" maxcols=1024 maxrows=1024 fs=' ' --q | awk
'{print "gdalbuildvrt ./dem_wcs_bb_"NR".vrt ./dem_wcs.vrt -te " $1, $2, $3, $4
}' | sh
gdalbuildvrt vrtmosaic.vrt dem_wcs_bb_*.vrt
r.in.gdal input=vrtmosaic.vrt output=test_wcs
-Original Message-
From: grass-user-boun...@lists.osgeo.org
[mailto:grass-user-boun...@lists.osgeo.org] On Behalf Of Blumentrath, Stefan
Sent: 5. januar 2014 22:48
To: Martin Zbinden; grass-user@lists.osgeo.org
Subject: Re: [GRASS-user] WCS import into GRASS GIS - select region extent only
Hi Martin and Markus,
The GDAL VRT-format is really very handy for many kinds of data extraction, so
I simply applied one of my most important learning-strategies (not only for
GIS): "try-and-error"
However, I have to admit, that I did not really understand your last question,
Martin:
What do you mean by a "bunch of vrt" if you only have 4 maps? Did you set -te
according to the size limits of the server and therefor get maybe hundrets of
tiles? And if r.in.gdal allows you to import more than WCS-limits for HEIGHT
and WIDTH, did you try to import your entire region of interest through the
second vrt?
Alternatively you could try something like this (Unix shell script for GRASS
6.4 (in GRASS 7 r.tileset seems to be still broken, see
http://trac.osgeo.org/grass/ticket/1825):
This solution assumes that you set your GRASS region to your entire area of
interest...
source_prj=`eval g.proj -fj`
r.tileset sourceproj="$source_prj" maxcols=1024 maxrows=1024 fs=' ' --q | awk
'{print "gdalbuildvrt ./dem_wcs_bb_"NR".vrt ./dem_wcs.vrt -te " $1, $2, $3, $4
}' | sh gdalbuildvrt vrtmosaic.vrt dem_wcs_bb_*.vrt r.in.gdal
input=vrtmosaic.vrt output=test_wcs
Please note that the WCS I used for testing probably does not have any
tile-limits because I am having no issues (on UBUNTU 12.04 and GRASS 7 / GRASS
6.4.4).
As for the r-flag Markus mentioned: personally I would not miss it in
r.external because data is not really duplicated, or are there drawbacks when
registering big external datasets?
In r.in.gdal it might be useful, but does the module change the resolution or
would it do so if a bounding box is applied to the input-raster?
Cheers
Stefan
-Original Message-
From: martin.zbin...@gmail.com [mailto:martin.zbin...@gmail.com] On Behalf Of
Martin Zbinden
Sent: 5. januar 2014 00:10
To: Blumentrath, Stefan
Subject: Re: [GRASS-user] WCS import into GRASS GIS - select region extent only
Wow, great! This trick with a second VRT works with gdal_translate and
r.in.gdal, with the second even in case of regions bigger than MAXSIZE. How did
you find that out, if I may ask you?
What do you think: I have to handle at least 4 raster maps. What is the most
efficient way to organize this bunch of VRTs?
I think, I could base on the function _createXML() from r.in.wms2 an generate
the VRT myself. The second VRT introduces the following important lines
depending on TE parameter:
2.5920e+06, 2.e+00, 0.e+00,
1.1825e+06, 0.e+00,
-2.e+00
Cheers and good night
Martin
--
Martin Zbinden
Riedacker 523
3154 Rüschegg Heubach
+41 78 628 28 82
2014/1/4 Blumentrath, Stefan :
> Actually, a second vrt does the trick:
>
> gdalbuild vrt dem_wcs_bb.vrt dem_wcs.vrt -te 381464.564874
> 7282848.63791 422460.634623 7310018.58218
>
> This worked for me.
> I shall add the solution to the Wiki.
>
> Cheers
> Stefan
>
> -Original Message-
> From: martin.zbin...@gmail.com [mailto:martin.zbin...@gmail.com] On
> Behalf Of Martin Zbinden
> Sent: 4. januar 2014 23:32
> To: Blumentrath, Stefan; nete...@osgeo.org
> Cc: grass-user@lists.osgeo.org
> Subject: Re: [GRASS-user] WCS import into GRASS GIS - select region
> extent only
>
> 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
Re: [GRASS-user] WCS import into GRASS GIS - select region extent only
Hi Martin and Markus,
The GDAL VRT-format is really very handy for many kinds of data extraction, so
I simply applied one of my most important learning-strategies (not only for
GIS): "try-and-error"
However, I have to admit, that I did not really understand your last question,
Martin:
What do you mean by a "bunch of vrt" if you only have 4 maps? Did you set -te
according to the size limits of the server and therefor get maybe hundrets of
tiles? And if r.in.gdal allows you to import more than WCS-limits for HEIGHT
and WIDTH, did you try to import your entire region of interest through the
second vrt?
Alternatively you could try something like this (Unix shell script for GRASS
6.4 (in GRASS 7 r.tileset seems to be still broken, see
http://trac.osgeo.org/grass/ticket/1825):
This solution assumes that you set your GRASS region to your entire area of
interest...
source_prj=`eval g.proj -fj`
r.tileset sourceproj="$source_prj" maxcols=1024 maxrows=1024 fs=' ' --q | awk
'{print "gdalbuildvrt ./dem_wcs_bb_"NR".vrt ./dem_wcs.vrt -te " $1, $2, $3, $4
}' | sh
gdalbuildvrt vrtmosaic.vrt dem_wcs_bb_*.vrt
r.in.gdal input=vrtmosaic.vrt output=test_wcs
Please note that the WCS I used for testing probably does not have any
tile-limits because I am having no issues (on UBUNTU 12.04 and GRASS 7 / GRASS
6.4.4).
As for the r-flag Markus mentioned: personally I would not miss it in
r.external because data is not really duplicated, or are there drawbacks when
registering big external datasets?
In r.in.gdal it might be useful, but does the module change the resolution or
would it do so if a bounding box is applied to the input-raster?
Cheers
Stefan
-Original Message-
From: martin.zbin...@gmail.com [mailto:martin.zbin...@gmail.com] On Behalf Of
Martin Zbinden
Sent: 5. januar 2014 00:10
To: Blumentrath, Stefan
Subject: Re: [GRASS-user] WCS import into GRASS GIS - select region extent only
Wow, great! This trick with a second VRT works with gdal_translate and
r.in.gdal, with the second even in case of regions bigger than MAXSIZE. How did
you find that out, if I may ask you?
What do you think: I have to handle at least 4 raster maps. What is the most
efficient way to organize this bunch of VRTs?
I think, I could base on the function _createXML() from r.in.wms2 an generate
the VRT myself. The second VRT introduces the following important lines
depending on TE parameter:
2.5920e+06, 2.e+00, 0.e+00,
1.1825e+06, 0.e+00,
-2.e+00
Cheers and good night
Martin
--
Martin Zbinden
Riedacker 523
3154 Rüschegg Heubach
+41 78 628 28 82
2014/1/4 Blumentrath, Stefan :
> Actually, a second vrt does the trick:
>
> gdalbuild vrt dem_wcs_bb.vrt dem_wcs.vrt -te 381464.564874
> 7282848.63791 422460.634623 7310018.58218
>
> This worked for me.
> I shall add the solution to the Wiki.
>
> Cheers
> Stefan
>
> -Original Message-
> From: martin.zbin...@gmail.com [mailto:martin.zbin...@gmail.com] On
> Behalf Of Martin Zbinden
> Sent: 4. januar 2014 23:32
> To: Blumentrath, Stefan; nete...@osgeo.org
> Cc: grass-user@lists.osgeo.org
> Subject: Re: [GRASS-user] WCS import into GRASS GIS - select region
> extent only
>
> 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/
___
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] WCS import into GRASS GIS - select region extent only
On 1/4/14, Blumentrath, Stefan wrote: > Actually, a second vrt does the trick: > > gdalbuild vrt dem_wcs_bb.vrt dem_wcs.vrt -te 381464.564874 7282848.63791 > 422460.634623 7310018.58218 > > This worked for me. > I shall add the solution to the Wiki. This is probably better than a -r flag for the raster import. The issue is that if the resolution of the current region does not match that of the raster map or layer to be imported, resampling is needed. But which? It depends on the input data (for available methods in GRASS, see the 'Interpolation' page in the Wiki). Hence, reflecting about it again, such a convenient flag is not easily integrated... Markus ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] WCS import into GRASS GIS - select region extent only
Actually, a second vrt does the trick: gdalbuild vrt dem_wcs_bb.vrt dem_wcs.vrt -te 381464.564874 7282848.63791 422460.634623 7310018.58218 This worked for me. I shall add the solution to the Wiki. Cheers Stefan -Original Message- From: martin.zbin...@gmail.com [mailto:martin.zbin...@gmail.com] On Behalf Of Martin Zbinden Sent: 4. januar 2014 23:32 To: Blumentrath, Stefan; nete...@osgeo.org Cc: grass-user@lists.osgeo.org Subject: Re: [GRASS-user] WCS import into GRASS GIS - select region extent only 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 : > 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-laye > r-with-gdal-translate > > Cheers > Stefan > > -Original Message- > From: grass-user-boun...@lists.osgeo.org > [mailto:grass-user-boun...@lists.osgeo.org] On Behalf Of Blumentrath, > Stefan > Sent: 4. januar 2014 22:17 > To: Martin Zbinden; grass-user@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. > &SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&FORMA > T=GeoTIFF&COVERAGE=layername&BBOX=381464.564874,7282848.63791,7310018. > 58218,422460.634623&CRS=EPSG:32633&WIDTH=800&HEIGHT=600 tra> > > 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.zbin...@gmail.com [mailto:martin.zbin...@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 : >> 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"): >> >> >> http://your_wcs_server/wcs.dtm? >> layername >> >> >> 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),
Re: [GRASS-user] WCS import into GRASS GIS - select region extent only
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 : > 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-boun...@lists.osgeo.org > [mailto:grass-user-boun...@lists.osgeo.org] On Behalf Of Blumentrath, Stefan > Sent: 4. januar 2014 22:17 > To: Martin Zbinden; grass-user@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. > &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 > > 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.zbin...@gmail.com [mailto:martin.zbin...@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 : >> 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"): >> >> >> http://your_wcs_server/wcs.dtm? >> layername >> >> >> 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-boun...@lists.osgeo.org >> [mailto:grass-user-boun...@lists.osgeo.org] On Behalf Of Martin >> Zbinden >> Sent: 3. januar 2014 21:46 >> To: grass-user@lists.osgeo.org >&g
Re: [GRASS-user] WCS import into GRASS GIS - select region extent only
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-boun...@lists.osgeo.org [mailto:grass-user-boun...@lists.osgeo.org] On Behalf Of Blumentrath, Stefan Sent: 4. januar 2014 22:17 To: Martin Zbinden; grass-user@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. &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 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.zbin...@gmail.com [mailto:martin.zbin...@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 : > 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"): > > > http://your_wcs_server/wcs.dtm? > layername > > > 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-boun...@lists.osgeo.org > [mailto:grass-user-boun...@lists.osgeo.org] On Behalf Of Martin > Zbinden > Sent: 3. januar 2014 21:46 > To: grass-user@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
Re: [GRASS-user] WCS import into GRASS GIS - select region extent only
Hi Stefan, On Sat, Jan 4, 2014 at 10:40 PM, Blumentrath, Stefan wrote: > Hi Markus, > > I should not write emails offline ;-), so sorry for reposting parts of your > answer... > And thanks for the hints. I was not aware that the BBOX had no effect in the > XML. ... this is what I understood from FrankW's email. > Hopefully, the HEIGHT and WIDTH parameter work, as Martin seemed to have > problems to fetch the data for his region due to tile-size limits... Please add your findings in this new Wiki page: http://grasswiki.osgeo.org/wiki/WCS > If you consider enhancing r.external (or r.in.gdal) in the direction of a > projwin option, maybe a "limit import to current region flag" would be more > easy to use and maybe also more consistent (like e.g in v.in.ogr). Indeed, this is a long-term wish for me :-) The code in http://svn.osgeo.org/gdal/trunk/gdal/apps/gdal_translate.cpp should give some inspiration... Best Markus ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] WCS import into GRASS GIS - select region extent only
Hi Markus, I should not write emails offline ;-), so sorry for reposting parts of your answer... And thanks for the hints. I was not aware that the BBOX had no effect in the XML. Hopefully, the HEIGHT and WIDTH parameter work, as Martin seemed to have problems to fetch the data for his region due to tile-size limits... If you consider enhancing r.external (or r.in.gdal) in the direction of a projwin option, maybe a "limit import to current region flag" would be more easy to use and maybe also more consistent (like e.g in v.in.ogr). Cheers Stefan -Original Message- From: neteler.os...@gmail.com [mailto:neteler.os...@gmail.com] On Behalf Of Markus Neteler Sent: 4. januar 2014 21:49 To: Blumentrath, Stefan Cc: grass-user@lists.osgeo.org Subject: Re: [GRASS-user] WCS import into GRASS GIS - select region extent only Hi, just for fun I have made a test with r.external (of GRASS GIS 7, it potentially works with G6 as well): # create ASCII file with this content, named e.g. "wcs_geoserver_LL.wcs": http://demo.opengeo.org/geoserver/wcs? Img_Sample or download this file from http://svn.osgeo.org/gdal/trunk/autotest/gdrivers/data/geoserver.wcs Start GRASS in latlong location, then GRASS 7.0.svn (latlong):~ > r.external input=wcs_geoserver_LL.wcs output=wcs_data --o Projection of input dataset and current location appear to match Reading band 1 of 3... r.external complete. Link to raster map created. Reading band 2 of 3... r.external complete. Link to raster map created. Reading band 3 of 3... r.external complete. Link to raster map created. Imagery group created Interestingly the content of wcs_geoserver_LL.wcs is expanded on the fly, i.e. the file is rewritten. BTW I had to run r.external twice on the same file to avoid some GDAL errors (perhaps the demo server is simply stressed). GRASS 7.0.svn (latlong):~ > g.region rast=wcs_data.1 -p projection: 3 (Latitude-Longitude) zone: 0 datum: wgs84 ellipsoid: wgs84 north: 54:06:50.76N south: 20:42:18.72N west: 130:51:06.048W east: 62:00:19.44W nsres: 0:03:21.123813 ewres: 0:04:12.132867 rows: 598 cols: 983 cells: 587834 Voilà, the WCS layer is loaded (then I used the RGB viewer in the wxGUI). Eventually I then played around with the bounding box (BBOX=minx,miny,maxx,maxy) and "injected" it in above small file &bbox=-103.93,33.74-93.4940.61 so that it becomes: http://demo.opengeo.org/geoserver/wcs? Img_Sample &bbox=-103.93,33.74-93.4940.61 BUT, as stated in http://lists.osgeo.org/pipermail/gdal-dev/2008-March/016552.html it does not have any effect. I don't know if it would make sense to implement a "projwin" parameter in r.external. BTW: Another file/WCS server to play with is http://svn.osgeo.org/gdal/trunk/autotest/gdrivers/data/srtmplus.wcs Not really helpful but at least we know that r.external works with WCS servers. Markus ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
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. &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 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.zbin...@gmail.com [mailto:martin.zbin...@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 : > 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"): > > > http://your_wcs_server/wcs.dtm? > layername > > > 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-boun...@lists.osgeo.org > [mailto:grass-user-boun...@lists.osgeo.org] On Behalf Of Martin > Zbinden > Sent: 3. januar 2014 21:46 > To: grass-user@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/G
Re: [GRASS-user] WCS import into GRASS GIS - select region extent only
Hi, just for fun I have made a test with r.external (of GRASS GIS 7, it potentially works with G6 as well): # create ASCII file with this content, named e.g. "wcs_geoserver_LL.wcs": http://demo.opengeo.org/geoserver/wcs? Img_Sample or download this file from http://svn.osgeo.org/gdal/trunk/autotest/gdrivers/data/geoserver.wcs Start GRASS in latlong location, then GRASS 7.0.svn (latlong):~ > r.external input=wcs_geoserver_LL.wcs output=wcs_data --o Projection of input dataset and current location appear to match Reading band 1 of 3... r.external complete. Link to raster map created. Reading band 2 of 3... r.external complete. Link to raster map created. Reading band 3 of 3... r.external complete. Link to raster map created. Imagery group created Interestingly the content of wcs_geoserver_LL.wcs is expanded on the fly, i.e. the file is rewritten. BTW I had to run r.external twice on the same file to avoid some GDAL errors (perhaps the demo server is simply stressed). GRASS 7.0.svn (latlong):~ > g.region rast=wcs_data.1 -p projection: 3 (Latitude-Longitude) zone: 0 datum: wgs84 ellipsoid: wgs84 north: 54:06:50.76N south: 20:42:18.72N west: 130:51:06.048W east: 62:00:19.44W nsres: 0:03:21.123813 ewres: 0:04:12.132867 rows: 598 cols: 983 cells: 587834 Voilà, the WCS layer is loaded (then I used the RGB viewer in the wxGUI). Eventually I then played around with the bounding box (BBOX=minx,miny,maxx,maxy) and "injected" it in above small file &bbox=-103.93,33.74-93.4940.61 so that it becomes: http://demo.opengeo.org/geoserver/wcs? Img_Sample &bbox=-103.93,33.74-93.4940.61 BUT, as stated in http://lists.osgeo.org/pipermail/gdal-dev/2008-March/016552.html it does not have any effect. I don't know if it would make sense to implement a "projwin" parameter in r.external. BTW: Another file/WCS server to play with is http://svn.osgeo.org/gdal/trunk/autotest/gdrivers/data/srtmplus.wcs Not really helpful but at least we know that r.external works with WCS servers. Markus ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] WCS import into GRASS GIS - select region extent only
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"): http://your_wcs_server/wcs.dtm? layername 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-boun...@lists.osgeo.org [mailto:grass-user-boun...@lists.osgeo.org] On Behalf Of Martin Zbinden Sent: 3. januar 2014 21:46 To: grass-user@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&coverage=blw.eros_z&FORMAT=image/tiff&CRS=EPSG:2056&BBOX=2596000,1181000,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_example [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@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
[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&coverage=blw.eros_z&FORMAT=image/tiff&CRS=EPSG:2056&BBOX=2596000,1181000,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_example [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@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
