Merging GeoJSON files

There are at least two tools to merge GeoJSON files auto-magically with CLI: mapshaper and @mapbox/geojson-merge. TBH, my favorite one is mapshaper.

Installing mapshaper

C:\Users\USER>npm install -g mapshaper
C:\Users\USER\AppData\Roaming\npm\mapshaper -> C:\Users\USER\AppData\Roaming\npm\node_modules\mapshaper\bin\mapshaper
C:\Users\USER\AppData\Roaming\npm\mapshaper-gui -> C:\Users\USER\AppData\Roaming\npm\node_modules\mapshaper\bin\mapshaper-gui
C:\Users\USER\AppData\Roaming\npm\mapshaper-xl -> C:\Users\USER\AppData\Roaming\npm\node_modules\mapshaper\bin\mapshaper-xl

  • mapshaper@0.4.158
    added 50 packages from 74 contributors in 38s

C:\Users\USER>

Installing @mapbox/geojson-merge

npm install global @mapbox/geojson-merge

npm install -g mapshaper

I used to get E418, an ERROR of 418 when trying to install any package with NPM. Thanks God that I found this solution: change your registry into https:// instead of http://.

registry=https://registry.npmjs.org


log:

Microsoft Windows [Version 10.0.17134.1069]
(c) 2018 Microsoft Corporation. All rights reserved.

C:\Users\USER>npm install -g mapshaper
npm ERR! code E418
npm ERR! 418 I’m a teapot – GET http://registry.npmjs.org/mapshaper – Received invalid host

npm ERR! A complete log of this run can be found in:
npm ERR! C:\Users\USER\AppData\Roaming\npm-cache_logs\2020-02-18T12_29_07_339Z-debug.log

C:\Users\USER>notepad .npmrc

C:\Users\USER>npm install -g mapshaper
C:\Users\USER\AppData\Roaming\npm\mapshaper -> C:\Users\USER\AppData\Roaming\npm\node_modules\mapshaper\bin\mapshaper
C:\Users\USER\AppData\Roaming\npm\mapshaper-gui -> C:\Users\USER\AppData\Roaming\npm\node_modules\mapshaper\bin\mapshaper-gui
C:\Users\USER\AppData\Roaming\npm\mapshaper-xl -> C:\Users\USER\AppData\Roaming\npm\node_modules\mapshaper\bin\mapshaper-xl

  • mapshaper@0.4.158
    added 50 packages from 74 contributors in 38s

RBI 1000

#RBI atau #RupaBumiIndonesia memiliki penomoran yang unik, dan dapat dibeli peta resminya dari #BIG. Namun, terkadang RBI 10 000, yaitu peta rupa bumi berskala 1:10,000, sering belum ada di pemerintah daerah tingkat 2 (kabupaten / kota), dan menunggu dibeli dahulu, padahal citra satelit dan peta #GoogleMap sudah dapat diperoleh dengan mudah di Internet, walaupun bukan resmi dari pemerintah.

Nah, untuk memudahkan nantinya menggabungkan peta dengan RBI, ada baiknya pemerintah daerah juga menyiapkan peta yang nantinya gampang digabungkan / di-overlay dengan peta RBI.

Usulan: RBI 5000, RBI 2500, RBI 1000

Di sinilah saya mengusulkan tiga jenis “cakupan peta RBI” berikutnya: RBI 5000, RBI 2500, dan RBI 1000.

Dibilang CAKUPAN / EXTENTS artinya, bukanlah sebuah peta, namun luas dan pembagian batasnya peta dan penomoran petanya saja.

RBI 5000 dan RBI 2500 dan juga RBI 1000 tentunya tidak official, hanya pendekatan versi saya sendiri, mengikuti pola penomoran dan pembagian wilayah RBI yang sudah ada.

Contoh adalah pada gambar di atas, yaitu proposed RBI 1000 untuk daerah di sekitar Kota Pekanbaru. Terlihat RBI 1000 memiliki 11 angka, dengan contoh sebagai berikut:

0816-5225-311

4 angka pertama adalah RBI 250 000, yaitu: 0816
5 angka pertama adalah RBI 100 000, yaitu: 0816-5
6 angka pertama adalah RBI 50 000, yaitu: 0816-52
7 angka pertama adalah RBI 25 000, yaitu: 0816-522
8 angka pertama adalah RBI 10 000, yaitu: 0816-5225

sedangkan – 3 1 1, tiga angka berikutnya, adalah penambahan saya sendiri, sesuai dengan skala yang sering digunakan pada peta, yaitu skala 5000, 2500 dan 1000.

Jadi, dengan memegang sebuah peta RBI 1000, seseorang bisa mengatakan wilayah tersebut ada pada RBI berapa saja pada skala lebih kecil (2500, 5000, 10000, 25000, 50000, 100000, 250000)

Kenapa RBI?

Karena Peta Rupa Bumi Indonesia adalah dari BIG (Badan Informasi Geospasial) Indonesia, dahulunya Bakosurtanal.

Kartografer seluruh Indonesia sudah terbiasa dengan peta ini.

Kenapa tidak peta batas administrasi?

Peta batas administrasi pemerintahan (misalnya Peta Kabupaten / Peta Kota) memiliki skala yang berbeda-beda tergantung luas wilayah dan ukuran kertanya, sehingga sedikit susah saat di-overlay, ditempelkan diatas peta lain saat rapat koordinasi atau rapat lainnya (non-online) menggunakan printed map, peta kertas.

Memang lebih unggul peta wilayah ini, untuk memetakan wilayah administrasi, karena memang itu tujuan dibuatkan peta tersebut. Namun, saat ada pembangunan jalan, drainase, penanggulangan banjir, dan hal-hal detail lainnya, sangat tidak praktis lagi, karena peta ini biasanya berskala kecil, misalnya 1: 100 ribu atau 1: 250 ribu.

Nah, di sinilah bergunanya pembagian wilayah berdasarkan kode peta RBI ini.

Kenapa tidak peta batas wilayah administrasi pemerintah lebih kecil, misal: peta kelurahan / peta kecamatan?

Terkadang, (sering) kita tidak memiliki petanya. Kalaupun ada, tidak terkini.

Lebih sering, peta tersebut berskala macam-macam, tidak seragam, sehingga saat ditempelkan / di-overlay tidak bisa, karena tidak satu skala, skalanya berbeda.

Di sinilah menangnya peta RBI, karena skalanya sudah pre-determined sebelumnya.

Apa kekurangan penomoran RBI ini?

Tentu ada kekurangannya, yaitu pembagian wilayah hanya berdasarkan posisinya di bumi, yaitu dalam derjat. Sering sebuah wilayah berada di dua atau lebih batas RBI.

Silahkan perhatikan peta administrasi wilayah Kabupaten Kampar dan Kota Pekanbaru berikut:

