2

I have many CSV file in the form of lon, lat, Z that in would like to import into gdal for further processing. I need to do this on the fly in a python script so any step involving gdal_grid or anything else from the command line is not an option.

It seemed that creating a vrt xml would solve the problem but I am unable to get gdal to load the data.

Loading data this way seems to be such a common task I can't believe there is not a way to do this.

1
  • I was already able to successfully create a CSV file with the corresponding vrt file that checked out with ogrinfo. However, I was unable to load it with a gdal python script.
    – denson
    Commented Jun 13, 2014 at 18:39

2 Answers 2

1

Step 1: Create CSV file "xyztest.csv"

x,y,z
1,2,3
4,5,6
7,8,9

Step 2: Write a VRT file "xyz.vrt" that maps the previous file

<OGRVRTDataSource>

    <OGRVRTLayer name="xyztest">
        <SrcDataSource>xyztest.csv</SrcDataSource> 
    <SrcLayer>xyztest</SrcLayer> 
    <GeometryType>wkbPoint</GeometryType> 
        <LayerSRS>WGS84</LayerSRS>
    <GeometryField encoding="PointFromColumns" x="x" y="y" z="z"/> 
    </OGRVRTLayer>

</OGRVRTDataSource>

Step 3: Test with ogrinfo

C:\temp>ogrinfo xyz.vrt -al
INFO: Open of `xyz.vrt'
      using driver `VRT' successful.

Layer name: xyztest
Geometry: Point
Feature Count: 3
Extent: (1.000000, 2.000000) - (7.000000, 8.000000)
Layer SRS WKT:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9108"]],
    AUTHORITY["EPSG","4326"]]
x: String (0.0)
y: String (0.0)
z: String (0.0)
OGRFeature(xyztest):1
  x (String) = 1
  y (String) = 2
  z (String) = 3
  POINT (1 2 3)

OGRFeature(xyztest):2
  x (String) = 4
  y (String) = 5
  z (String) = 6
  POINT (4 5 6)

OGRFeature(xyztest):3
  x (String) = 7
  y (String) = 8
  z (String) = 9
  POINT (7 8 9)

All done. It may still be necessary to explicitly define that the input vrt is 3D. For example ogr2ogr conversion into shapefile creates a "POINT" type shapefile unless "-SHPT POINTZ" is used as a parameter.

3
  • Question: How do you actually write the VRT file? In my case, I'm trying to convert a csv of points into a vrt to be used for gdal_grid, all inside a python script. I've looked seemingly everywhere, but all tutorials say is "here is a csv, and here's the resulting vrt," but don't actually say how to get the resulting vrt. Perhaps you can help? Thx
    – user20408
    Commented Dec 11, 2014 at 16:42
  • By hand, if you mean that. I believe I first just read and followed the example at gdal.org/drv_csv.html and once it worked I started to make my own modifications.
    – user30184
    Commented Dec 11, 2014 at 16:58
  • Well, not by hand. It looks like downloading the FW tool kit or whatever is how you obtain the VRT from CSV. Only thing is, I need to automate the conversion to VRT for a large amount of csvs, so doing it by hand is not viable. I tried os.system("ogr2ogr -f 'VRT' %s %s -nlt 'POINT'" %(tempvrt, rbcsv)) but get the error: "VRT driver does not support data source creation." Perhaps it's not possible to do in a script.
    – user20408
    Commented Dec 11, 2014 at 17:03
0

If you have at least GDAL 1.11, you can open CSV files with the XYZ driver. The file must conform to the rules described, such as increasing X values, or the names used for the columns. In lieu of what your file looks like, here's a working example:

$ cat > my.csv
lon, lat, Z
0,0,3
1,0,6
0,1,9
1,1,2

$ gdalinfo my.csv
Driver: XYZ/ASCII Gridded XYZ
Files: my.csv
Size is 2, 2
Coordinate System is `'
Origin = (-0.500000000000000,-0.500000000000000)
Pixel Size = (1.000000000000000,1.000000000000000)
Corner Coordinates:
Upper Left  (  -0.5000000,  -0.5000000)
Lower Left  (  -0.5000000,   1.5000000)
Upper Right (   1.5000000,  -0.5000000)
Lower Right (   1.5000000,   1.5000000)
Center      (   0.5000000,   0.5000000)
Band 1 Block=2x1 Type=Byte, ColorInterp=Undefined
  Min=2.000 Max=9.000

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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