Convert raster to CSV with lat, lon, and value columns

Question

I'm looking for a process to convert ASCII gridded data (in this case 60min/1 degree gridded ASCII population data (GPW) from SEDAC: https://sedac.ciesin.columbia.edu/data/set/gpw-v4-population-count-adjusted-to-2015-unwpp-country-totals-rev11/data-download) to a CSV with three columns: lon, lat, and value using open source tools.

This is the data header:

ncols 360
nrows 180
xllcorner -180
yllcorner -90
cellsize 1.0000000000001
NODATA_value -9999

so the data is structured as space-separated values, with each value occupying a position within 360 columns and 180 rows that represents its coordinates via that position:

-9999 -9999 -9999 -9999 -9999 -9999 0.000936705 0.002529013 0.001377391 0.001723122 0.0004472999 ... 

I only want to include positive values in the lon/lat CSV.

I've been using GDAL for other tasks and have been looking in the documentation, but I'm not seeing this type of data manipulation in its scope.

Answer 1

You can do this with rioxarray:

Note: I used data for the column name as values is a pandas specific attribute.

import rioxarray
import pandas
rds = rioxarray.open_rasterio(
 "gpw_v4_population_count_adjusted_to_2015_unwpp_country_totals_rev11_2020_1_deg.tif",
)
rds = rds.squeeze().drop("spatial_ref").drop("band")
rds.name = "data"
df = rds.to_dataframe().reset_index()
df[df.data>=0.0].to_csv("out.csv", index=False)

• interwebjill: This solution backtracked on my process a bit by using the GeoTiff downloaded from SEDAC instead of the ASCII file. It's a great solution...much better than the brute force python method I was writing for the ASCII file. I wonder if I should change the title of the question to make it more general so others can find it.
• snowman2: Glad to hear that it worked - it should work the same with an ASCII file.
• snowman2: Title could be: Convert raster to CSV with lat, lon, and value columns
• interwebjill: Will rioxarray translate between raster and geojson point features as well? I didn't see that in the documentation. On a lark I tried df[df.data>0.0].to_geojson("out.geojson", index=False) (didn't work).
• snowman2: This may help: gis.stackexchange.com/questions/220997/…

Answer 2

Even easier using GDAL from the command line, nothing to install!

1. Create a text file with your longitudes and latitudes, call it xy.txt:

-7, 41
-0.5, 51.4

2. Create the header of your shiny CSV file

echo Longitude,Latitude,Pop > final.csv

3. Then pipe this through to gdallocationinfo:

cat xy.txt | tr -d , | gdallocationinfo -wgs84 -valonly \
 gpw_v4_population_count_rev11_2010_2pt5_min.tif > z.txt

3. Paste xy.txt and z.txt together and append to output final.csv

paste -d, xy.txt z.txt >> final.csv

You can then use ogr2ogr to convert the CSV to GeoJSON (see here)

ogr2ogr -f GeoJSON out.geojson final.csv \
 -oo X_POSSIBLE_NAMES=Longitude \
 -oo Y_POSSIBLE_NAMES=Latitude \
 -oo KEEP_GEOM_COLUMNS=NO

From ogrinfo -al out.geojson I get:

INFO: Open of `out.geojson'
 using driver `GeoJSON' successful.
Layer name: final
Geometry: Point
Feature Count: 2
Extent: (-7.000000, 41.000000) - (-0.500000, 51.400000)
Layer SRS WKT:
GEOGCRS["WGS 84",
 DATUM["World Geodetic System 1984",
 ELLIPSOID["WGS 84",6378137,298.257223563,
 LENGTHUNIT["metre",1]]],
 PRIMEM["Greenwich",0,
 ANGLEUNIT["degree",0.0174532925199433]],
 CS[ellipsoidal,2],
 AXIS["geodetic latitude (Lat)",north,
 ORDER[1],
 ANGLEUNIT["degree",0.0174532925199433]],
 AXIS["geodetic longitude (Lon)",east,
 ORDER[2],
 ANGLEUNIT["degree",0.0174532925199433]],
 ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Pop: String (0.0)
OGRFeature(final):0
 Pop (String) = 4.11418008804321
 POINT (-7 41)
OGRFeature(final):1
 Pop (String) = 12097.869140625
 POINT (-0.5 51.4)