[GRASS-user] How to best access hundreds of orthophotos -- mapserver as WMS server?

GRASS GRASS user
1

Hi

I have hundreds, probably thousands, of orthophotos which I would like to access as easy as possible. They are in MrSid format, but I will convert them to geotiff.

My problem is that it is sometimes quite difficult to find the one covering the region I am working in and it involves quite some trial and error to get the right one (which is time intensive and annoying...).

I thought about setting up mapserver as a WMS server, but I didn't manage (might have to spend some m ore time on that). What I would like to know:

1) Is mapserver as WMS server the best solution or is there a more effective one?

2) If mapserver (WMS server) is the way to go - can somebody give me an example .map file which is achieving that?

3) the whole layer is obviously in several single tiles - can mapserver handle this as a single layer or do I( have to do something special?

Sorry for the slightly oftopic questions - but 1) should be ontopic.

Thanks

Rainer

2

Hi Rainer,
Rainer M Krug píše v Pá 14. 12. 2007 v 09:22 +0200:

Hi

I have hundreds, probably thousands, of orthophotos which I would like
to access as easy as possible. They are in MrSid format, but I will
convert them to geotiff.

My problem is that it is sometimes quite difficult to find the one
covering the region I am working in and it involves quite some trial and
error to get the right one (which is time intensive and annoying...).

I thought about setting up mapserver as a WMS server, but I didn't
manage (might have to spend some m ore time on that). What I would like
to know:

1) Is mapserver as WMS server the best solution or is there a more
effective one?

IMHO yes - but I would advice WCS to you. WMS does render "map", where
WCS returns raw data (coverage) to you.

2) If mapserver (WMS server) is the way to go - can somebody give me an
example .map file which is achieving that?

Some example attached - but please, take it as skeleton. I do not
expect, single line would stay untouched after you start to modify it.

have a look at
http://mapserver.gis.umn.edu/docs/reference/mapfile
http://mapserver.gis.umn.edu/docs/howto/wcs_server
http://mapserver.gis.umn.edu/docs/howto/tileindex

3) the whole layer is obviously in several single tiles - can mapserver
handle this as a single layer or do I( have to do something special?

see the tileindex link above

Sorry for the slightly oftopic questions - but 1) should be ontopic.

Thanks

Rainer

good luck :slight_smile:

j
--
Jachym Cepicky
e-mail: jachym.cepicky@gmail.com
URL: http://les-ejk.cz
GPG: http://www.les-ejk.cz/pgp/jachym_cepicky-gpg.pub

(attachments)

mapfile.map (2.47 KB)

3

Jachym Cepicky wrote:

Hi Rainer,
Rainer M Krug píše v Pá 14. 12. 2007 v 09:22 +0200:

Hi

I have hundreds, probably thousands, of orthophotos which I would like to access as easy as possible. They are in MrSid format, but I will convert them to geotiff.

My problem is that it is sometimes quite difficult to find the one covering the region I am working in and it involves quite some trial and error to get the right one (which is time intensive and annoying...).

I thought about setting up mapserver as a WMS server, but I didn't manage (might have to spend some m ore time on that). What I would like to know:

1) Is mapserver as WMS server the best solution or is there a more effective one?

IMHO yes - but I would advice WCS to you. WMS does render "map", where
WCS returns raw data (coverage) to you.

thanks for pointing that out - I will consider it. But can I access the WCS server from GRASS and QGIS as well?

2) If mapserver (WMS server) is the way to go - can somebody give me an example .map file which is achieving that?

Some example attached - but please, take it as skeleton. I do not
expect, single line would stay untouched after you start to modify it.

I'll do my best - but is always eaier to understand what's going on when having a template to start with.

have a look at
http://mapserver.gis.umn.edu/docs/reference/mapfile
http://mapserver.gis.umn.edu/docs/howto/wcs_server
http://mapserver.gis.umn.edu/docs/howto/tileindex

Thanks - I'll look at them as well.

3) the whole layer is obviously in several single tiles - can mapserver handle this as a single layer or do I( have to do something special?

see the tileindex link above

Sorry for the slightly oftopic questions - but 1) should be ontopic.

Thanks

Rainer

good luck :slight_smile:

That sounds as if I would need it ...

j

------------------------------------------------------------------------

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

4

> Rainer Krug writes:
>> I have hundreds, probably thousands, of orthophotos which I would
>> like to access as easy as possible.

...

>> My problem is that it is sometimes quite difficult to find the one
>> covering the region I am working in and it involves quite some
>> trial and error to get the right one (which is time intensive and
>> annoying...).

an idea: write a little script to make a vector coversheet index.
maybe output that with the HTMLMAP driver for an interactive web brower
pull up.

