1

I have a shapefile, made up of many individual polygons that are essentially arranged in a pixel structure (this file: https://data.gov.uk/dataset/vessel-density-grid-2014). I am interested in the Annual_Ave field data.

I also have a regular array of points that I wish to interpolate against this data. The problem that I have is that the vector data is at a different angle to my array of points (my points are parallel to lines of constant latitude and longitude). The first image roughly illustrates what I mean, the second image is a close up of the actual data to illustrate the problem. I want to be able to perform this interpolation in Python, using GDAL/OGR, and am happy to store the vector data in a different format to facilitate this.

I tried converting the vector polygon layer to a raster, using the ArcGIS polygon to raster function, in the hope that I could extract the latitude and longitude points of each pixel and use this to interpolate my own data. However, when I import the resulting raster into Python using GDAL, the rotational information does not seem to be present so the geotransform is incorrect.

Example (raster file created with ArcGIS Polygon to raster, Cell Size=0.004):

import gdal
in_ds = gdal.Open(<path to raster file>)
in_ds.GetGeoTransform()
Out[39]: 
(-10.036726874521381, 0.004, 0.0, 63.932379486900615, 0.0, 
-0.003999999999999999)

As you can see, the rotation field is 0 even though the polygons are clearly rotated as can be seen in the file. This means that the latitude and longitude values that can be generated from the geotransform data do not match up with those for the pixels in the file.

The other thing that I considered was loading all of the polygons into Python, getting the latitude and longitudes of all the points from the polygon centres, ordering them and interpolating my data points against this. The issue I think with this is that it will be pretty slow, as it is a large vector dataset containing some 280,000 polygons.

Ideally, I think that it would be preferable to save the vector layer as a rotated, interpolated raster as it is essentially raster data that they have stored as a vector for convenience.

enter image description here

enter image description here

  • You appear to be using the term "vector field" in a novel, idiosyncratic way. This makes it difficult to understand what you are asking. Could you explain your meaning? I would guess you have some kind of tessellation of a region by parallelograms, and possibly there are values defined on these parallelograms--making it a form of raster data--but those are purely impressions from the graphics. – whuber Apr 11 '17 at 22:52
  • Apologies for the confusion @whuber: I am quite new to GIS and still learning the appropriate terminology. Within the vector shape file (Vessel_Density_Grid_2014.shp) there are approx 280,000 features (each is a polygon geometry). You are correct, these are parallelograms that tessellate and each cover a different area. There are 18 fields contained in the file, by "vector field" I am referring to the data values corresponding to one of these fields, the "Annual_Ave" field. – EngStan Apr 12 '17 at 14:07
  • It would be good to know whether these are exactly parallelograms or only approximately so. If they are exact, you could convert your data to a raster format, one polygon per cell, and thereby have access to (fast) raster-based solutions. BTW, here is the standard meaning of "vector field": google.com/search?q="vector+field". – whuber Apr 12 '17 at 14:09
  • Thanks for the quick reply. I did try converting the data to a raster, using a cell size of 0.004 which lined up well with the individual polygons. Do you know how I can interpolate against the resulting raster? The issue I had is that I want to do it in Python, but the geotransform corresponding to the raster data does not include the rotation of the polygons. This meant that I could not get the actual pixel coordinates – EngStan Apr 12 '17 at 14:15
  • I would do this using a local coordinate system in which each integral coordinate (i,j) designates a specific parallelogram. This amounts to an affine transformation of the original coordinates that maps parallelograms to squares. By applying that transformation (which is easily computed) to the interpolation points, you would reduce the problem to the ordinary task of interpolating values on a square grid onto arbitrary points. This can be done in myriad ways; how you do it in Python will depend on what code you have available, what the data values mean, and your accuracy needs. – whuber Apr 12 '17 at 21:53

Your Answer

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

Browse other questions tagged or ask your own question.