6

I'm trying to find an optimum vector file format to be used to effectively store data (compressed binary?) and select subregions.

I have a quite large ascii point file containing 3 columns (x and y projected coordinates and the point ID). A test file is about 20MB, containing approximately 600000 points, but I plan to create even larger files with significantly more points. After writing those files to ascii, I'm searching for a way (i.e. file format) to store the data more effectively while still allowing for the following two tasks:

  • the files shall be loaded into QGIS (v2.14.17)
  • subregions shall be selected (e.g. via GDAL/OGR v2.1 "ogr2ogr -spat") and rewritten in the original ascii format

Do you have any suggestions, which file format shall be used here?

EDIT: Just to add some more information the data looks like

CoordX,CoordY,ID  
1200400,-399600,1
1201200,-399600,2
1202000,-399600,3
1202800,-399600,4
...

I've created a VRT file as follows:

<OGRVRTDataSource>
    <OGRVRTLayer name="points">
        <SrcDataSource>points.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <LayerSRS>EPSG:3031</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="CoordX" y="CoordY"/>
        <Field name="X" src="CoordX" type="integer" />
        <Field name="Y" src="CoordY" type="integer" />
        <Field name="ID" src="ID" type="integer" />
    </OGRVRTLayer>
</OGRVRTDataSource>

As Luke suggested (What format to use for highly compressed vector data?) I tried SpatiaLite with compressed geometry: ogr2ogr -f SQLite -dsco SPATIALITE=YES -lco COMPRESS_GEOM=YES -lco COMPRESS_COLUMNS=x,y,id points.sqlite points.vrt but the result was a 90MB file (compared to 14MB of the now slightly modified CSV-Input).

Playing around a bit showed, that leaving away the -lco COMPRESS_* did not have any effect on the 90MB file size. Instead, leaving away -dsco SPATIALITE=YES resulted in a 26MB file (I know, this is not SpatiaLite any more). I'm working on Ubuntu 14.04. GDAL 2.1.0 was installed using the UbuntuGIS-PPA. Is there any additional package (besides libspatialite5 v4.1.1) necessary to enable compression in SpatiaLite?

  • Spatialite compressed geometry perhaps. See also gis.stackexchange.com/q/31213/2856 – user2856 Aug 8 '17 at 13:07
  • In your sqlite database, you have stored the points coordinates twice: once in the geometry and once in the x and y fields. You only need them in the geometry field, leaving only one integer column for id. If you extract your data back to CSV, you may use GEOMETRY=AS_XY to recover x and y. – AndreJ Aug 9 '17 at 5:56
  • I agree, but without the two coordinates as attribute columns, the size is still 86MB. Here I just found that "geometry compression only affects LINESTRINGs and POLYGONs, not POINTs", which is a bit more clear as in the OGR Vector Formats page. Looks like SpatiaLite won't help here. – Ludwig Aug 9 '17 at 8:04
  • 1
    Keep in mind that the spatial index you want also needs some space. You can try gdal.org/frmt_netcdf_vector.html and ESRI FileGDB as alternatives. NetCDF vector is supported by GDAL, but probably not yet by QGIS. – AndreJ Aug 9 '17 at 9:47
3

Great AndreJ, that did the trick! I was working with NetCDF for a long time but never thought about using it for vector data.
ogr2ogr -f NetCDF points.nc points.shp produced a NetCDF file which was still 17MB, i.e. more than the input. But this NetCDF was indeed possible to be loaded into QGIS. I didn't find any option to tell ogr2ogr to use compression within the netcdf. However, now that I know the structure how a NetCDF should look to be recognized as vector data by GDAL/OGR, I wrote a Python script to read the data and rewrite it as NetCDF similar to the ogr2ogr-output but including compression. It worked like a charm!! Now the the ~600000 points just occupy 111kB. It needs some time to load in QGIS (I guess it is converted back into a temporary shapefile in the background) but finally it works!

#!/usr/bin/python
import sys
import os
import re
from netCDF4 import Dataset

