Hi Andrea,

Great news. I took some measures if you want to compare: GDAL 2.1-dev with 
OpenJPEG driver is using 1.4 GB of memory and the conversion of one jp2 into 
tiled tiff takes 16 seconds. You should have no need to accept anything worse. 
Or did you run commands in parallel? The input JPEG2000 images are INSPIRE 
compliant and OpenJPEG handles them quite well. Perhaps it would be a time to 
add GDAL OpenJPEG driver to the supported GDAL image formats 
http://docs.geoserver.org/stable/en/user/data/raster/gdal.html#data-gdal as an 
open source alternative for JP2ECW and JP2KAK?

Reasonable tiling structure would be nice to have for large output files. I 
would immediately convert the output tiff which has this wide stripes “Band 1 
Block=24000x8 Type=Byte, ColorInterp=Red” into 256x256 tiled GeoTIFF before 
using it for anything serious. However, I guess that WCS users would like to 
get an output that is immediately usable without further processing.

Tiling is also defined in the document “OGC® Web Coverage Service 2.0 Interface 
Standard - GeoTIFF Coverage Encoding Extension” 
https://portal.opengeospatial.org/files/?artifact_id=51142
Mapserver seems to implement this extension which brings support for these 
selections:
    GEOTIFF:COMPRESSION=compression: The compression method used for the 
returned image. Valid options are: None, PackBits, Deflate, Huffman, LZW and 
JPEG.
    GEOTIFF:JPEG_QUALITY=1-100: When the compression method JPEG is chosen, 
this value defines the quality of the algorithm.
    GEOTIFF:PREDICTOR=None|Horizontal|FloatingPoint: The predictor value for 
the LZW and Deflate compression methods.
    GEOTIFF:INTERLEAVE=Band|Pixel: Defines wether the image shall be band or 
pixel interleaved.
    GEOTIFF:TILING=true|false: Defines if the output image shall be internally 
tiled.
    GEOTIFF:TILEWIDTH=tilewidth, GEOTIFF:TILEHEIGHT=tileheight: Define the size 
of the internal tiles. Must be positive integer divisible by 16.

I have been wondering what would be some practical maximum file size that WCS 
service should be able to deliver. 1-2 GB of high resolution imagery does not 
cover especially large area so probably the service should be made to survive 
such requests. However, above that the processing and data transfer might take 
too long time and lead to timeouts, truncated files and other problems. I 
wonder also how the servers would handle many concurrent large requests.  But 
coverages like nationwide DEMs can be very large and WCS does not directly 
support paging yet so I do not know how a user could collect the whole coverage 
from WCS.

Mapserver and Rasdaman define the RectifiedGrid for EPSG:4326 in the same way 
than in the examples of the document “OGC® GML Application Schema - Coverages - 
GeoTIFF Coverage Encoding Profile”
https://portal.opengeospatial.org/files/?artifact_id=54813

I paste the example below. In short, it defines that firs axis is “lat” and the 
second is “long”. First offset vector is increasing values along the second 
axis (long) and the second offset vector is decreasing values along the firs 
axis (lat).

<gml:RectifiedGrid dimension="2" gml:id="grid_grey">
<gml:limits>
<gml:GridEnvelope>
<gml:low>0 0</gml:low>
<gml:high>39 29</gml:high>
</gml:GridEnvelope>
</gml:limits>
<gml:axisLabels>lat long</gml:axisLabels>
<gml:origin>
<gml:Point gml:id="grid_origin_grey"
srsName="http://www.opengis.net/def/crs/EPSG/0/4326";>
<gml:pos>0.0030991877286375 0.0009432310483165</gml:pos>
</gml:Point>
</gml:origin>
<gml:offsetVector srsName="http://www.opengis.net/def/crs/EPSG/
0/4326">0 0.000089831528393</gml:offsetVector>
<gml:offsetVector srsName="http://www.opengis.net/def/crs/EPSG/
0/4326">-0.000089831528393 0</gml:offsetVector>
</gml:RectifiedGrid>

I will study more the pixel is center/pixel is area thing. First tests that I 
made seem to prove that Mapserver 7.0 and Geoserver 2.7 are producing 
pixel-by-pixel identical results with the same WCS 2.0 GetCoverage from the 
same source image and with the same area subset which is great. Probably it 
means both developer teams have made a similar interpretation about the pixel 
is point/area and how subsetting should be performed.

-Jukka Rahkonen-



Andrea Aime wrote:


On Fri, Sep 25, 2015 at 1:39 PM, Rahkonen Jukka (MML) 
<[email protected]<mailto:[email protected]>>
 wrote:

Hi,



Here is a procedure for creating the issue:



Download 4 ortophotos:

http://kartat.kapsi.fi/files/orto/etrs-tm35fin/mavi_v_25000_50/2015/L44/02m/1/L4422A.jp2

