0

I want to use an old sofware to analyse my data with a new angle. This software is not compatible with the .tif format so I need to transform all my data into a .csv.

The expected .csv structure is the following:

ADM0_CODE pixel_area lat lon layer1_name layer2_name ... layerN_name

I've prepared the data so that all this information (including the 4 first specific columns) are on a list of .tif files with exact same number of pixels, crs etc.

I've tried many things involving for loops that are never going to finish and that are eating up all my computer memory.

Is there a canonic way to do it with gdal and/or rasterio ?

The maps are covering half he world with 1km grid cells so each individual file is about 1.9 Gb and there are 29.

5
  • I can help you to extrac values from a raster and the area, Is it ok for you?
    – Helios
    Mar 1, 2022 at 21:09
  • Do you want to calculate the total area?
    – Helios
    Mar 1, 2022 at 21:13
  • I already have computed the pixel_area of every pixel and stored it in a raster Mar 1, 2022 at 21:14
  • I can show you a simple example about how to extrat the raster values with geopandas, Do you agree?
    – Helios
    Mar 1, 2022 at 21:30
  • I mean, yes of course, if you have any solution you can post an answer Mar 1, 2022 at 21:39

2 Answers 2

1

Here's a rasterio method that reads your csv and uses the dataset.sample method for each of the input rasters.

Note sample accepts an iterable of coordinates and returns a generator (hence the next(sample([[lat, lon]])) in the code below).

This is pretty slow as it's reading the coordinates and values one by one. It may be quicker to just use gdal_translate to convert to x,y,z format then paste the z values to your input csv. But I've done it this way rather than reading all the coordinates into memory at once because of the size.

import csv
from pathlib import Path

import rasterio

dir = Path('path/to/raster_dir')

incsvfile = 'path/to/in.csv'
outcsvfile = 'path/to/out.csv'
rasters = ['test1.tif', 'test2.tif', 'test3.tif']

with open(incsvfile) as incsvfile, open(outcsvfile, 'w', newline='') as outcsvfile:
    csvreader = csv.reader(incsvfile)
    csvwriter = csv.writer(outcsvfile)

    fieldnames = next(csvreader)
    csvwriter.writerow(fieldnames)

    datasets = [rasterio.open(dir/raster) for raster in rasters]

    for adm0_code, pixel_area, lat, lon in csvreader:
        values = []
        for src in datasets:
            values.append(next(src.sample([[float(lon), float(lat)]]))[0])

        csvwriter.writerow([adm0_code, pixel_area, lat, lon] + values)

And here's a way using the stack_vrts function I showed you previously, it's slightly faster.

import xml.etree.ElementTree as ET

import rasterio as rio
from rasterio.shutil import copy as riocopy
from rasterio.io import MemoryFile


def stack_vrts(srcs, band=1):
    vrt_bands = []
    for srcnum, src in enumerate(srcs, start=1):
        with rio.open(src) as ras, MemoryFile() as mem:
            riocopy(ras, mem.name, driver='VRT')
            vrt_xml = mem.read().decode('utf-8')
            vrt_dataset = ET.fromstring(vrt_xml)
            for bandnum, vrt_band in enumerate(vrt_dataset.iter('VRTRasterBand'), start=1):
                if bandnum == band:
                    vrt_band.set('band', str(srcnum))
                    vrt_bands.append(vrt_band)
                    vrt_dataset.remove(vrt_band)
    for vrt_band in vrt_bands:
        vrt_dataset.append(vrt_band)

    return ET.tostring(vrt_dataset).decode('UTF-8')


dir = Path('path/to/raster_dir')

incsvfile = 'path/to/in.csv'
outcsvfile = 'path/to/out.csv'
rasters = ['test1.tif', 'test2.tif', 'test3.tif']

with open(incsvfile) as incsvfile, open(outcsvfile, 'w', newline='') as outcsvfile:
    csvreader = csv.reader(incsvfile)
    csvwriter = csv.writer(outcsvfile)

    fieldnames = next(csvreader)
    csvwriter.writerow(fieldnames)

    with rio.open(stack_vrts([dir/raster for raster in rasters])) as src:

        for adm0_code, pixel_area, lat, lon in csvreader:
            values = next(src.sample([[float(lon), float(lat)]]))
            csvwriter.writerow([adm0_code, pixel_area, lat, lon] + list(values))
0
# Import modules

import rasterio 
from rasterio import plot
import geopandas as gpd

# Read data (raster and points)

raster = rasterio.open(inputpath_raster)
points = gpd.read_file(inputpath_shapefile_or_geojson_points)

# Plot raster and points

rasterio.plot.show(raster)
print(points)

# Get dataframe with points and coords

coord_list = [(x,y) for x,y in zip(points['geometry'].x , points['geometry'].y)]

points['value'] = [x for x in raster.sample(coord_list)]
points['x'] = points['geometry'].x
points['y'] = points['geometry'].y
points["value"] = points["value"].astype("int64") # or float64
points

# Save csv

points.to_csv(outputpath_csv)
1
  • thanks for the solution. It could work if I had a list of point but here I want instead to download each pixel of the initial raster. There size and the amont of raster to work on mak me look for something else than a for loop solution. I've edited to question with my current code to reflect this idea. Mar 6, 2022 at 11:31

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.