I am not sure, but I think v.patch does not clean topology so overlap
may be ok. Also I am not sure about using v.edit to reasign a category
number. Also the "for MAP in ``" method would need to be changed if
there are thousands of maps?

this is just the idea. non functional...

# store current region
g.region save=old_region
# create new empty map
v.in.ascii -e out=map_index
v.db.addtable map_index columns="image varchar(128)"

i=0
for MAP in `g.mlist patt='ortho_*'` ; do
  i=`expr $i + 1`
  g.region rast="$MAP"
  v.in.region out="${MAP}_box"

? v.edit to reset the region box centroid's cat to $i ?
  # or a series of v.category del then add steps

  # merge into master index
  g.remove vect=map_index_old
  g.rename vect=map_index,map_index_old
  v.patch in=map_index_old,"${MAP}_box" out=map_index

  # write map name to DB entry
  v.db.update map_index column=image value="'$MAP'" where="cat = $i"
done

#reset to original region
g.region old_region

v.label map=map_index column=image
d.labels
  #or
d.vect display=attr attrcol=image xref=center

Hamish

      ____________________________________________________________________________________
Be a better friend, newshound, and
know-it-all with Yahoo! Mobile. Try it now. http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ

5

Hamish <hamish_b@yahoo.com> writes:

[...]

> Also the "for MAP in ``" method would need to be changed if there are
> thousands of maps?

  I don't think so. Since the `for' command is internal to the
  Shell, there's no command line to be passed to any other
  process, and so there're no command line length limits in force.

  Still, using COMMAND | while read VARN ; do ...; done looks more
  clean to me.

> this is just the idea. non functional...

[...]

> i=0
> for MAP in `g.mlist patt='ortho_*'` ; do

[...]

> done

[...]

6

Hamish <hamish_b@yahoo.com> writes:

[...]

> an idea: write a little script to make a vector coversheet index.
> maybe output that with the HTMLMAP driver for an interactive web brower
> pull up.

> I am not sure, but I think v.patch does not clean topology so overlap
> may be ok.

  Wouldn't it be better to save non-overlapping parts and an
  intersection as separate polygons? A separate attribute will be
  needed to store the list of IDs of the covering images.

[...]

> this is just the idea. non functional...

> # store current region
> g.region save=old_region

  It makes me wonder each time I see such a fragile construct like
  this, will GRASS ever support overriding the region
  ``temporarily''? (E. g., via a command line argument, or an
  environment variable.)

> # create new empty map
> v.in.ascii -e out=map_index
> v.db.addtable map_index columns="image varchar(128)"

> i=0
> for MAP in `g.mlist patt='ortho_*'` ; do
> i=`expr $i + 1`
> g.region rast="$MAP"
> v.in.region out="${MAP}_box"

  Oh, perhaps `v.in.region' could be patched to print the category
  of the vector feature it creates (if requested)? I guess, this
  usage pattern is quite common for `v.in.region'?

  And, could there be something a bit more clever than storing the
  region of the image as a vector feature? A combination of
  downsampling, `r.mapcalc' (to obtain the fill mask), and
  `r.contour', may be?

> ? v.edit to reset the region box centroid's cat to $i ?
> # or a series of v.category del then add steps