Terlihat, Kabupaten Kampar terletak pada 4 buah RBI, se-cuil ada di RBI 250 000 nomor peta 0817 dan 0716, sebahagian besar di 0816, dan ada sebahagian kecil di 0815.

Akibatnya, pemerintah terkadang harus membeli beberapa peta RBI untuk wilayahnya, walaupun sebenarnya wilayahnya kecil, namun berada pada beberapa RBI.

Kebetulan Kota Pekanbaru berada sepenuhnya pada RBI 0816.

Kelemahan berikutnya adalah nomornya yang panjang dan susah di-ingat. RBI 250,000 s/d 10,000 (peta resmi dari BIG), biasanya disiasati dengan memberikan nama kota terbesar atau nama most populated place yang termasuk dalam peta tersebut.

Kelemahan lainnya diwariskan oleh sistem koordinat referensi WGS 84 (yang dalam derjat), yang digunakan oleh RBI, yaitu luas wilayah nanti akan berbeda tergantung posisinya dari khatulistiwa. Namun, ini dapat sedikit diakali dengan menggunakan UTM ataupun TM-3.

Apa kelebihan penomoran RBI 1000 ini?

Pemerintah Daerah Tingkat II Kabupaten / Kota tentu akan sangat terbantu dengan pembagian wilayah ini, karena nanti akan sangat gampang digabungkan kembali dengan peta RBI resmi dari pemerintah / BIG. Wilayah yang dipetakan oleh RBI 1000 dengan nomor demikian ada pada peta resmi RBI dengan nomor-nomor sebelumnya.

Karena batas RBI1000 sudah pre-determined berdasarkan aturan RBI yang bersifat nasional, pemda dapat dengan mudah meminta organisasi/perusahaan manapun membuatkannya, cukup dengan memberitahu nomor nya saja.

RBI 1000 juga cocok untuk peta skala 1:1000 pada kertas A2 (59,4 x 42cm), karena ukuran peta tidak sampai 40 x 40 cm.

Terlihat bahwa RBI 10 000, dengan luas cakupan sekitar 2 137 hektar (di Pekanbaru), cocok untuk peta kelurahan,

Proposed RBI 5000 dengan luas cakupan sekitar 534 hektar, cocok untuk RW atau kelurahan yang lebih detil. “RBI 5000” usulan ini merupakan pecahan seperempat dari RBI 10 000. Tiap-tiap RBI 10000 akan memiliki 4 buah RBI 5000. RBI 5000 ini sudah terlihat persil tanah warga, namun belum bisa diberi label . Sila lihat ukuran skalanya berikut:

Catatan:
Pada peta skala 1:100 000, tiap 1 cm pada peta adalah 1 km di dunia nyata.
RBI 10,000, maka 1 cm pada peta adalah 100 meter aslinya,
RBI 5,000, maka 1 cm pada peta adalah 50 meter aslinya,
RBI 2,500, maka 1 cm adalah 25 meter aslinya,
RBI 1,000, maka 1 cm adalah 10 meter aslinya.

RBI 2500, dengan luas cakupan sekitar 136 hektar, cocok untuk peta RW, dan penomoran blok rumah.

RBI 1000, cocok untuk penomoran ruko atau perumahan, karena tanah perumahan di Pekanbaru, biasanya ukurannya adalah 9 x 12 m, jadi di peta masih tergambarkan sekitar 1 cm x 1 cm.

40 x 40 cm?

Ukuran 40×40 cm muat di kertas A2 dan A1 ataupun lebih besar lagi, kertas A0.

Namun, ukuran ini tidak standar dengan kemampuan printer yang ada di pemda, yang biasanya maksimum A3.

Namun, biasanya teknisi kantor pemda sudah tahu cara memprint peta ukuran A0 dengan printer kertas A4 menggunakan metode poster printing.

RBI 1000 dan SAS Planet

RBI 1000 pada SAS Planet untuk ArcGIS Imagery 20 terlihat berukuran 1296 x 1296 pixels.

pada resolusi aslinya, 1296 x 1296 pixels, sudah dapat dilihat dengan jelas rumah dari citra satelit zoom 19. SAS Planet menggunakan zoom dari 1, sedangkan aslinya zoom WebGIS itu dari 0

Tantangan

Tidak semua pemerintah daerah pandai / mampu memecah RBI ini hingga skala sebesar 1000.

Di beberapa tulisan saya, sudah saya buatkan contoh pemecahan RBI ini menggunakan JavaScript + LeafletJS.

Namun, menemukan pegawai Pemda yang mengerti JavaScript (dan mau bekerja untuk ini) juga cukup sulit.

Jika ada pemerintah daerah yang serius ingin memetakan wilayahnya, sebaiknya bekerja sama dengan kampus.

Sebenarnya, saya bisa saja membuatkan RBI 1000 untuk seluruh indonesia, lalu pemda dapat mengunduhnya, namun kita terkendala dengan ukuran file di sini.

D:\GIS\BIG\RBI>dir *indonesia.geojson

28/02/2019  06.28         1.034.654 rbi100Indonesia.geojson
28/02/2019  05.11           164.182 rbi250Indonesia.geojson
28/02/2019  15.27        19.123.853 rbi25Indonesia.geojson
28/02/2019  08.42         4.448.427 rbi50Indonesia.geojson

GeoJSON di atas merupakan polygon batas RBI di Indonesia.

Terlihat RBI 25000 saja sudah berukuran 19MB.

RBI 10000 seharusnya berukuran 9 * 19 MB, sekitar 9 * 20 MB = 180 MB, sudah cukup besar kalau dibuka di QGIS.

Solusi adalah dengan membagi RBI tersebut per provinsi, misal:

D:\GIS\BIG\RBI>dir rbi*riau*

28/02/2019  16.03         5.795.249 rbi10RiauDaratan.geojson
28/02/2019  15.34           408.225 rbi25RiauDaratan.geojson
28/02/2019  16.39        23.527.769 rbi5RiauDaratan.geojson

Alternatif lain, selain membagi per provinsi, yaitu dengan membagi RBI per RBI 250 ribunya.

06/01/2020  19.15        47.442.569 rbi1000_0816-XXXX-XXX.geojson
06/01/2020  19.03         4.842.713 rbi2500_0816-XXXX-XX.geojson
06/01/2020  16.27         1.243.961 rbi5000_0816-XXXX-X.geojson

Terlihat di atas, bahwa RBI 250000 dengan nomor 0816, saat dibuatkan RBI 5000 berukuran 1MB, sedangkan RBI 2500 nya menjadi 5MB, danRBI 1000 nya menjadi 50MB (9 kali-lipat)

Esri Shape File (SHP)

Membuka berkas GeoJSON ini di QGIS untuk ukuran 50MB cukup memberatkan. Solusi adalah dengan meng-convert-nya menjadi SHP dahulu menggunakan OGR2OGR.

ogr2ogr -nlt POLYGON -skipfailures rbi50Indonesia.shp rbi50Indonesia.geojson

