[GIS] Georeferencing uneven/irregularly gridded rasters using GDAL

coordinate systemgdalgeoreferencingraster

I'm trying to take a set of NetCDFs and turn them into georeferenced datasets using GDAL.

However, the only Geospatial data is a set of 2-D lat/lon arrays. The projection type is not lat/lon. THe arrays give the lat/lon of every individual cell.

I can't use a geotransform to reference the data because of the irregular gridding.

How can I create a dataset (such as a GeoTiff) that properly references the data?

Ideally I would then warp to some even grid.

I tried using an even spread of Ground Control Points, but the resulting GeoTiff did not display correctly… Despite giving GCPs up to 90 deg lat, ArcMap claims the GTiff has an extent ending at 70 deg latitude…

I was able to solve the issue using gdalwarp -geoloc as suggested. I created 3 VRTs (1 with cell values, 1 with latitude, 1 with longitude), and wrote the lat/lon vrts in a geoarrays in the data vrt metadata. Then using gdalwarp -geoloc worked like a charm.

The issue with extent was unrelated and due simply to the large difference between the original coordinate system and lat/lon. It was solved by explicitely stating the extent instead of having GDAL try to guess what it should be. (i.e. I added -te -180 0 180 90 to the gdalwarp line)

Best Answer

gdalwarp -geoloc allows you to use the complete 2d-array of latlon as georeference. With that, you can use any target CRS to reproject your data to a commonly used projection.

See my answer to this question for an example:

How to match a raster NetCDF data with a vector layer in QGIS?

Related Solutions

[GIS] NumPy to GeoTiff for use with Gdal

Can be done fairly easy. First, create the dataset:

drv = gdal.GetDriverByName("GTiff")
ds = drv.Create("name.tif", width, height, 6, gdal.GDT_Float32)

You then need to set the geotransform using ds.SetGeoTransform where the argument is a six element tuple: (upper_left_x, x_resolution, x_skew, upper_left_y, y_skew, y_resolution)

Set the projection using ds.SetProjection, argument is the way string, then Write the array for the first band:

ds.GetRasterBand(1).WriteArray(arr)

That should be enough to get you going. If the purpose is just to reproject the data, then you could look at using an in memory dataset

QGIS – Saving Mapzen Global Terrain with Extent Layer

No. Mapzen Global Terrain tiles will be designed for use in Web Mercator projection, which only goes to ±85.05.. degrees of latitude. From Wikipedia:

Because the Mercator projects the poles at infinity, a map using the Web Mercator projection cannot show the poles. Services such as Google Maps cut off coverage at 85.051129° north and south. This is not a limitation for street maps, which is the primary purpose for such services. The value 85.051129° is the latitude at which the full projected map becomes a square...

So it has one useful property: all tiles are square, including the one at zoom 0. But, the point at ±90° would be stretched to infinity in order to be rendered as a line along the top or bottom of the map. So the projection just cuts off before that point. That is why Antarctica looks so large, when it's really about the same size as Australia.

If you need estimates of elevation at the poles, you will need to look elsewhere. However if you only ever intend to do stuff on webmaps projected in Web Mercator, there's no point obtaining the remaining data since you won't be able to see it.

Related Question