xyz2grd $(gmtinfo -I- *.xyz) ...<\/code><\/strong><\/p>\nupdate: xyz2grd supports GeoTIFF via its GDAL driver (now?)! For example -Gfile.tif=gd:GTiff<\/code>. See eg http:\/\/gmt.soest.hawaii.edu\/doc\/5.2.1\/grdconvert.html<\/a>. So the way below is overly complicated. I do not know how to apply the advanced settings of GDAL though like compression and prediction etc.<\/strong><\/p>\n![\"zensus2011_berlin\"](\"http:\/\/hannes.enjoys.it\/blog\/wp-content\/uploads\/zensus2011_berlin-250x194.png\")
<\/a>![\"zensus2011_irgendwo\"](\"http:\/\/hannes.enjoys.it\/blog\/wp-content\/uploads\/zensus2011_irgendwo-250x194.png\")
<\/a>![\"zensus2011_bremen\"](\"http:\/\/hannes.enjoys.it\/blog\/wp-content\/uploads\/zensus2011_bremen-250x194.png\")
<\/a><\/p>\nThe Statistische \u00c4mter des Bundes und der L\u00e4nder offer a 100 meter grid of Germany’s population density<\/a>: csv_Bevoelkerung_100m_Gitter.zip<\/a> (110MB). Datensatzbeschreibung_Bevoelkerung_100m_Gitter.xlsx<\/a> provides additional information.<\/p>\nLet’s turn that dataset into a GeoTIFF so we can use it in our GIS. We will use free and open-source tools from GMT<\/a> and GDAL<\/a>. GDAL loves to interpolate values but our data is discrete\/regular. We do not want any kind of interpolation. So xyz2grd<\/a> from GMT is the best choice for turning the xyz data into a “continuous” GIS format (tell me if not). <\/p>\nGetting to know the data<\/h2>\n
Inside the zip is a 1.3GB file Zensus_Bevoelkerung_100m-Gitter.csv with about 36 million lines.<\/p>\n
Gitter_ID_100m;x_mp_100m;y_mp_100m;Einwohner
\n100mN26840E43341;4334150;2684050;-1
\n100mN26840E43342;4334250;2684050;-1
\n…
\n100mN27407E44044;4404450;2740750;3
\n100mN27407E44045;4404550;2740750;31
\n100mN27407E44046;4404650;2740750;13
\n100mN27407E44047;4404750;2740750;14
\n100mN27407E44048;4404850;2740750;10<\/p><\/blockquote>\n
Datensatzbeschreibung_Bevoelkerung_100m_Gitter.xlsx says the coordinates are in ETRS89-LAEA Europe – EPSG:3035<\/em>.<\/p>\nFirst we need to find out the geographic extends of the data, you could use your favourite cli tools for that, I wrote a quick .vrt file and used ogrinfo on that:<\/p>\n
$ cat Zensus_Bevoelkerung_100m-Gitter.csv.vrt<\/code><\/p>\n\n\r\n<OGRVRTDataSource>\r\n <OGRVRTLayer name="Zensus_Bevoelkerung_100m-Gitter">\r\n <LayerSRS>EPSG:3035<\/LayerSRS>\r\n <SrcDataSource>Zensus_Bevoelkerung_100m-Gitter.csv<\/SrcDataSource>\r\n <GeometryType>wkbPoint<\/GeometryType>\r\n <GeometryField encoding="PointFromColumns" x="x_mp_100m" y="y_mp_100m" \/>\r\n <\/OGRVRTLayer>\r\n<\/OGRVRTDataSource>\r\n<\/pre>\n<\/blockquote>\n$ ogrinfo -al Zensus_Bevoelkerung_100m-Gitter.csv.vrt<\/code><\/p>\nINFO: Open of `Zensus_Bevoelkerung_100m-Gitter.csv.vrt’ using driver `VRT’ successful.<\/p>\n
Layer name: Zensus_Bevoelkerung_100m-Gitter
\nGeometry: Point
\nFeature Count: 35785840
\nExtent: (4031350.000000, 2684050.000000) – (4672550.000000, 3551450.000000)
\n…<\/p><\/blockquote>\n
Creating a netcdf\/grd file from the xyz data with xyz2grd<\/h2>\n
xyz2grd wants xyz, nothing else. The Gitter_ID_100m column is redundant in any case, you can calculate it yourself from the x and y fields if needed. So first let’s convert it to a “x y z” format with awk. The separator is ;<\/p>\n
awk 'FS=\";\" {print $2\" \"$3\" \"$4}' Zensus_Bevoelkerung_100m-Gitter.csv > Zensus_Bevoelkerung_100m-Gitter.xyz<\/code><\/p>\nNow we can write our xyz2grd commandline.<\/p>\n
We have our extends:
\n-R4031350\/4672550\/2684050\/3551450 <\/code><\/p>\nWe know the spacing is 100 units (meters):
\n-I100<\/code><\/p>\nThere is one header line:
\n-h1<\/code><\/p>\nAnd of course we know that we want a “classic” netcdf4 chunk size, whatever that means. We knew that right away, not after googling helplessly for an hour and eventually finding the hint on some mailing list<\/a>. Not knowing this might have lead to QGIS only seeing NaN values for z, R’s ncdf\/raster saying “Error in substr(w, 1, 3) : invalid multibyte string at ‘<89>HDF” and GDAL “0ERROR 1: nBlockYSize = 130, only 1 supported when reading bottom-up dataset”.
\n--IO_NC4_CHUNK_SIZE=c<\/code><\/p>\nThe resulting commandline:
\nxyz2grd -Vl -R4031350\/4672550\/2684050\/3551450 -I100 -h1 --IO_NC4_CHUNK_SIZE=c -GZensus_Bevoelkerung_100m-Gitter.cdf Zensus_Bevoelkerung_100m-Gitter.xyz<\/code><\/p>\nYou can inspect the file with grdinfo and gdalinfo now if you want.<\/p>\n
Creating a GeoTIFF from the netcdf\/grd file<\/h2>\n
Let’s turn it into a GeoTIFF with gdal_translate<\/a>. We will need a bunch of commandline parameters.<\/p>\nThe spatial reference system is EPSG:3035, so:
\n-a_srs EPSG:3035<\/code><\/p>\nValues of -1 mean “no data”:
\n-a_nodata -1<\/code><\/p>\nTIFF is uncompressed by default, we want good lossless compression:
\n-co COMPRESS=DEFLATE<\/code><\/p>\nThe resulting commandline:
\ngdal_translate -co COMPRESS=DEFLATE -a_srs EPSG:3035 -a_nodata -1 Zensus_Bevoelkerung_100m-Gitter.cdf Zensus_Bevoelkerung_100m-Gitter.tif<\/code><\/p>\nThe resulting file is about 8 Megabytes and should work in any reasonable GIS. Have fun!<\/p>\n
Downloads<\/h2>\n
http:\/\/hannes.enjoys.it\/opendata\/Zensus_Bevoelkerung_100m-Gitter.tif<\/a>
\nTODO: What license is this now?<\/p>\n![\"zensus2011_hamburg\"](\"http:\/\/hannes.enjoys.it\/blog\/wp-content\/uploads\/zensus2011_hamburg.png\")
\n![\"zensus2011_altersheim\"](\"http:\/\/hannes.enjoys.it\/blog\/wp-content\/uploads\/zensus2011_altersheim-250x188.jpg\")
<\/a>![\"zensus2011_schreberg\u00e4rten\"](\"http:\/\/hannes.enjoys.it\/blog\/wp-content\/uploads\/zensus2011_schreberg\u00e4rten-250x188.jpg\")
<\/a><\/p>\n","protected":false},"excerpt":{"rendered":"tl;dr: GMT is documented for people who use it since the 80s. update 2: You can use gmtinfo to calculate the extents, using -I- will make it output the -R parameter! xyz2grd $(gmtinfo -I- *.xyz) … update: xyz2grd supports GeoTIFF via its GDAL driver (now?)! For example -Gfile.tif=gd:GTiff. See eg http:\/\/gmt.soest.hawaii.edu\/doc\/5.2.1\/grdconvert.html. So the way below […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[9,22,29,10],"tags":[],"class_list":["post-448","post","type-post","status-publish","format-standard","hentry","category-cartography","category-commandline","category-gis","category-open-data"],"_links":{"self":[{"href":"https:\/\/hannes.enjoys.it\/blog\/wp-json\/wp\/v2\/posts\/448","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/hannes.enjoys.it\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/hannes.enjoys.it\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/hannes.enjoys.it\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/hannes.enjoys.it\/blog\/wp-json\/wp\/v2\/comments?post=448"}],"version-history":[{"count":30,"href":"https:\/\/hannes.enjoys.it\/blog\/wp-json\/wp\/v2\/posts\/448\/revisions"}],"predecessor-version":[{"id":942,"href":"https:\/\/hannes.enjoys.it\/blog\/wp-json\/wp\/v2\/posts\/448\/revisions\/942"}],"wp:attachment":[{"href":"https:\/\/hannes.enjoys.it\/blog\/wp-json\/wp\/v2\/media?parent=448"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/hannes.enjoys.it\/blog\/wp-json\/wp\/v2\/categories?post=448"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/hannes.enjoys.it\/blog\/wp-json\/wp\/v2\/tags?post=448"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}