### read argument 
if(len(sys.argv) < 2 or not re.search("\.csv$",sys.argv[1])):
    print "Usage: ",sys.argv[0],"<input>.csv"
    exit(1)

input_csv=sys.argv[1]
basnam=os.path.splitext(os.path.basename(input_csv))[0]
output_nc="%s.nc"%(basnam)

### read data file
(x,y,id)=([],[],[])
fcsv = open(input_csv,'r')
lines = fcsv.readlines()
for line in lines:
    cols = line.split(',')
    if(cols[0]=="CoordX"): continue
    x.append(int(cols[0]))
    y.append(int(cols[1]))
    id.append(int(cols[2]))
fcsv.close()
print "read",len(id),"records"

### Writing NetCDF file
nc_out = Dataset(output_nc, "w", format="NETCDF4_CLASSIC")
nc_out.Conventions = "CF-1.6"
nc_out.ogr_layer_name = basnam
nc_out.featureType = "point"
## dimension
nc_out.createDimension('record', None)

## variables
# coordinates
nc_out_var = nc_out.createVariable('x','f8',('record'),zlib=True)
nc_out_var.standard_name = "projection_x_coordinate"
nc_out_var.units = "m"
nc_out.variables['x'][:]    = x
nc_out_var = nc_out.createVariable('y','f8',('record'),zlib=True)
nc_out_var.standard_name = "projection_y_coordinate"
nc_out_var.units = "m"
nc_out.variables['y'][:]    = y

# projection
nc_out_var = nc_out.createVariable('polar_stereographic','c')
nc_out_var.grid_mapping_name = "polar_stereographic" ;
nc_out_var.straight_vertical_longitude_from_pole = 0. ;
nc_out_var.false_easting = 0. ;
nc_out_var.false_northing = 0. ;
nc_out_var.latitude_of_projection_origin = -90. ;
nc_out_var.standard_parallel = -71. ;
nc_out_var.long_name = "CRS definition" ;
nc_out_var.longitude_of_prime_meridian = 0. ;
nc_out_var.semi_major_axis = 6378137. ;
nc_out_var.inverse_flattening = 298.257223563 ;
nc_out_var.spatial_ref = "PROJCS[\"WGS_84_Antarctic_Polar_Stereographic\",GEOGCS[\"GCS_WGS_1984\",DATUM[\"WGS_1984\",SPHEROID[\"WGS_84\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]],PROJECTION[\"Polar_Stereographic\"],PARAMETER[\"latitude_of_origin\",-71],PARAMETER[\"central_meridian\",0],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]" ;

# fields
nc_out_var = nc_out.createVariable('CoordX','i4',('record'),zlib=True)
nc_out_var.long_name = "Field X"
nc_out_var.grid_mapping = "polar_stereographic"
nc_out.variables['CoordX'][:]    = x

nc_out_var = nc_out.createVariable('CoordY','i4',('record'),zlib=True)
nc_out_var.long_name = "Field Y"
nc_out_var.grid_mapping = "polar_stereographic"
nc_out.variables['CoordY'][:]    = y

nc_out_var = nc_out.createVariable('ID','i4',('record'),zlib=True)
nc_out_var.long_name = "Field ID"
nc_out_var.grid_mapping = "polar_stereographic"
nc_out.variables['ID'][:]    = id

nc_out.close()

print "Finished writing",output_nc

I consider the topic as solved, however an option to create a compressed NetCDF using ogr2ogr would be desirable.

  • Nice work, your solution sounds promising. It would be great and useful for others if you could include your Python script in your answer. – Stefan Aug 9 '17 at 12:17
  • NetCDF is not in the filetype filter of ´Add vector Layer", but it seems to work with All files. It does NOT load from the QGIS browser, because it is classified as raster there. – AndreJ Aug 10 '17 at 15:31
  • 1
    Nevertheless, it exists in the OGR Vector Formats for GDAL>=2.1, hence the GDAL/OGR can import NetCDF as vector layer too. – Ludwig Aug 10 '17 at 19:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.