From b2e5d4a9b7adb698cbde2ab36d283098b2038d2a Mon Sep 17 00:00:00 2001 From: Pritam Date: Fri, 10 Jul 2020 11:07:26 +0530 Subject: [PATCH 1/7] Rasterio style Tutorial --- docs/src/quickstart.md | 132 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 132 insertions(+) create mode 100644 docs/src/quickstart.md diff --git a/docs/src/quickstart.md b/docs/src/quickstart.md new file mode 100644 index 00000000..097c10ad --- /dev/null +++ b/docs/src/quickstart.md @@ -0,0 +1,132 @@ +# ArchGDAL quickstart +>This tutorial is highly inspired by the [Python Quickstart](https://rasterio.readthedocs.io/en/latest/quickstart.html), and aims to provide a similar introductory guide to the ArchGDAL package. + +> The output of each code block is represented as a series of comments after each code block. For example, +>```Julia +> println("Hello World") +># Hello World +>``` + +This guide uses a GeoTIFF image of the band 4 of a Landsat 8 image, covering a part of the United State's Colorado Plateau. This image was acquired on July 12, 2016, and can be download at [this link](https://landsatonaws.com/L8/037/034/LC08_L1TP_037034_20160712_20170221_01_T1). +![example.tif](https://user-images.githubusercontent.com/4471859/87169013-a32ee600-c2cf-11ea-9d09-e82446812282.png) + +## Opening the dataset +Assuming that the ArchGDAL package is already installed, first import the package, and assign an alias for the package. +```Julia +using ArchGDAL; const ag = ArchGDAL +``` +Assuming that the dataset exists in the following path: `./example.tiff`, the dataset can be opened in reading mode as follows - +```Julia +dataset = ag.read("example.tiff") +# GDAL Dataset (Driver: GTiff/GeoTIFF) +# File(s): +# example.tiff + +# Dataset (width x height): 7731 x 7861 (pixels) +# Number of raster bands: 1 +# [GA_ReadOnly] Band 1 (Gray): 7731 x 7861 (UInt16) +``` +A dataset can be opened using ArchGDAL's `read(...)` function. This function takes in the path of a dataset as a string, and returns a GDAL Dataset object in read-only mode. + +ArchGDAL can read a variety of file types using appropriate drivers. The driver used for opening the dataset can be queried as such. +```Julia +ag.getdriver(dataset) +# Driver: GTiff/GeoTIFF +``` + +## Dataset Attributes +A raster dataset can have multiple bands of information. To get the number of bands present in a dataset - +```Julia +ag.nraster(dataset) +# 1 +``` +In our case, the raster dataset only has a single band. + +Other properties of a raster dataset, such as the width and height can be accessed using the following functions +```Julia +ag.width(dataset) +# 7731 +ag.height(dataset) +# 7861 +``` + +## Dataset Georeferencing +All raster datasets contain embedded georeferencing information, using which the raster image can be located at a geographic location. The georeferencing attributes are stored as the dataset's Geospatial Transform. +```Julia +gt = ag.getgeotransform(dataset) +# 6-element Array{Float64,1}: +# 358485.0 +# 30.0 +# 0.0 +# 4.265115e6 +# 0.0 +# -30.0 +``` +The x-y coordinates of the top-left corner of the raster dataset are denoted by `gt[1]` and `gt[4]` values. In this case, the coordinates . The cell size is represented by `gt[2]` and `gt[6]` values, corresponding to the cell size in x and y directions respectively, and the `gt[3]` and `gt[5]` values define the rotation. In our case, the top left pixel is at an offset of 358485.0m and 4.265115e6m in x and y directions respectively, from the origin. + +The Origin of the dataset is defined using what is called, a Coordinate Reference System (CRS in short). +```Julia +p = ag.getproj(dataset) +# PROJCS["WGS 84 / UTM zone 12N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32612"]] +``` +The `getproj(...)` function returns the projection of a dataset as a string representing the CRS in the OpenGIS WKT format. To convert it into other CRS representations, such as PROJ.4, `toPROJ4(...)` can be used. However, first the string representation of the CRS has to be converted into an `ArchGDAL.ISpatialRef` type using the `importWKT(...)` function. +```Julia +ag.toPROJ4(ag.importWKT(p)) +# "+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs" +``` + +## Reading Raster data +In the previous steps, we read the raster file as a dataset. To get the actual raster data, stored as bands, we can again use the `read(...)` function. +```Julia +ag.read(dataset, 1) +# 7731×7861 Array{UInt16,2}: +# 0x0000 0x0000 0x0000 0x0000 0x0000 … 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 … 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 … 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# ⋮ ⋱ ⋮ +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 … 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 … 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 +# 0x0000 0x0000 0x0000 0x0000 0x0000 … 0x0000 0x0000 0x0000 0x0000 +``` +The function takes in a dataset as the first argument, and the index of the band that is to be read, as the second argument. To read the whole dataset, `read(dataset)` can be used, in which case, a single multidimensional array will be returned, containing all the bands. The indices will correspond to `(cols, rows, bands)`. + +To get the individual band +```Julia +band = ag.getband(dataset, 1) +# [GA_ReadOnly] Band 1 (Gray): 7731 x 7861 (UInt16) +# blocksize: 512×512, nodata: -1.0e10, units: 1.0px + 0.0 +# overviews: +``` +Band specific attributes can be easily accessed through a variety of functions +```Julia +ag.width(band) +# 7731 +ag.height(band) +# 7861 +ag.getnodatavalue(band) +# -1.0e10 +ag.indexof(band) +# 1 +``` +The no-data value of the band can be obtained using the `getnodatavalue(...)` function, which in our case is -1.0e10. + +## Creating Data From 65e5f356bd6dabe36900b324ee8e0150f130faad Mon Sep 17 00:00:00 2001 From: Pritam Date: Thu, 16 Jul 2020 13:17:47 +0530 Subject: [PATCH 2/7] Writing raster dataset (incomplete) --- docs/src/quickstart.md | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/docs/src/quickstart.md b/docs/src/quickstart.md index 097c10ad..c03bf88c 100644 --- a/docs/src/quickstart.md +++ b/docs/src/quickstart.md @@ -129,4 +129,17 @@ ag.indexof(band) ``` The no-data value of the band can be obtained using the `getnodatavalue(...)` function, which in our case is -1.0e10. -## Creating Data +## Writing into Dataset +Creating dummy data to be written as raster dataset. +```Julia +using Plots +x = range(-4.0, 4.0, length=240) +y = range(-3.0, 3.0, length=180) +X = repeat(collect(x)', size(y)[1]) +Y = repeat(collect(y), 1, size(x)[1]) +Z1 = exp.(-2*log(2) .*((X.-0.5).^2 + (Y.-0.5).^2)./(1^2)) +Z2 = exp.(-3*log(2) .*((X.+0.5).^2 + (Y.+0.5).^2)./(2.5^2)) +Z = 10.0 .* (Z2 .- Z1) +``` +The fictional field for this example consists of the difference of two Gaussian distributions and is represented by the array `Z`. Its contours are shown below. +![creating-dataset](https://user-images.githubusercontent.com/7526346/87633084-5bd5a900-c758-11ea-8fd3-548d039f1a43.png) From 71992901fb59a1b33a582a33cb924c768abffee5 Mon Sep 17 00:00:00 2001 From: Pritam Date: Thu, 16 Jul 2020 13:20:52 +0530 Subject: [PATCH 3/7] refactor --- docs/src/quickstart.md | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/docs/src/quickstart.md b/docs/src/quickstart.md index c03bf88c..4456527b 100644 --- a/docs/src/quickstart.md +++ b/docs/src/quickstart.md @@ -13,11 +13,11 @@ This guide uses a GeoTIFF image of the band 4 of a Landsat 8 image, covering a p ## Opening the dataset Assuming that the ArchGDAL package is already installed, first import the package, and assign an alias for the package. ```Julia -using ArchGDAL; const ag = ArchGDAL +using ArchGDAL; const AG = ArchGDAL ``` Assuming that the dataset exists in the following path: `./example.tiff`, the dataset can be opened in reading mode as follows - ```Julia -dataset = ag.read("example.tiff") +dataset = AG.read("example.tiff") # GDAL Dataset (Driver: GTiff/GeoTIFF) # File(s): # example.tiff @@ -30,30 +30,30 @@ A dataset can be opened using ArchGDAL's `read(...)` function. This function tak ArchGDAL can read a variety of file types using appropriate drivers. The driver used for opening the dataset can be queried as such. ```Julia -ag.getdriver(dataset) +AG.getdriver(dataset) # Driver: GTiff/GeoTIFF ``` ## Dataset Attributes A raster dataset can have multiple bands of information. To get the number of bands present in a dataset - ```Julia -ag.nraster(dataset) +AG.nraster(dataset) # 1 ``` In our case, the raster dataset only has a single band. Other properties of a raster dataset, such as the width and height can be accessed using the following functions ```Julia -ag.width(dataset) +AG.width(dataset) # 7731 -ag.height(dataset) +AG.height(dataset) # 7861 ``` ## Dataset Georeferencing All raster datasets contain embedded georeferencing information, using which the raster image can be located at a geographic location. The georeferencing attributes are stored as the dataset's Geospatial Transform. ```Julia -gt = ag.getgeotransform(dataset) +gt = AG.getgeotransform(dataset) # 6-element Array{Float64,1}: # 358485.0 # 30.0 @@ -66,19 +66,19 @@ The x-y coordinates of the top-left corner of the raster dataset are denoted by The Origin of the dataset is defined using what is called, a Coordinate Reference System (CRS in short). ```Julia -p = ag.getproj(dataset) +p = AG.getproj(dataset) # PROJCS["WGS 84 / UTM zone 12N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32612"]] ``` The `getproj(...)` function returns the projection of a dataset as a string representing the CRS in the OpenGIS WKT format. To convert it into other CRS representations, such as PROJ.4, `toPROJ4(...)` can be used. However, first the string representation of the CRS has to be converted into an `ArchGDAL.ISpatialRef` type using the `importWKT(...)` function. ```Julia -ag.toPROJ4(ag.importWKT(p)) +AG.toPROJ4(AG.importWKT(p)) # "+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs" ``` ## Reading Raster data In the previous steps, we read the raster file as a dataset. To get the actual raster data, stored as bands, we can again use the `read(...)` function. ```Julia -ag.read(dataset, 1) +AG.read(dataset, 1) # 7731×7861 Array{UInt16,2}: # 0x0000 0x0000 0x0000 0x0000 0x0000 … 0x0000 0x0000 0x0000 0x0000 # 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 @@ -111,20 +111,20 @@ The function takes in a dataset as the first argument, and the index of the band To get the individual band ```Julia -band = ag.getband(dataset, 1) +band = AG.getband(dataset, 1) # [GA_ReadOnly] Band 1 (Gray): 7731 x 7861 (UInt16) # blocksize: 512×512, nodata: -1.0e10, units: 1.0px + 0.0 # overviews: ``` Band specific attributes can be easily accessed through a variety of functions ```Julia -ag.width(band) +AG.width(band) # 7731 -ag.height(band) +AG.height(band) # 7861 -ag.getnodatavalue(band) +AG.getnodatavalue(band) # -1.0e10 -ag.indexof(band) +AG.indexof(band) # 1 ``` The no-data value of the band can be obtained using the `getnodatavalue(...)` function, which in our case is -1.0e10. From fdf58cee7fe98b70d4f450792b87d5180de44412 Mon Sep 17 00:00:00 2001 From: Pritam Date: Thu, 16 Jul 2020 13:21:09 +0530 Subject: [PATCH 4/7] Add page to make.jl --- docs/make.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/make.jl b/docs/make.jl index 257cd981..7afd7c27 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -9,6 +9,7 @@ makedocs( strict = true, pages = [ "index.md", + "ArchGDAL Quickstart Guide" => "quickstart.md" "GDAL Datasets" => "datasets.md", "Feature Data" => "features.md", "Raster Data" => "rasters.md", From 84686701f44392e754d47685e2b35604848ee9e0 Mon Sep 17 00:00:00 2001 From: Pritam Date: Thu, 20 Aug 2020 12:50:45 +0530 Subject: [PATCH 5/7] writing dataset from scratch --- docs/src/quickstart.md | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/docs/src/quickstart.md b/docs/src/quickstart.md index 4456527b..a1bed5ba 100644 --- a/docs/src/quickstart.md +++ b/docs/src/quickstart.md @@ -143,3 +143,38 @@ Z = 10.0 .* (Z2 .- Z1) ``` The fictional field for this example consists of the difference of two Gaussian distributions and is represented by the array `Z`. Its contours are shown below. ![creating-dataset](https://user-images.githubusercontent.com/7526346/87633084-5bd5a900-c758-11ea-8fd3-548d039f1a43.png) + +Once we have the data to write into a raster dataset, before performing the write operation itself, it is a good idea to define the geotransform and the coordinate reference system of the resulting dataset. For North-up raster images, only the *pixel-size/resolution* and the coordinates of the *upper-left* pixel is required. +```Julia +# Resolution +res = (x[end] - x[1]) /240.0 + +# Upper-left pixel coordinates +ul_x = x[1] - res/2 +ul_y = y[1] - res/2 + +gt = [ul_x, res, 0.0, ul_y, 0.0, res] +``` +The coordinate reference system of the dataset needs to be a string in the WKT format. +```Julia +crs = AG.toWKT(AG.importPROJ4("+proj=latlong")) +``` + +To write this array, first a dataset has to be created. This can be done using the `AG.create` function. The first argument defines the path of the resulting dataset. The following keyword arguments are hopefully self-explanatory. + +Once inside the ```do...end``` block, the `write!` method can be used to write the array(`Z`), in the 1st band of the opened `dataset`. Next, the geotransform and CRS can be specified using the `setgeotransform!` and `setproj!` methods. +```Julia +AG.create( + "./temporary.tif", + driver = AG.getdriver("GTiff"), + width=240, + height=180, + nbands=1, + dtype=Float64 +) do dataset + AG.write!(dataset, Z, 1) + + AG.setgeotransform!(dataset, gt) + AG.setproj!(dataset, crs) +end +``` \ No newline at end of file From 8b661cb32c41ceff2bc7a0262e1a328f513afa45 Mon Sep 17 00:00:00 2001 From: Pritam Date: Fri, 21 Aug 2020 08:42:23 +0530 Subject: [PATCH 6/7] Make documentation Documenter.jl friendly --- docs/src/quickstart.md | 104 +++++++++-------------------------------- 1 file changed, 22 insertions(+), 82 deletions(-) diff --git a/docs/src/quickstart.md b/docs/src/quickstart.md index a1bed5ba..eb48fedc 100644 --- a/docs/src/quickstart.md +++ b/docs/src/quickstart.md @@ -1,138 +1,78 @@ # ArchGDAL quickstart >This tutorial is highly inspired by the [Python Quickstart](https://rasterio.readthedocs.io/en/latest/quickstart.html), and aims to provide a similar introductory guide to the ArchGDAL package. -> The output of each code block is represented as a series of comments after each code block. For example, ->```Julia -> println("Hello World") -># Hello World ->``` This guide uses a GeoTIFF image of the band 4 of a Landsat 8 image, covering a part of the United State's Colorado Plateau. This image was acquired on July 12, 2016, and can be download at [this link](https://landsatonaws.com/L8/037/034/LC08_L1TP_037034_20160712_20170221_01_T1). ![example.tif](https://user-images.githubusercontent.com/4471859/87169013-a32ee600-c2cf-11ea-9d09-e82446812282.png) ## Opening the dataset -Assuming that the ArchGDAL package is already installed, first import the package, and assign an alias for the package. -```Julia +Assuming that the ArchGDAL package is already installed, first import the package, and preferably assign an alias to the package. +```@setup quickstart using ArchGDAL; const AG = ArchGDAL ``` -Assuming that the dataset exists in the following path: `./example.tiff`, the dataset can be opened in reading mode as follows - -```Julia -dataset = AG.read("example.tiff") -# GDAL Dataset (Driver: GTiff/GeoTIFF) -# File(s): -# example.tiff - -# Dataset (width x height): 7731 x 7861 (pixels) -# Number of raster bands: 1 -# [GA_ReadOnly] Band 1 (Gray): 7731 x 7861 (UInt16) +In this example, we'll be fetching the raster dataset from where it is hosted. +```@repl quickstart +dataset = AG.readraster("/vsicurl/http://landsat-pds.s3.amazonaws.com/c1/L8/037/034/LC08_L1TP_037034_20160712_20170221_01_T1/LC08_L1TP_037034_20160712_20170221_01_T1_B4.TIF") ``` -A dataset can be opened using ArchGDAL's `read(...)` function. This function takes in the path of a dataset as a string, and returns a GDAL Dataset object in read-only mode. +A dataset can be opened using ArchGDAL's `readraster(...)` method. This method takes in the path of a dataset as a string, and returns a GDAL Dataset object in read-only mode. ArchGDAL can read a variety of file types using appropriate drivers. The driver used for opening the dataset can be queried as such. -```Julia +```@repl quickstart AG.getdriver(dataset) -# Driver: GTiff/GeoTIFF ``` ## Dataset Attributes A raster dataset can have multiple bands of information. To get the number of bands present in a dataset - -```Julia +```@repl quickstart AG.nraster(dataset) -# 1 ``` In our case, the raster dataset only has a single band. Other properties of a raster dataset, such as the width and height can be accessed using the following functions -```Julia +```@repl quickstart AG.width(dataset) -# 7731 AG.height(dataset) -# 7861 ``` ## Dataset Georeferencing All raster datasets contain embedded georeferencing information, using which the raster image can be located at a geographic location. The georeferencing attributes are stored as the dataset's Geospatial Transform. -```Julia +```@repl quickstart gt = AG.getgeotransform(dataset) -# 6-element Array{Float64,1}: -# 358485.0 -# 30.0 -# 0.0 -# 4.265115e6 -# 0.0 -# -30.0 ``` The x-y coordinates of the top-left corner of the raster dataset are denoted by `gt[1]` and `gt[4]` values. In this case, the coordinates . The cell size is represented by `gt[2]` and `gt[6]` values, corresponding to the cell size in x and y directions respectively, and the `gt[3]` and `gt[5]` values define the rotation. In our case, the top left pixel is at an offset of 358485.0m and 4.265115e6m in x and y directions respectively, from the origin. The Origin of the dataset is defined using what is called, a Coordinate Reference System (CRS in short). -```Julia +```@repl quickstart p = AG.getproj(dataset) -# PROJCS["WGS 84 / UTM zone 12N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32612"]] ``` The `getproj(...)` function returns the projection of a dataset as a string representing the CRS in the OpenGIS WKT format. To convert it into other CRS representations, such as PROJ.4, `toPROJ4(...)` can be used. However, first the string representation of the CRS has to be converted into an `ArchGDAL.ISpatialRef` type using the `importWKT(...)` function. -```Julia +```@repl quickstart AG.toPROJ4(AG.importWKT(p)) -# "+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs" ``` ## Reading Raster data -In the previous steps, we read the raster file as a dataset. To get the actual raster data, stored as bands, we can again use the `read(...)` function. -```Julia -AG.read(dataset, 1) -# 7731×7861 Array{UInt16,2}: -# 0x0000 0x0000 0x0000 0x0000 0x0000 … 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 … 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 … 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# ⋮ ⋱ ⋮ -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 … 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 … 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 0x0000 -# 0x0000 0x0000 0x0000 0x0000 0x0000 … 0x0000 0x0000 0x0000 0x0000 +In the previous steps, we read the raster file as a dataset. To obtain the actual raster data, we can slice the array accordingly. +```@repl quickstart +dataset[:, :, 1] ``` -The function takes in a dataset as the first argument, and the index of the band that is to be read, as the second argument. To read the whole dataset, `read(dataset)` can be used, in which case, a single multidimensional array will be returned, containing all the bands. The indices will correspond to `(cols, rows, bands)`. +Since the dimensions of the dataset are `[cols, rows, bands]`, respectively, using `[:, :, 1]` returns all the columns and rows for the first band. Accordingly other complex slice operations can also be performed. To get the individual band -```Julia +```@repl quickstart band = AG.getband(dataset, 1) -# [GA_ReadOnly] Band 1 (Gray): 7731 x 7861 (UInt16) -# blocksize: 512×512, nodata: -1.0e10, units: 1.0px + 0.0 -# overviews: ``` Band specific attributes can be easily accessed through a variety of functions -```Julia +```@repl quickstart AG.width(band) -# 7731 AG.height(band) -# 7861 AG.getnodatavalue(band) -# -1.0e10 AG.indexof(band) -# 1 ``` The no-data value of the band can be obtained using the `getnodatavalue(...)` function, which in our case is -1.0e10. ## Writing into Dataset Creating dummy data to be written as raster dataset. -```Julia -using Plots +```@example quickstart x = range(-4.0, 4.0, length=240) y = range(-3.0, 3.0, length=180) X = repeat(collect(x)', size(y)[1]) @@ -145,7 +85,7 @@ The fictional field for this example consists of the difference of two Gaussian ![creating-dataset](https://user-images.githubusercontent.com/7526346/87633084-5bd5a900-c758-11ea-8fd3-548d039f1a43.png) Once we have the data to write into a raster dataset, before performing the write operation itself, it is a good idea to define the geotransform and the coordinate reference system of the resulting dataset. For North-up raster images, only the *pixel-size/resolution* and the coordinates of the *upper-left* pixel is required. -```Julia +```@example quickstart # Resolution res = (x[end] - x[1]) /240.0 @@ -156,14 +96,14 @@ ul_y = y[1] - res/2 gt = [ul_x, res, 0.0, ul_y, 0.0, res] ``` The coordinate reference system of the dataset needs to be a string in the WKT format. -```Julia +```@example quickstart crs = AG.toWKT(AG.importPROJ4("+proj=latlong")) ``` To write this array, first a dataset has to be created. This can be done using the `AG.create` function. The first argument defines the path of the resulting dataset. The following keyword arguments are hopefully self-explanatory. Once inside the ```do...end``` block, the `write!` method can be used to write the array(`Z`), in the 1st band of the opened `dataset`. Next, the geotransform and CRS can be specified using the `setgeotransform!` and `setproj!` methods. -```Julia +```@example quickstart AG.create( "./temporary.tif", driver = AG.getdriver("GTiff"), From 94fd2abf4ec94fc8d01d3d30d02dcfe0f2b8c400 Mon Sep 17 00:00:00 2001 From: Pritam Date: Fri, 21 Aug 2020 08:42:45 +0530 Subject: [PATCH 7/7] Fix typo --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index 7afd7c27..15bf2327 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -9,7 +9,7 @@ makedocs( strict = true, pages = [ "index.md", - "ArchGDAL Quickstart Guide" => "quickstart.md" + "ArchGDAL Quickstart Guide" => "quickstart.md", "GDAL Datasets" => "datasets.md", "Feature Data" => "features.md", "Raster Data" => "rasters.md",