perintah di atas adalah untuk mengconvert GeoJSON kita menjadi SHP

Ada pun versi ESRI Shapefilenya:

06/01/2020  19.19        22.021.826 rbi1k_0816-XXXX-XXX.dbf
06/01/2020  19.19               143 rbi1k_0816-XXXX-XXX.prj
06/01/2020  19.20         3.605.488 rbi1k_0816-XXXX-XXX.qix
06/01/2020  19.19        16.920.676 rbi1k_0816-XXXX-XXX.shp
06/01/2020  19.19           995.428 rbi1k_0816-XXXX-XXX.shx

Di suatu sisi, ESRI Shapefile lebih cepat diload, karena bisa dibuatkan berkas QIX nya, yaitu berkas index spatialnya…

ESRI Shapefile juga memisahkan antara database file nya dengan shapenya, PRJ (Projection), SHX (shape index), dan QIX (EXTENTS/tree index)

Beberapa contoh, dapat didownload dari Google Drive saya,

https://drive.google.com/drive/u/1/folders/1sgeQNPltJnPqDzWxbl4FoYbNpqDYgevf

Download VIIRS_SNPP_DayNightBand_ENCC Indonesia PowerPoint Jan 1 2018 – Dec 31 2018

Buat Anda yang ingin melihat langit malam Indonesia hasil satelit Suomi-NPP, dapat mengunduh citra berukuran 1920×1080 yang sudah saya jadikan slides Ms. PowerPoint. Uniknya, kode untuk membuat sudah saya sertakan juga pada VBE nya, jadi Anda bisa melihat langkah-langkahnya, jika Anda mengerti atau mau belajar Visual Basic.

https://drive.google.com/open?id=1J7lOvAz8jh3EFUaGpFR6kXCU_XR15Zps

Ukuran berkas adalah 122 MB

1920 x 1080 pixels

Ukuran 1920×1080 saya pilih karena awalnya adalah saya ingin membuat WallPaper pada Laptop yang resolusi layarnya adalah 16:9. Mendapatkan extent yang pas untuk EPSG:4326 rada-rada triki juga. Setelah beberapa kali uji coba, saya dapatkan Extents yang pas, sebagai berikut:

Driver: JPEG/JPEG JFIF
Size is 1920, 1080

Origin = (86.504122784888892,16.734790559000000)

Pixel Size = (0.031663178490741,-0.031663178490741)

Image Structure Metadata:
COMPRESSION=JPEG
INTERLEAVE=PIXEL
SOURCE_COLOR_SPACE=YCbCr


Corner Coordinates:
Upper Left ( 86.5041228, 16.7347906)
Lower Left ( 86.5041228, -17.4614422)
Upper Right ( 147.2974255, 16.7347906)
Lower Right ( 147.2974255, -17.4614422)

Center ( 116.9007741, -0.3633258)

Yang perlu diperhatikan adalah hanya Corner Coordinates, yaitu EXTENTS batas gambar kita.

Artinya, kalau ingin membuat gambar 16:9, dan ingin memasukkan Indonesia dalam peta tersebut, dapat menggunakan 4 titik koordinat tersebut.

Origin adalah koordinat kiriatas gambar …

Center adalah koordinat tengah gambar …

Query ke NASA

BBOX=-17.4614422,86.5041228,16.7347906,147.2974255

&CRS=EPSG:4326

&LAYERS=VIIRS_SNPP_DayNightBand_ENCC

&FORMAT=image/jpeg

&WIDTH=1920
&HEIGHT=1080

WLD file

Jika sudah dapat EXTENT ataupun BBOX (bounding box), dan dapat ukuran pixel (width dan height), kita dapat membuat World Filenya.

Format WorldFile adalah

baris#1 => resolusi x

(147.2974255 – 86.5041228 )/ 1920

0.03166317848958333

baris#4=> negatif resolusi y

(16.7347906 – -17.4614422) / 1080

ingat, pengurangan dengan negatif adalah penambahan

(16.7347906 + 17.4614422) / 1080

0.03166317851851852

baris#5 dan #6 adalah origin

Origin adalah koordinat kiri atas gambar

Hasil WLD filenya

0.03166317848958333
0
0
-0.03166317851851852
86.5041228
16.7347906 

Download Whole Year

Saya coba download per hari untuk ukuran 1920×1080 ini ke server nasa, dengan BBOX tersebut, dapat ukuran file berkisar 229kB hingga 464kB per hari.

Berhubung ukurannya yang relatif kecil (karena JPEG), saya coba download untuk resolusi lain. Enaknya adalah kalau kita sudah dapat BBOX atau EXTENTS atau CornerCoordinate (all they are related, sama tapi beda format saja), kita bisa mendownload untuk berbagai resolusi, selagi perbandingan nya adalah 16:9 juga …

2560 x 1440 dalah resolusi monitor 27 inch

dan hasil dari NASA masih berukuran 500an kiloByte …

masih kecil …

terkadang lebih baik menghitung dengan python atau JavaScript console, karena karburator eh calc.exe Kalkulator di Windows menggunakan koma, sehingga harus mengganti titik ke koma dulu, lalu hasilnya diganti lagi ke titik. Abaikan Error di atas karena main copy paste saja, ada spasi sebelumnya … hahahaha

3840 x 2160 4K UHDTV

Komputer gamers dengan monitor 4k , ultra high definition TV, adalah 3840 x 2160 pixels.

5120 x 2880 Retina 5K

buat yang hobi makan apel sempong …

7680 x 4320 8K UHDTV Super Hi-Vision

Nah, tadi saya coba download ukuran TV 8K ini … ternyata ukurannya cukup besar, 8 Megaan per hari

Setelah dilihat, nasib kita di Indonesia adalah selalu terlindung oleh awan, saya rasa, terlalu mahal ukuran harddisk segitu …

Gambarnya cukup jelas, tetapi kalau dikalikan 365 hari, saya rasa besar juga ukuran akhir Ms. PowerPoint nya … 8 * 365 adalah sekitar 3 Giga

Makanya saya putuskan pake resolusi layar komputer saya saja … hahaha

Kode VBA – Buat file TXT yang berisi hari dalam setahun

' buat file yang isinya hari dalam tahun 2018
Sub buatfile2018()

    Dim x As Date
    x = DateSerial(2018, 1, 1)
    Open "D:\RS\NASA\VIIRS_SNPP_DayNightBand_ENCC\Indonesia\1920x1080\2018.txt" For Output As #1
    For i = 0 To 365
        Print #1, Format(DateSerial(Year(x), Month(x), Day(x) + i), "yyyy-mm-dd")
    Next i
    Close #1
End Sub

Kode di atas adalah untuk membuat sebuah file 2018.txt yang ber isi semua hari dalam tahun 2018.

Kita butuh file yang berisi tanggal yang valid untuk nantinya diminta ke NASA citra hari tersebut.