> # merge into master index
> g.remove vect=map_index_old
> g.rename vect=map_index,map_index_old
> v.patch in=map_index_old,"${MAP}_box" out=map_index

  May there `v.patch -a' be used instead? Like:

v.patch -a in="${MAP}_box" output=map_index

> # write map name to DB entry
> v.db.update map_index column=image value="'$MAP'" where="cat = $i"
> done

> #reset to original region
> g.region old_region

[...]

7

On Sunday 16 December 2007 04:12:09 am Ivan Shmakov wrote:

>>>>> Hamish <hamish_b@yahoo.com> writes:

[...]

> an idea: write a little script to make a vector coversheet index.
> maybe output that with the HTMLMAP driver for an interactive web brower
> pull up.
>
> I am not sure, but I think v.patch does not clean topology so overlap
> may be ok.

  Wouldn't it be better to save non-overlapping parts and an
  intersection as separate polygons? A separate attribute will be
  needed to store the list of IDs of the covering images.

[...]

> this is just the idea. non functional...
>
> # store current region
> g.region save=old_region

  It makes me wonder each time I see such a fragile construct like
  this, will GRASS ever support overriding the region
  ``temporarily''? (E. g., via a command line argument, or an
  environment variable.)

> # create new empty map
> v.in.ascii -e out=map_index
> v.db.addtable map_index columns="image varchar(128)"
>
> i=0
> for MAP in `g.mlist patt='ortho_*'` ; do
> i=`expr $i + 1`
> g.region rast="$MAP"
> v.in.region out="${MAP}_box"

  Oh, perhaps `v.in.region' could be patched to print the category
  of the vector feature it creates (if requested)? I guess, this
  usage pattern is quite common for `v.in.region'?

  And, could there be something a bit more clever than storing the
  region of the image as a vector feature? A combination of
  downsampling, `r.mapcalc' (to obtain the fill mask), and
  `r.contour', may be?

Hi,

What you are asking for is starting to sound like a Mapserver driven WMS.
Here is how I have dealt with this problem in the past:

* get and compile the latest version of GDAL
* get and compile the latest version of GRASS
* get and compile the latest version of Mapserver

1. put all of your imagery somewhere, in a consistant format.
2. add overviews to the imagery (internal tiling) with gdaladdo
3. make a shapefile tile index with gdaltindex
4. make a spatial index on this shapefile with shptree (from the mapserver
source code)
5. create a very simple mapfile pointing to the tile index shapefile. I will
work up an example and post it back
6. fill in appropriate WMS metadata into the mapfile

7. use r.in.wms to import imagery into GRASS
-or-
7. once GDAL 1.5 is released it will be possible to an r.in.gdal and point it
to an XML file defining a WMS resource. A mini-script might be able to
construct the XML file on the fly...

Unfortunately you will still have to manage the images that are imported into
GRASS with GRASS commands, unless you constantly overwrite a temporary set of
files by always using the same names and the --o flag.

Once you have mapserver setup WMS is quite fast and flexible - it will even
project the data on the fly to other coordinate systems.

Dylan

> ? v.edit to reset the region box centroid's cat to $i ?
> # or a series of v.category del then add steps
>
> # merge into master index
> g.remove vect=map_index_old
> g.rename vect=map_index,map_index_old
> v.patch in=map_index_old,"${MAP}_box" out=map_index

  May there `v.patch -a' be used instead? Like:

v.patch -a in="${MAP}_box" output=map_index

> # write map name to DB entry
> v.db.update map_index column=image value="'$MAP'" where="cat = $i"
> done
>
> #reset to original region
> g.region old_region

[...]

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

8

Hamish:

> an idea: write a little script to make a vector coversheet index.
> maybe output that with the HTMLMAP driver for an interactive web
> browser pull up.
>
> I am not sure, but I think v.patch does not clean topology so
> overlap may be ok.

Ivan:

Wouldn't it be better to save non-overlapping parts and an
intersection as separate polygons? A separate attribute will be
needed to store the list of IDs of the covering images.

see overlapping polygon example at:
  http://grass.ibiblio.org/screenshots/vector.php

> # store current region
> g.region save=old_region

  It makes me wonder each time I see such a fragile construct like
  this, will GRASS ever support overriding the region
  ``temporarily''? (E. g., via a command line argument, or an
  environment variable.)

there is, see GRASS_REGION and WIND_OVERRIDE shell variables:
  http://grass.ibiblio.org/grass63/manuals/html63_user/variables.html

but saving & restoring g.region isn't need for this task, it is just
being nice leaving the place as you found it. The script is not
dependent on the starting region. And I fail to see how that's fragile?
No more fragile than trusting anything else written to the disk.

> v.in.region out="${MAP}_box"

Oh, perhaps `v.in.region' could be patched to print the category
of the vector feature it creates (if requested)?

ie have a new cat= option added to set the cat number. (default "1")
right now it always creates as cat 1.

I guess, this usage pattern is quite common for `v.in.region'?

no idea. I've used the module a lot and this is first time I've had to
worry about that.

And, could there be something a bit more clever than storing the
region of the image as a vector feature?

the user can get as complicated as they like..

A combination of downsampling, `r.mapcalc' (to obtain the fill mask),
and `r.contour', may be?

r.contour is not so good for ortho photos or high resolution land use
map with many small areas. But again, up to the user....

> v.patch in=map_index_old,"${MAP}_box" out=map_index

May there `v.patch -a' be used instead? Like:

v.patch -a in="${MAP}_box" output=map_index

yes, that's the correct way.

Hamish

      ____________________________________________________________________________________
Never miss a thing. Make Yahoo your home page.
http://www.yahoo.com/r/hs

9

Hamish <hamish_b@yahoo.com> writes:

>>> an idea: write a little script to make a vector coversheet index.

[...]

>>> I am not sure, but I think v.patch does not clean topology so
>>> overlap may be ok.

> Ivan:

>> Wouldn't it be better to save non-overlapping parts and an
>> intersection as separate polygons? A separate attribute will be
>> needed to store the list of IDs of the covering images.

> see overlapping polygon example at:
> http://grass.ibiblio.org/screenshots/vector.php

  Thanks, I'll check it.

>>> # store current region
>>> g.region save=old_region

>> It makes me wonder each time I see such a fragile construct like
>> this, will GRASS ever support overriding the region
>> ``temporarily''? (E. g., via a command line argument, or an
>> environment variable.)

> there is, see GRASS_REGION and WIND_OVERRIDE shell variables:
> http://grass.ibiblio.org/grass63/manuals/html63_user/variables.html

  This is exactly what I've asked for. Thanks!

> but saving & restoring g.region isn't need for this task, it is just
> being nice leaving the place as you found it. The script is not
> dependent on the starting region. And I fail to see how that's
> fragile?

  Consider this script being interrupted in the middle (say, by
  the user hitting ^C.) Will the region be restored?

> No more fragile than trusting anything else written to the disk.

  It looks to me that the following will be much more robust:

export WIND_OVERRIDE=tmp_region
...

# no need to store current region
# g.region save=old_region

...

  g.region rast="$MAP"

...

# no need to reset to original region, either
# g.region old_region

[...]

>> And, could there be something a bit more clever than storing the
>> region of the image as a vector feature?

> the user can get as complicated as they like..

>> A combination of downsampling, `r.mapcalc' (to obtain the fill mask),
>> and `r.contour', may be?

> r.contour is not so good for ortho photos or high resolution land use
> map with many small areas. But again, up to the user....

  Actually, I meant using `r.contour' on the downsampled fill
  mask. (Meanwhile, I've silently switched to my own needs, which
  are about the satellite data.)

[...]

10

Dylan Beaudette <dylan.beaudette@gmail.com> writes:

[...]

>>> an idea: write a little script to make a vector coversheet index.
>>> maybe output that with the HTMLMAP driver for an interactive web brower
>>> pull up.

[...]

>> And, could there be something a bit more clever than storing the
>> region of the image as a vector feature? A combination of
>> downsampling, `r.mapcalc' (to obtain the fill mask), and
>> `r.contour', may be?

> Hi,

> What you are asking for is starting to sound like a Mapserver driven
> WMS.

  Actually, I'm seeking for a some kind of cataloguing solution
  for GRASS. However, setting up Mapserver is also on my TODO
  list.

> Here is how I have dealt with this problem in the past:

[...]

> 5. create a very simple mapfile pointing to the tile index
> shapefile. I will work up an example and post it back

  TIA.

> 6. fill in appropriate WMS metadata into the mapfile

> 7. use r.in.wms to import imagery into GRASS
> -or-
> 7. once GDAL 1.5 is released it will be possible to an r.in.gdal and
> point it to an XML file defining a WMS resource. A mini-script might
> be able to construct the XML file on the fly...

> Unfortunately you will still have to manage the images that are
> imported into GRASS with GRASS commands, unless you constantly
> overwrite a temporary set of files by always using the same names and
> the --o flag.

> Once you have mapserver setup WMS is quite fast and flexible - it
> will even project the data on the fly to other coordinate systems.

  Yes, I know. (Though I hadn't succeed in getting it to show me
  any rasters last time I've tried it; the software packages were
  from Debian etch, if that matters.)

11

Ivan Shmakov wrote:

> Also the "for MAP in ``" method would need to be changed if there are
> thousands of maps?

  I don't think so. Since the `for' command is internal to the
  Shell, there's no command line to be passed to any other
  process, and so there're no command line length limits in force.

  Still, using COMMAND | while read VARN ; do ...; done looks more
  clean to me.

The while/read approach saves keeping the list of maps in memory. Once
upon a time, that might have mattered; with modern systems, you would
need millions of maps before it became an issue.

--
Glynn Clements <glynn@gclements.plus.com>

12

Glynn Clements <glynn@gclements.plus.com> writes:

>>> Also the "for MAP in ``" method would need to be changed if there
>>> are thousands of maps?

>> I don't think so. Since the `for' command is internal to the Shell,
>> there's no command line to be passed to any other process, and so
>> there're no command line length limits in force.

>> Still, using COMMAND | while read VARN ; do ...; done looks more
>> clean to me.

> The while/read approach saves keeping the list of maps in
> memory. Once upon a time, that might have mattered; with modern
> systems, you would need millions of maps before it became an issue.

  ``while read'' saves one from using one level of backquotes,
  thus, in my opinion, making the code a bit more readable.

  Consider:

for r in `g.mlist pattern="example-\`date +%Y-%m-%d\`-*"`; do
    ...
done

g.mlist pattern="example-`date +%Y-%m-%d`-*" | while read r; do
    ...
done

  Besides, the body of a ``while read'' loop is executed ``in
  parallel'' to the command to the left of the pipe. On the
  contrary, the ``for'' pattern will wait until the command inside
  of the backquotes completes.

  There're whitespace issues, though hardly relevant to the
  `g.mlist' case.