http://kartat.kapsi.fi/files/orto/etrs-tm35fin/mavi_v_25000_50/2015/L44/02m/1/L4422B.jp2

http://kartat.kapsi.fi/files/orto/etrs-tm35fin/mavi_v_25000_50/2015/L44/02m/1/L4422C.jp2

http://kartat.kapsi.fi/files/orto/etrs-tm35fin/mavi_v_25000_50/2015/L44/02m/1/L4422D.jp2



Convert JPEG2000 files into GeoTIFFs
gdal_translate -of GTiff -co tiled=yes L4422A.jp2 L4422A.tif
gdal_translate -of GTiff -co tiled=yes L4422B.jp2 L4422B.tif
gdal_translate -of GTiff -co tiled=yes L4422C.jp2 L4422C.tif
gdal_translate -of GTiff -co tiled=yes L4422D.jp2 L4422D.tif
Wow, they really put my machine to the test, gdal was using 6GB+ of memory 
during
the translation (no kakadu).



Create store and publish the layer and send GetCoverage:

http://localhost:8080/geoserver/wcs?service=WCS&version=2.0.1&request=describecoverage&coverageid=topp__orto_nls



DescribeCoverage from this layer is attached. There are some odd things in the 
document:



1) The Envelope in boundedBy has coordinates in wrong order N-E though they 
should be E-N



<gml:boundedBy>

<gml:Envelope srsName="http://www.opengis.net/def/crs/EPSG/0/3067"; 
axisLabels="E N" uomLabels="m m" srsDimension="2">

<gml:lowerCorner>6750000.0 404000.0</gml:lowerCorner>

<gml:upperCorner>6762000.0 416000.0</gml:upperCorner>

</gml:Envelope>

</gml:boundedBy>

Fixed as part of https://osgeo-org.atlassian.net/browse/GEOS-7142, I now get:

gml:Envelope srsName="http://www.opengis.net/def/crs/EPSG/0/3067"; axisLabels="E 
N" uomLabels="m m" srsDimension="2">
<gml:lowerCorner>404000.0 6750000.0</gml:lowerCorner>
<gml:upperCorner>416000.0 6762000.0</gml:upperCorner>
</gml:Envelope>
</gml:boundedBy>

Full describe results attached.



2) If may be correct to use the axisOrder rule +2 +1 but other WCS servers 
which I have tried do not use it (MapServer, Rasdaman)



<gml:GridFunction>

<gml:sequenceRule axisOrder="+2 +1">Linear</gml:sequenceRule>

<gml:startPoint>0 0</gml:startPoint>

</gml:GridFunction>



This is used to indicate the grid pixels have one orientation, but the axis 
order is the opposite.
What do MapServer/Rasdaman do to describe a coverage in 
http://www.opengis.net/def/crs/EPSG/0/4326, which has flipped axis in WCS 2.0?


3) It is too difficult for me to follow it the mentioned GridFunction together 
with these origin and offsetVectors lead to a correct result or not



<gml:origin>

<gml:Point gml:id="p00_topp__orto_nls" 
srsName="http://www.opengis.net/def/crs/EPSG/0/3067";>

<gml:pos>6761999.75 404000.25</gml:pos>

</gml:Point>

</gml:origin>

<gml:offsetVector srsName="http://www.opengis.net/def/crs/EPSG/0/3067";>0.0 
0.5</gml:offsetVector>

<gml:offsetVector srsName="http://www.opengis.net/def/crs/EPSG/0/3067";>-0.5 
0.0</gml:offsetVector>



So, first vector is increasing 404000.25 (which is Easting) in 0.5 m steps, the 
other is decreasing 6761999.75 (which is Northing in 0.5 m steps and at some 
moment the whole coverage is turned around with the GridFunction.



I attach also WCS 2.0.1 GetCoverage from Mapserver 7.0 so it is possible to 
compare how it is expressing boundedBy, origin, and offSet vectors. At least 
for me this GML without GridFunction is easier to understand. Rasdaman is 
creating GML pretty much alike MapServer.



There seems to be once again half a pixel shifts in the origins from Geoserver 
vs. Mapserver.

GeoServer is as far as I know trying to repect the OGC rule of "pixel is 
center", which is not how raster data is saved
unfortunately (normally raster data is encoded as pixel is corner instead).

I've then tried running:

http://localhost:8080/geoserver/wcs?service=WCS&version=2.0.0&request=GetCoverage&coverageid=finmosaic

and now GeoServer is busy writing the result to a temp file (we cannot stream 
it directly), at least by the looks
of the following trace (which I got using jstack, it's not a recorded error, 
it's a way to check what a Java app
is busy working on):

"btpool0-2" prio=10 tid=0x00007f0930006800 nid=0x1c6a runnable 
[0x00007f09743ab000]
   java.lang.Thread.State: RUNNABLE
        at java.io.RandomAccessFile.writeBytes0(Native Method)
        at java.io.RandomAccessFile.writeBytes(RandomAccessFile.java:520)
        at java.io.RandomAccessFile.write(RandomAccessFile.java:550)
        at 
javax.imageio.stream.FileCacheImageOutputStream.write(FileCacheImageOutputStream.java:141)
        at 
it.geosolutions.imageioimpl.plugins.tiff.TIFFNullCompressor.encode(TIFFNullCompressor.java:104)
        at 
it.geosolutions.imageioimpl.plugins.tiff.TIFFImageWriter.writeTile(TIFFImageWriter.java:1905)
        at 
it.geosolutions.imageioimpl.plugins.tiff.TIFFImageWriter.write(TIFFImageWriter.java:2920)
        at 
it.geosolutions.imageioimpl.plugins.tiff.TIFFImageWriter.write(TIFFImageWriter.java:2614)
        at 
org.geotools.gce.geotiff.GeoTiffWriter.writeImage(GeoTiffWriter.java:427)
        at org.geotools.gce.geotiff.GeoTiffWriter.write(GeoTiffWriter.java:271)

This is going to take a while, the mosaic is 24k by 24k... looking at it, we 
probably need some tweaks in the
code to generate such a large output, I can see there are extra operations in 
the chain that could be
avoided for this case, and we could force a "reasonable" tiling structure in 
the output too, now it's
following the input one,  which is... really bad.

Here is a prelim gdalinfo against the temp file, to give you an idea:


gdalinfo /tmp/imageio3762159434441867746.tmp
Driver: GTiff/GeoTIFF
Files: /tmp/imageio3762159434441867746.tmp
Size is 24000, 24000
Coordinate System is:
PROJCS["ETRS89 / TM35FIN(E,N)",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4258"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",27],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","3067"]]
Origin = (404000.000000000000000,6762000.000000000000000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  404000.000, 6762000.000) ( 25d13'33.38"E, 60d58'52.54"N)
Lower Left  (  404000.000, 6750000.000) ( 25d13'54.90"E, 60d52'24.88"N)
Upper Right (  416000.000, 6762000.000) ( 25d26'51.34"E, 60d59' 2.39"N)
Lower Right (  416000.000, 6750000.000) ( 25d27'10.18"E, 60d52'34.68"N)
Center      (  410000.000, 6756000.000) ( 25d20'22.45"E, 60d55'43.79"N)
Band 1 Block=24000x8 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=24000x8 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=24000x8 Type=Byte, ColorInterp=Blue
  NoData Value=0


Cheers
Andrea

--
==
GeoServer Professional Services from the experts! Visit
http://goo.gl/it488V for more information.
==

Ing. Andrea Aime
@geowolf
Technical Lead

GeoSolutions S.A.S.
Via Poggio alle Viti 1187
55054  Massarosa (LU)
Italy
phone: +39 0584 962313<tel:%2B39%200584%20962313>
fax: +39 0584 1660272<tel:%2B39%200584%201660272>
mob: +39  339 8844549<tel:%2B39%20%C2%A0339%208844549>

http://www.geo-solutions.it
http://twitter.com/geosolutions_it


AVVERTENZE AI SENSI DEL D.Lgs. 196/2003

Le informazioni contenute in questo messaggio di posta elettronica e/o nel/i 
file/s allegato/i sono da considerarsi strettamente riservate. Il loro utilizzo 
è consentito esclusivamente al destinatario del messaggio, per le finalità 
indicate nel messaggio stesso. Qualora riceviate questo messaggio senza esserne 
il destinatario, Vi preghiamo cortesemente di darcene notizia via e-mail e di 
procedere alla distruzione del messaggio stesso, cancellandolo dal Vostro 
sistema. Conservare il messaggio stesso, divulgarlo anche in parte, 
distribuirlo ad altri soggetti, copiarlo, od utilizzarlo per finalità diverse, 
costituisce comportamento contrario ai principi dettati dal D.Lgs. 196/2003.



The information in this message and/or attachments, is intended solely for the 
attention and use of the named addressee(s) and may be confidential or 
proprietary in nature or covered by the provisions of privacy act (Legislative 
Decree June, 30 2003, no.196 - Italy's New Data Protection Code).Any use not in 
accord with its purpose, any disclosure, reproduction, copying, distribution, 
or either dissemination, either whole or partial, is strictly forbidden except 
previous formal approval of the named addressee(s). If you are not the intended 
recipient, please contact immediately the sender by telephone, fax or e-mail 
and delete the information in this message that has been received in error. The 
sender does not give any warranty or accept liability as the content, accuracy 
or completeness of sent messages and accepts no responsibility  for changes 
made after they were sent or for other risks which arise as a result of e-mail 
transmission, viruses, etc.

-------------------------------------------------------
------------------------------------------------------------------------------
_______________________________________________
Geoserver-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/geoserver-users

Reply via email to