for /f %i in (2018.txt) do wget %i

Perintah for ini di CMD.exe berguna untuk melakukan perintah secara repetitif untuk setiap baris file 2018.txt yang sekarang diassign variabel i

WGET saya gunakan karena mendukung -O, yaitu menyimpan dengan nama tertentu.

VBA – Menyiapkan 365 slides

Sub tambahslide()
    Dim objPresentation As Presentation
    Dim objSlide As Slide
    Set objPresentation = ActivePresentation
    Set objSlide = objPresentation.Slides.Item(1)
    For i = 1 To 364
        ActivePresentation.Slides(1).Duplicate
    Next
End Sub

Perintah di atas adalah untuk membuat 365 slide, yang merupakan hasil duplikasi slide pertama nya saja …. saya melakukan perulangan 364 kali menggandakan slide pertama, dan hasilnya saya jadi punya 365 slide.

VBA – Memasukkan gambar satelit otomatis ke PowerPoint

'baca file
Sub baca2018lalubuatkanSlidenya()


    Dim DataLine As String
    
    ' Ms. PowerPoint
    Dim objPresentation As Presentation
    Dim objSlide As Slide
    Dim objImage As Shape
    Dim objText As Shape
    
    Set objPresentation = ActivePresentation

    Open "D:\RS\NASA\VIIRS_SNPP_DayNightBand_ENCC\Indonesia\1920x1080\2018.txt" For Input As #1
    Dim i As Integer
    i = 1
    While Not EOF(1)
        Line Input #1, DataLine ' read in data 1 line at a time
    
        ' decide what to do with dataline,
        ' depending on what processing you need to do for each case
        'Debug.Print DataLine

        Set objSlide = objPresentation.Slides.Item(i)
    
    
        Set objImage = objSlide.Shapes.AddPicture("D:\RS\NASA\VIIRS_SNPP_DayNightBand_ENCC\Indonesia\1920x1080\2018\" & DataLine & ".jpg", msoCTrue, msoCTrue, 0, 0, 1440, 810)
        Set objText = objSlide.Shapes.AddTextbox(msoTextOrientationHorizontal, 0, 0, 200, 50)
        objText.TextFrame.TextRange.Text = DataLine
        objText.TextFrame2.TextRange.Characters.Font.Fill.ForeColor.RGB = RGB(255, 255, 0)
        objText.TextFrame.TextRange.Font.Bold = msoTrue
        i = i + 1
    Wend

    Close #1


End Sub

Wew … jangan takut dulu …

Kode di atas adalah untuk membaca file 2018.txt kita tadi yang berisi 365 hari valid di 2018. Kita butuh nama harinya, karena hasil WGET tadi juga saya buatkan pakai nama hari.

Tinggal diputar, looping WHILE not EndOfFile

Looping ini artinya setiap baris 2018.txt ini, artinya setiap hari kita di tahun 2018

Lalu, kita tambahkan gambar yang berpola hari.jpg ke dalam slide kita, 365 kali..

Nah, agar cantik, kita tambahkan juga tanggal hari nya (TextFrame) … Tanpa teks ini kita tidak tahu nanti tanggal berapa citra satelitnya …

Tadi awalnya ngga pake kolor … ngga pake kolor default nya anunya item … dan kadang ngga nampak tulisannya

Makanya saya tambahkan kolor ijo … eh RGB Red+Green = kuning

Kolor kuning bukan karena saya anak UI ya … tapi karena biar lebih jelas saja untuk malam hari… lagian ngasal aja tadi itu

Voila … 365 gambar satelit setahun masuk secara otomatis ke slide kita, berikut dengan tanggal harinya …

PPTM – PPTX Macro Enabled

Berhubung kemungkinan Anda yang baca hingga di sini ngga banyak, dan hanya mereka yang belum muntah dengan kode, dan ingin belajar …

SELAMAT!!!

kode VBA tadi sudah saya masukkan ke dalam PowerPoint nya, dan karena ada kode VBA nya, harus Macro Enabled PPTX … atau lebih dikenal dengan PPTM

Anda dapat melihat kode saya ini dengan Alt+F11 pada PowerPoint nya

Hati hati dengan PPTM

File PPTM bisa saja disusupi kode jahat, ya … jadi jangan langsung di enabled, tapi Alt+F11 dulu lalu pelajari kodenya. Kalau kodenya tidak jahat, baru di enable …

Hati-hati dengan macro / VBA, karena Anda bisa saja menghapus atau merusak data Anda dengannya … saya pernah kehilangan data berharga hanya karena bikin makro … lalu di eksekusi oleh Abang saya tanpa sengaja karena saya yang salah, ter assign pula Ctrl+S … yang pastinya biasanya adalah untuk menyimpan ….

Namun, dibalik kekuatannya, kalau Anda mampu menggunakannya dengan baik, akan sangat menbantu pekerjaan. Bayangkan kalau saya harus mengetik nama hari setahun, 365 baris, tentu dengan sedikit kodingan menjadi lebih simpel. (lihat kode kita di atas, membuat file TXT)

Begitupun dengan mendownload dan menyimpan , bisa diotomatisasi pakai perintah FOR

Membuat slides 365 kali sebenarnya cukup enter 364 kali …

Nah disinilah kelebihan VBA, bisa lebih presisi … capek juga kan ngitung enter 364 kali …

VBA juga terlihat lebih presisi dalam meletakkan gambar. Awalnya, saya harus menggeser geser agar benar-benar pas gambarnya di koordinat 0,0 dengan mouse.

Dengan VBA, otomatis …

Dan VBA melakukannya 365 kali secara otomatis …

Download NASA/NOAA Suomi-National Polar Orbiting Partnership (S-NPP) Satellite – Visible Infrared Imaging Radiometer Suite (VIIRS) Corrected Reflectance (True Color)

GIBS Nasa memberikan layanan percuma untuk mengunduh hampir-seketika (near-realtime) data hasil koreksi reflektan dari instrumen VIIRS dari satelit Suomi. Data yang puluhan Giga disarikan menjadi layanan ubin pada WorldView-nya EarthData -nya NASA.

Nah, karena sudah berbentuk layanan ubin, kita bisa unduh otomatis, dan kalau Anda memiliki WGET.EXE dan juga sudah menginstal ImageMagick, Anda dapat menggunakan script batch file saya ini untuk mendownload eh meng-unduh (bahasa Indonesia yang baku untuk download adalah unduh) semua ubin, lalu di-montage (saya tidak tahu padanan kata montage dalam Bahasa Indonesia)

POC (Proof of Concept)

Sila lihat skrip saya di sini:

https://drive.google.com/drive/folders/14iNuMAgYXbtb4q99HIRpgnNv0H4CK-Wa?usp=sharing

Pastikan Anda Memiliki WGET.exe

karena batch file saya mendownload ubin menggunakan perintah wget

Pastikan Anda memiliki ImageMagick

karena batch file saya dalam menggabungkan ubin menjadi citra penuh menggunakan perintah montage dari ImageMagick

Cara Penggunaan

buat sebuah folder dengan format yyyy-mm-dd, misalnya 2019-12-01, lalu pindahkan file *.cmd dan *.txt ke folder tersebut.

Empat buah file wget*.cmd adalah mendownload secara paralel secara terpisah dari masing-masing subdomain GIBS, jadi biar lebih cepat saja downloadnya.

Eksekusi wget3_0.cmd, lebih baik pakai tombol keyboard enter daripada double click, karena nanti bisa kegeser. Lalu pilih wget3_a.cmd lalu enter lagi. Pilih wget3_b.cmd lalu enter lagi. Terakhir pilih wget3_c.cmd lalu tekan tombol enter juga di keyboard.

Jika semua ubin telah terdownload, dengan ditandai anda telah memiliki 50 buah ubin dengan format matriks-colom-baris.jpg,

Baru eksekusi (pilih, lalu enter) file domontage3.cmd

Hasil

Anda akan otomatis membuat file berformat m3_yyyy-mm-dd.jpg, yang merupakan hasil penggabungan 50 buah ubin tadi. Ukuran citra akhir seharusnya 5120 x 2560 pixels.

Penjelasan Sedikit tentang Matrix Set 3

Zoom atau Matrix set 3, sudah berukuran 5120 x 2560 pixels, jadi sudah sangat jelas kalau hanya untuk bahan ajar. Untuk penelitian, mungkin Anda harus download hingga Matrix 7. Tapi ingat, setiap penambahan level matriks, ukuran menjadi 4 kali lipat lebih banyak. m4 adalah 10240 x 5120 pixels. m5 adalah 20512 x 10240 pixels, begitu seterusnya …

Banyak unduhan ubin juga berkali 4 tiap penambahan level matrix. Matrix 4 berukuran 5 baris 10 kolom, yaitu 50 ubin.

Matrix 5 akan berukuran 10 baris 20 kolom, yaitu 200 ubin.

Matrix 6 akan berukuran 20 baris dan 40 kolom, yaitu 800 ubin.

Matrix 7 akan berukuan 40 baris dan 80 kolom, yaitu 3200 ubin.

Masing-masing ubin berukuran 512×512 pixels.

Dan ingat, walaupun satelit ini melintasi bumi berkali-kali, namun ada beberapa daerah terutama semakin dekat ke khatulistiwa yang tidak tercitrakan, karena lebar swath satelit tidak cukup berhimpitan, dan sepertinya diisi otomatis oleh algoritma NASA, dan tentu kita tertolong dengan ini … Tapi algoritma tersebut belum 100% seprulna, sehingga masih kelihatan sedikit blur.

Ada lagi sedikit, tentang gambar yang tidak smooth tersambung, karena sebenarnya itu adalah hasil blending dari waktu yang tidak sama. Terlihat awan terpotong atau bahkan cahaya nya berbeda.

Nah, satu lagi yang sedikit mengganggu adalah bagian kutub tertentu yang tidak terpetakan, karena memang menjauhi matahari dan seperti kebanyakan satelit, orbitnya sinkron dengan matahari. Bagian bumi yang menjauhi matahari jadi susah untuk dilihat…

Going Deeper …

Jika Anda sukses mengikuti semua hingga sini, WOW, Anda sudah jago dalam mengunduh citra satelit dari WorldView nya NASA.

Congratulate Yourself, Mate! ngga banyak orang yang sudah sampai ke sini …

Sedikit permasalahan adalah kalau Anda ingin mengolah lagi gambar tersebut di QGIS. Anda butuh World File, yaitu file berekstensi WLD ataupun JGW ataupun PGW jika file Anda berakhiran JPG atau PNG respektifli.

Cara cepat adalah tinggal download saja file WLD nya dari Google Drive saya, karena EXTENTS dan WIDTH dan HEIGHT gambar sama, makanya WLD nya juga akan pasti sama.

Biar Anda tidak ditanya mulu CRS nya apa (kode EPSG nya) oleh QGIS, makanya kita perlu juga buatkan berkas AUX. Berkas tersebut berisi kode Coordinate Reference Systems. Terkadang, karena NASA sudah menggunakan WGS 84 (EPSG:4326), otomatis default EPSG adalah ini.

Evolusi Bumi

Nah, ada yang menarik dari citra hasil satelit Suomi ini, karena jatuhnya yang mengelilingi bumi, posisinya sedikit bergeser dari hari-ke-hari. Untuk melihat lebih kontras, saya demokan di GoogleDrive setiap 3 bulan. Sila lihat hasil unduh bulan 3, bulan 6, bulan 9 dan bulan 12.

Terdapat ada wilayah bumi yang tidak terpetakan oleh satelit.

Desember (desi=10, ember=bulan), Matahari di Selatan

Pada bulan September ke puncaknya bulan Desember, hingga Maret nantinya, maka kutub Utara akan tidak terpetakan, karena utara bumi menjauhi matahari pada fase itu. Bagian kutub utara bumi akan tidak akan mendapati siang hari, akan terus malam hari selama beberapa hari.

Di beberapa tempat di bagian utara bumi akan lebih lama malam daripada siangnya. Jika Anda muslim, dan dapat waktu Ramadhan saat bulan Desember, maka Anda akan sahur, imsak, lalu sholat subuh, lalu sholat zuhur, lalu sholat ashar lalu berbuka dan sholat Maghrib, serasa sholat subuh 2 + 4 + 4 + 3 rakaat … Sholat Jum’at nya ngga boleh lama khutbah, takutnya sudah Ashar bahkan Maghrib …

Bumi akan berada pada posisi “tegak lurus” pada Maret dan September. Pada saat itu, terjadi “equinox”, yaitu di mana semua tempat di bumi sama panjang siang dan malamnya.

Juni, Matahari di Utara

Pada bulan Juni, matahari akan berada di bagian utara bumi, sehingga kutub Utara akan siang selama beberapa hari tanpa malam. Di tempat lain sekitarnya akan lebih lama siang daripada malam.

Sebenarnya, bagian utara bumi mendekati matahari, dengan teleng 23,5 derjat. Begitupun dengan bagian selatan bumi yang menjauhi matahari.

Akibatnya, bagian selatan juga tidak terpetakan oleh Saomi …

[QGIS] [Layout] Rounded Corner LINE

Terkadang, kita ingin membuat garis batas peta kita 1 cm dari ujung kertas. Nah, untuk membuat garis 1 cm dari ujung masing-masing kertas tentu kita harus mengetahui dulu ukuran kertasnya, lalu dibuatkan kotak 1 cm dari atas kiri, dengan panjang kotak adalah ukuran kertas dikurang 2 cm (karena digunakan sebagai jarak 1 cm di kiri dan kanan), dan begitu pun tinggi kotak adalah tinggi kertas dikurangi 2 cm, karena 1 cm untuk batas atas dan bawah.

Contoh: Kertas A2

Ukuran kertas A2 adalah 594 x 420 mm. Ukuran kotak kita adalah (594-20) x (420 – 20) mm, dengan posisi awal kotak adalah 10 mm atas, 10 mm bawah.

Kok malah 10 mm? Karena 1 centimeter adalah 10 milimeter!

Orang Indonesia biasa dengan ukuran senti, jarang menggunakan milimeter, kecuali beberapa kalangan anak teknik!

Rounded Corner

Untuk membuat kotak dengan ujung tumpul bulat, atau istilah sumaterabarat nya rounded corner, pada QGIS cukup ditambahkan sebuah kotak persegi panjang, dengan sudut radius diubah menjadi 5 milimeter, seperti terlihat pada gambar.

Seperti contoh kita di atas, kita posisikan X dengan angka 10mm, dan Y dengan angka 10 mm, sehingga kotak kita posisi awalnya adalah 10 milimeter dari atas kiri kertas.

Nah, sekarang yang jadi masalah adalah Lebar dan Tinggi. Karena kertas kita nantinya adalah A2, dengan lebar kertas adalah 594 mm, maka lebar kotak persegi panjang kita adalah 574 mm.

Tinggi kertas A2 kita adalah 420mm, makanya tinggi kotak persegi panjang kita adalah 400 mm!

Gampang sekali, kan …!

Double Rounded-Corner Line

Pada contoh gambar di atas, kalau dilihat dengan teliti, terlihat sebenarnya sudah ada 2 buah kotak persegi panjang, dengan garis pertama (tebal) adalah 10 mm pada kertas, dan garis kedua (tipis) berjarak 9 mm dari ujung kertas, sehingga menjadi garis terluar.

Silahkan lihat ukuran kotak kita sedikit lebih aneh, karena jarak dari ujung kertas adalah 9 mm, makanya berdasarkan persamaan berikut:

lebar kotak = lebar kertas – (2 * jarakdarikertas)

lebar kotak adalah =594 – ( 2 * 9) = 594 – 18 = 576 mm

tinggi kotak = tinggi kertas – (2 * jarakdarikertas)

tinggi kotak = 420 – (2 * 9) = 420 – 18 = 402 mm

Radius

Radius kotak kedua juga harus diubah, tidak lah 5 mm lagi, tetapi 6 mm, biar mecing dengan kotak pertama.

~

VBA menambahkan citra hasil satelit otomatis ke slides powerpoint

ternyata rempong juga ya memasukkan satu per satu citra hasil download satelit ke slides nya Ms. PowerPoint … Alhamdulillah nemu VBA nya …

Sila dilihat sampah VBA saya di sini …

Jangan langsung dicopy-paste ya … pelajari dahulu …

Kelebihan metode ini adalah PRESISI ukuran gambar, karena gambar biasanya saya POSISIKAN lagi (dari drag & drop) … dan tentu ini agak rempong …

Dengan kodingan, kita bisa langsung posisikan dengan presisi

Using Python ImageIO to generate GIF file from tiMe lapsed images

NASA WorldView way too many frames were selected…

If you love remote sensing, you’ll probably in some position willing to create an animated GIF file from some of your rasters. Well, there area plenty of options in on-the-shelf commercial software, but code it directly with python is some how a great accomplishment — at least for me –.

You’ll need to install ImageIO first to generate a GIF file, simply by typing pip install imageio


D:\RS\NASA\Himawari\Riau>pip install imageio
Collecting imageio
  Downloading https://files.pythonhosted.org/packages/1a/de/f7f985018f462ceeffada7f6e609919fbcc934acd9301929cba14bc2c24a/imageio-2.6.1-py3-none-any.whl (3.3MB)
     |¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦| 3.3MB 297kB/s
Requirement already satisfied: pillow in d:\python36\lib\site-packages (from imageio) (6.2.0)
Collecting numpy
  Downloading https://files.pythonhosted.org/packages/b0/ee/5ff445dd43b9820e5494d21240e689d3b7cb52bc93f4f164eba84206cd0d/numpy-1.17.4-cp36-cp36m-win_amd64.whl (12.7MB)
     |¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦| 12.7MB 327kB/s
Installing collected packages: numpy, imageio
Successfully installed imageio-2.6.1 numpy-1.17.4

then, create a simple python script like this … please change filenames array into your own files


#####=====#####=====#####=====#####=====#####=====#####=====
# file: imageio_creategif.py
#####=====#####=====#####=====#####=====#####=====#####=====
import imageio

#filenames = ['riau_20191201_0000.png','riau_20191201_0010.png','riau_20191201_0020.png']
# 0240 dan 1440 rusak
filenames = [
  'riau_20191130_2200.png','riau_20191130_2210.png','riau_20191130_2220.png','riau_20191130_2230.png','riau_20191130_2240.png','riau_20191130_2250.png','riau_20191130_2300.png','riau_20191130_2310.png','riau_20191130_2320.png','riau_20191130_2330.png','riau_20191130_2340.png','riau_20191130_2350.png'
  ,'riau_20191201_0000.png','riau_20191201_0010.png','riau_20191201_0020.png','riau_20191201_0030.png','riau_20191201_0040.png','riau_20191201_0050.png','riau_20191201_0100.png','riau_20191201_0110.png','riau_20191201_0120.png','riau_20191201_0130.png','riau_20191201_0140.png','riau_20191201_0150.png','riau_20191201_0200.png','riau_20191201_0210.png','riau_20191201_0220.png','riau_20191201_0230.png','riau_20191201_0250.png','riau_20191201_0300.png','riau_20191201_0310.png','riau_20191201_0320.png','riau_20191201_0330.png','riau_20191201_0340.png','riau_20191201_0350.png','riau_20191201_0400.png','riau_20191201_0410.png','riau_20191201_0420.png','riau_20191201_0430.png','riau_20191201_0440.png','riau_20191201_0450.png','riau_20191201_0500.png','riau_20191201_0510.png','riau_20191201_0520.png','riau_20191201_0530.png','riau_20191201_0540.png','riau_20191201_0550.png','riau_20191201_0600.png','riau_20191201_0610.png','riau_20191201_0620.png','riau_20191201_0630.png','riau_20191201_0640.png','riau_20191201_0650.png','riau_20191201_0700.png','riau_20191201_0710.png','riau_20191201_0720.png','riau_20191201_0730.png','riau_20191201_0740.png','riau_20191201_0750.png','riau_20191201_0800.png','riau_20191201_0810.png','riau_20191201_0820.png','riau_20191201_0830.png','riau_20191201_0840.png','riau_20191201_0850.png','riau_20191201_0900.png','riau_20191201_0910.png','riau_20191201_0920.png','riau_20191201_0930.png','riau_20191201_0940.png','riau_20191201_0950.png','riau_20191201_1000.png','riau_20191201_1010.png','riau_20191201_1020.png','riau_20191201_1030.png','riau_20191201_1040.png','riau_20191201_1050.png','riau_20191201_1100.png','riau_20191201_1110.png','riau_20191201_1120.png','riau_20191201_1130.png','riau_20191201_1140.png',
]
with imageio.get_writer('riau_20191201.gif', mode='I') as writer:
    for filename in filenames:
        image = imageio.imread(filename)
        writer.append_data(image)
        
#####=====#####=====#####=====#####=====#####=====#####=====

please being prepared to create a huge GIF file …

that previous python script is my real python file to generate this 88MB gif of an animated full day Himawari-8 scene on December 1st, 2019, local time … well Riau province is GMT+7 or UTC+7 … that’s why I need to also using previous day images. Watch thoroughly the minutes and the date in the files’ name!

Google Drive links for my messy Himawari-8 Riau

https://drive.google.com/open?id=12ifqMMogUntZezrm0kQXEAGK3LhANrmN

or if you prefer Indonesia:

https://drive.google.com/drive/u/1/folders/1oSKOLJtxuuiEC1gXgWOK49HC6luWoaqc

or directly download this 167MB Indonesian GIF

https://drive.google.com/open?id=1bYDK6QWuyx0HG8LPtTteswJ9mTIMIm_N

JavaScript untuk membuat WLD dan AUX.XML file NASA GIBS

GIBS NASA menggunakan EPSG:4326, yang sedikit berbeda dengan Pseudo Mercator / EPSG:3857 atau yang lebih dikenal dengan Google Map. Zoom pada GIBS pun dinamakan Matrixset m, jadi bukan z.

Ukuran gambar dibuat dengan perbandingan 2:1, untuk lebar dua kali lipat tinggi.

Ukuran tile pun 512×512 pixels, tidak seperti kebanyakan tile server lain yang 256×256 pixels.

Gambar baru “benar-benar” full saat m = 3, dimana terdapat 5 baris dan 10 kolom

Leaflet JS membuat Polygon untuk Matrik 3

Untuk membuat polygon batas masing-masing tile tersebut, kita gunakan saja Leaflet, karena dukungan GIS nya memang sangat mumpuni. Berikut petikan awal snippet JavaScript saya:

// fungsi untuk membuat polygon berdasarkan nomor row dan col
// untuk matrik 3
// matrik 3 adalah 10 col dan 5 row
function m3(row,col){
  //ingat L.polygon adalah latlong
  return new L.polygon(
    [
      [ 90-36*row , -180 + 36*col],
      [ 90-36*row , -180 + 36*(col+1)],
      [ 90-36*(row+1) , -180 + 36*(col+1)],
      [ 90-36*(row+1) , -180 + 36*(col)]
    ]
  );
}

Fungsi di atas untuk membuat sebuah polygon di Leflet, untuk spesifik nomor baris dan kolom tertentu.

Perhatikan bahwa Leaflet menggunakan LatLng, latitude dulu baru , baru longitude!

Untuk zoom level 3, kita tahu ada 10 kolom gambar untuk 360 derjat longitude. Dari sini, kita bisa tahu bahwa lebar sebuah ubin adalah 36 derjat.

Untuk zoom level 3, kita tahu ada 5 baris gambar yang menutupi bumi 90 derjat utara dan 90 derjat selatan, total 180 derjat latitude bumi.

180 dibagi 5 adalah …. 36 derjat juga …

Extents (batas) gambar dari NASA adalah seluruh bumi, dari -180 sampai 180 untuk longitude, dan dari 90 utara hingga minus 90 untuk selatan…

Longitude bertambah terus 36 derjat setiap kolom ubinnya, dari -180 (ujung barat bumi) hingga 180 (timur).

Latitude berkurang terus 36 derjat untuk penambahan baris ubin, mulai positif 90 di kutub utara, hingga minus 90 untuk kutub selatan.

Sekarang Anda tahu darimana datang angka-angka koordinat polygon nya!

LeafletJS Download GeoJSON untuk semua Matrik 3

Silahkan pelajari petikan kode berikut:

function downloadm3(){
	// Create an empty GeoJSON collection
	var collection = {
		"type": "FeatureCollection",
		"features": []
	};
	
	//sekarang loop tiap row dan col
  // m3 ada 5  baris 10 kolom
  var row=5,col=10;
	for(r=0;r<row;r++){
		for(c=0;c<col;c++){
			g = m3(r,c).toGeoJSON();//dari leaflet polygon langsung ke GeoJSON saja
      //nah, berhubung properties nya masih kosong, kita tambahin saja
      g.properties={m:3,r:r,c:c};
			collection.features.push(g);
		}
	}
	
	saveAs(new Blob([JSON.stringify(collection)], {type: "text/plain;charset=utf-8"}),'m3.geojson')
	
}  

Nah, jika fungsi tersebut dipanggil, Anda akan disuruh unduh eh download sebuah file bernama m3.geojson, yaitu geojson untuk semua batas ubin GIBS nasa untuk matriks 3.

Untuk sukses menjalankannya, Anda butuh FileSaver.min.js

Logika dibelakang kode ini adalah seperti berikut:

Kita buatkan dulu sebuah koleksi untuk menampung fitur kita. Ini merupakan aturan pembuatan FeatureCollection pada GeoJSON.

Lalu, kita push masing-masing ubin kita ke koleksi tadi.

Untuk lebih memudahkan kita identifikasi nomor baris dan kolom, kita tambahkan tabel atribut untuk GeoJSON kita. Pada GeoJSON, tabel atribut namanya adalah PROPERTIES!

Untuk memudahkan mendownload file GeoJSON, saya gunakan FileSaver JS, sehingga dengan fungsi saveAs nya kita bisa langsung bisa simpan berkas dari memori javascript kita.

Lanjut …

Setelah sukses untuk GeoJSON zoom level 3, kita coba untuk membuat zoom level berikutnya.

//m3 adalah 10x5, col nya adalah 360 / 10, masing masing 36 derjat
//m4 adalah 20x10, col nya adalah 360/ 20, masing-masing 18 derjat
//m5 adalah 40x20, col nya adalah 360/ 40, masing-masing berarti 9 deg
// m6  80x40     360/80 = 9/2 deg = 4.5 deg
// 7  160    36/16   36/2^4
function mx(m,row,col){
  //ingat L.polygon adalah latlong
  var deg = 36/2**(m-3);
  return new L.polygon(
    [
      [ 90-deg*row , -180 + deg*col],
      [ 90-deg*row , -180 + deg*(col+1)],
      [ 90-deg*(row+1) , -180 + deg*(col+1)],
      [ 90-deg*(row+1) , -180 + deg*(col)]
    ]
  );

}

NASA GIBS mengikuti gaya tile server pada umumnya, yaitu pada zoom level berikutnya eh matriks berikutnya, adalah matriks sekarang dipecah 2 sehingga menjadi 4 ubin,

m3 adalah 10 kolom x 5 baris.

Seperti dibahas tadi, masing-masing panjang dan tinggi nya adalah 36 derjat.

Untuk zoom berikut, dibagi 2,

m4 berisi 20 kolom dan 10 baris, karena merupakan hasil pembagian dua dari m3.

Masing-masing ubin di m4 jadinya berukuran 18 derjat x 18 derjat.

Dapatkah Anda melihat pola berikut?

m3 -> 36

m4 ->18

m5 -> 9

Berapakah derjat untuk angka m berikutnya?

var deg = 36/2**(m-3);

Jawabannya adalah petikan kode di atas …

dua buah asterik pada javascript adalah pangkat.

2**3 = 2 ^ 3 = 8

perhatikan tanda pangkat di JavaScript

Download GeoJSON untuk zoom berapa saja

Setelah kita tahu cara membuat polygon setiap zoom berapa saja , saatnya kita buatkan fungsi untuk mendownload GeoJSON untuk zoom berapa saja.

Kodenya sebagai berikut:


/**
 * fungsi untuk mendownload spesifik m atau zoom level NASA
*/
function dlm(m){
  // Create an empty GeoJSON collection
  var collection = {
    "type": "FeatureCollection",
    "features": []
  };
  
  //sekarang loop tiap row dan col
  // m3 ada 5  baris 10 kolom
  //var row=5,col=10;
  // m4 20
  var col=2**(m-3) * 10;
  var row=col/2;
  
  
  for(r=0;r<row;r++){
    for(c=0;c<col;c++){
      g = mx(m,r,c).toGeoJSON();//dari leaflet polygon langsung ke GeoJSON saja
      //nah, berhubung properties nya masih kosong, kita tambahin saja
      g.properties={m:m,r:r,c:c};
      collection.features.push(g);
    }
  }
  
  saveAs(new Blob([JSON.stringify(collection)], {type: "text/plain;charset=utf-8"})
    ,'m'+m+'.geojson')
  
}    
    

Jika Anda menjalankan dlm(8), maka Anda akan disuguhkan download file m8.geojson, yang berisi SEMUA polygon untuk zoom 8.

Masih belum muntah?

Membuat file WLD

Kelebihan kita telah memiliki polygon untuk setiap ubin, kita bisa tahu extents atau bound di Leaflet.

Perhatikan kode konsol JavaScript berikut:

var mx311 = mx(3,1,1)

JSON.stringify(mx311.getBounds())
"{"_southWest":{"lat":18,"lng":-144},"_northEast":{"lat":54,"lng":-108}}"

Setiap fitur polygon di Leaflet memiliki fungsi getBounds(), yang sangat berguna bagi kita untuk mengetahui EXTENTS gambar ubin kita.

Nah, kalau kita sudah memiliki EXTENTS, kita bisa buatkan World Filenya

function bound2WLD(bound){
  //baris pertama adalah resolusi x
  var a=[];
  a.push((bound._northEast.lng-bound._southWest.lng)/512);
  a.push(0);
  a.push(0);
  a.push( -1* (bound._northEast.lat - bound._southWest.lat)/512);
  a.push(bound._southWest.lng);
  a.push(bound._northEast.lat);
  return a;

} 

Nah, berhubung citra ubin kita adalah 512 x 512 pixel, makanya resolusi adalah dibagi 512.

Kesimpulan sementara

Polygon -> bound -> WLD

Download WLD otomatis

untuk sementara, kita coba download semua file WLD untuk zoom level 3 saja dahulu … nanti kalau bener, baru kita download semua …

//download semua file WLD untuk zoom 3
function wld3zip(){
  //coba pake JSZIP
  var zip = new JSZip();
  //untuk zoom 3, row adalah 5, cols adalah 10
  var row=5,col=10;
  for(r=0;r<row;r++){
    for(c=0;c<col;c++){
      zip.file((3+'-'+r+'-'+c+'.wld') , bound2WLD(mx(3,r,c).getBounds()).join('\r\n'));
    }
  }
  
  zip.generateAsync({type:"blob"})
    .then(function (blob) {
    saveAs(blob, "m3.zip");
  });
	
}

Nah, berhubung kita mendownlod banyak file, sebaiknya dizip saja, menggunakan JSZIP. Hal ini memudahkan kita mendownload 50 file, karena dijadikan 1 saja, m3.zip.

Download WLD untuk zoom berapa saja

Setelah sukses download semua WLD untuk zoom level 3, kita lanjutkan ke level berapa saja ….

Kodenya mirip dengan sebelumnya, hanya ada beberapa yang perlu kita ubah …


//fungsi untuk mendownload semua file WLD untuk zoom level m
// @input m zoom
// @output ZIP file containing all WLD files
function wldzip(m){
  // pake JSZIP
  var zip = new JSZip();
  var col=2**(m-3) * 10;
  var row=col/2;
	for(r=0;r<row;r++){
		for(c=0;c<col;c++){
			zip.file((m+'-'+r+'-'+c+'.wld') , bound2WLD(mx(m,r,c).getBounds()).join('\r\n'));
		}
	}
  
  zip.generateAsync({type:"blob"})
    .then(function (blob) {
    saveAs(blob, "m"+m+".zip");
  });
	
}

Perhatikan bahwa kita perlu mengubah juga nama file WLD nya ya dan juga nama zip download nya, dan juga mx(3,r,c) sekarang menjadi mx(m,r,c).

Download ZIP file PNG.AUX.XML

Berikut kode untuk membuat file auxiliary dari semua ubin yang kita unduh. Kita butuh file png.aux.xml ini agar tidak perlu ditanya terus CRS / Kode EPSG untuk gambar kita.


//download zip file containing all png.aux.xml 
function zipAUX(m){
  var zip = new JSZip();
  var col=2**(m-3) * 10;
  var row=col/2;
	for(r=0;r<row;r++){
		for(c=0;c<col;c++){
			zip.file(
        (m+'-'+r+'-'+c+'.png.aux.xml') //filename
        , '<PAMDataset>\r\n' 
          +'  <SRS>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"]]</SRS>\r\n'
          + '</PAMDataset>'
      );
		}
	}
  
  zip.generateAsync({type:"blob"})
    .then(function (blob) {
    saveAs(blob, "PNG_AUX_XML_m"+m+".zip");
  });

}

Logikanya mirip dengan sebelumnya …

Voila .. sekarang kita sudah punya semuanya … ubin (download via wget), WLD dan juga AUX nya

Sekarang, gambar ubin apapun dari NASA, bisa kita unduh dan tampil cantik di QGIS kita!

Mabok?

temui saya di Gedung Baru FST UIN Suska … kita bahas secaran offline, biasanya lebih nancap!

GIS Lecturer, in Bahasa (Indonesian)

%d blogger menyukai ini: