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:low0 0</gml:low>
gml:high39 29</gml:high>
</gml:GridEnvelope>
</gml:limits>
gml:axisLabelslat long</gml:axisLabels>
<gml:Point gml:id=“grid_origin_grey”
srsName=“http://www.opengis.net/def/crs/EPSG/0/4326”>
gml:pos0.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) <jukka.rahkonen@…6847…> 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:
DescribeCoverage from this layer is attached. There are some odd things in the document:
- The Envelope in boundedBy has coordinates in wrong order N-E though they should be E-N
<gml:Envelope srsName=“http://www.opengis.net/def/crs/EPSG/0/3067” axisLabels=“E N” uomLabels=“m m” srsDimension=“2”>
gml:lowerCorner6750000.0 404000.0</gml:lowerCorner>
gml:upperCorner6762000.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:lowerCorner404000.0 6750000.0</gml:lowerCorner>
gml:upperCorner416000.0 6762000.0</gml:upperCorner>
</gml:Envelope>
</gml:boundedBy>
Full describe results attached.
- 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:sequenceRule axisOrder=“+2 +1”>Linear</gml:sequenceRule>
gml:startPoint0 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?
- 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:Point gml:id=“p00_topp__orto_nls” srsName=“http://www.opengis.net/def/crs/EPSG/0/3067”>
gml:pos6761999.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:
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
fax: +39 0584 1660272
mob: +39 339 8844549
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.