texts:
I want to create a color-relief and a hillshade with gdaldem and use them with mapnik. System: Ubuntu 14.04 LTS Server, gdal-bin 1.10.

Here are my steps:

Download all SRTM tiles from viewfinderpanoramas:

wget -r -np -l 1 -A zip http://viewfinderpanoramas.org/Coverage%20map%20viewfinderpanoramas_org3.htm
Unpack all zip files
Fill all remaining voids (creates out of every .hgt a .hgt.tif file)

for hgtfile in *.hgt;do gdal_fillnodata.py $hgtfile $hgtfile.tif; done
Merge all .tifs into one huge tif. This file will is the raw DEM with full resolution and the start for any further steps:

gdal_merge.py -n 32767 -co BIGTIFF=YES -co TILED=YES -co COMPRESS=LZW -co PREDICTOR=2 -o ../raw.tif *.hgt.tif
Convert the raw file into Mercator projection, interpolate and shrink e.g. to 1 km/pixel:

gdalwarp -co BIGTIFF=YES -co TILED=YES -co COMPRESS=LZW -co PREDICTOR=2 -t_srs "+proj=merc +ellps=sphere +R=6378137 +a=6378137 +units=m" -r bilinear -tr 1000 1000 raw.tif warp-1000.tif
Create color-relief with transparency (alpha channel) and lzw compression:

gdaldem color-relief -co COMPRESS=LZW -co PREDICTOR=2 -alpha warp-1000.tif relief.txt relief-1000.tif
Create hillshade

gdaldem hillshade -z 7 -compute_edges -co TILED=YES -co COMPRESS=JPEG warp-1000.tif hillshade-1000.tif

I get these output file sizes:

warp-5000.tif: 113 MB
warp-1000.tif: 606 MB
warp-500.tif: 1.7 GB

relief-5000.tif: 17 MB
relief-1000.tif: 19 MB
relief-500.tif: >60 GB (hdd was full!)

hillshade-5000.tif: 2.7 MB
hillshade-1000.tif: 75 MB
hillshade-500.tif: >70 GB (hdd was full!)

Why does the file size suddenly explode?


texts: Randomly sample existing points from within polygons in SF | R

texts:
I have a set of existing sf points and want to randomly sample 1 point from each of a 1km2 grid cell over a landscape.

I can create a grid of 1km2 polygons, but I can't work out in sf how to then sample points within these polygons.

E.g.

library("sf")
polys <- st_read(system.file("shape/nc.shp", package="sf"))
polys <- polys[1:3,]
polys <- st_geometry(polys) %>% st_cast("POLYGON")

pts <- st_sample(polys, 30)

plot(polys)
plot(pts, add = TRUE)

Now I want to extract one element of pts for each polygon in polys, with the result looking like:

Is there a function to do this within the sf framework?


texts: GDAL driver.Create() TypeError for datatype argument

texts:
I have used the osgeo.gdal module to save numpy arrays as GeoTIFF files successfully in python for some time. Today I decided to write a simple module to handle the driver and file creation routines. Calling my module to save a numpy array gives the following error:

Traceback (most recent call last):
File "saveRaster.py", line 30, in <module>
    save_raster(destination,1,array,srs,gt)
File "saveRaster.py", line 10, in save_raster
    dataset_out = driver.Create(path, cols, rows, bands, dtype)
File "/Library/Frameworks/GDAL.framework/Versions/2.1/Python/2.7/site-packages/osgeo/gdal.py", line 1440, in Create
    return _gdal.Driver_Create(self, *args, **kwargs)
TypeError: in method 'Driver_Create', argument 5 of type 'int'

This is odd since argument 5 is the datatype argument, which should take a gdal data type such as gdal.GDT_Float32. However, if it try type(gdal.GDT_Float32) it returns <type 'int'>. Why then is gdal giving me a TypeError if an integer ought to be exactly what it is expecting?

Code below:

#! /usr/bin/env python

from osgeo import gdal,osr
import numpy as np

def save_raster(path, band_count, bands, srs, gt, format='GTiff', dtype=gdal.GDT_Float32):
    cols,rows = bands.shape
    # Initialize driver & create file
    driver = gdal.GetDriverByName(format)
    dataset_out = driver.Create(path, cols, rows, bands, dtype)
    dataset_out.SetGeoTransform(gt)
    dataset_out.SetProjection(srs)
    # Write the array to raster bands
    for b in range(band_count):
        out_band = dataset_out.GetRasterBand(b+1)
        out_band.WriteArray(bands[b])
    # Write file to disk
    dataset_out.FlushCache()

gt = [0,1,0,0,0,-1]

srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
srs = srs.ExportToWkt()

destination = '~/Desktop/arr.tif'

array = np.arange(0,25).reshape(5,5)

save_raster(destination,1,array,srs,gt)


texts: Choosing table structure for MSSQL Server to render GeoSpatial Data on OpenLayers Map?

texts:
I have billions of public data sets.
I created MSSQL Database with the following table Structure:

ID:nvarchar,
Name:nvarchar,
Latitude:nvarchar,
Longitude:nvarchar,
Geometry:geometry

To represent on Openlayer map I'm going to store data and to render storing data on Map for that I'm going to use GeoServer MSSQL DataStore.

What should my table structure be like so that I will not get any performance issues?


texts: Check field type and execute depending on it using ArcPy?

texts:
I had a script that converted data from a clients site to a usable GIS format until they changed the output.  I've tried two different directions with no luck.  Due to how an earlier part of the script works, CSV to Table, and the nature of the data being converted; the 'total_length' field that's imported might end up being a string or an integer (TableToTable conversion will not work).  I need it to look at the field type and if it's already an integer then populate the 'length_val' field with the number.  If it's a string then strip it down to the number (typically the greater than, less than or tilde off the front) and populate the 'length_val' field.

In other words:
If 'total_length' = STRING, then send only the number to 'length_val'
Otherwise, if 'total_length' = INTEGER, then calculate the field.

I've been using update cursors to do a lot of the rest of the work, but that doesn't have to be the case.

Row[3] = 'total_length'
Row[17] = 'length_val'

#---First Attempt---

with arcpy.da.UpdateCursor(aoi_table_gdb,aoi_fields) as cursor:
    totlen = "total_length"  #<--- This was just an idea
    fields = arcpy.ListFields(aoi_table_gdb)
    for field in fields:
        if field.name = '"'+ totlen +'"' and field.type = String:
            for row in cursor:
                if row[3] == None:
                    row[17] = None
                elif row[3].find(str('~ ')) > -1:
                    row[17] = row[3].lstrip(str('~ '))
                elif row[3].find(str('> ')) > -1:
                    row[17] = row[3].lstrip(str('> '))
                elif row[3].find(str('< ')) > -1:
                    row[17] = row[3].lstrip(str('< '))
                else:
                    row[17] = row[3]
        else:
            for row in cursor:
                row[17] = row[3]
        cursor.updateRow(row)

#---Second Attempt---

with arcpy.da.UpdateCursor(aoi_table_gdb,aoi_fields) as cursor:
    for row in cursor:
        #Length Value
        if row[3] == None:
            row[17] = None
            arcpy.AddWarning("A. Row: " + str(row[1]) + " - " + str(row[3]))
        elif row[3].isdigit():
            row[17] = row[3]
            arcpy.AddWarning("B. Row: " + str(row[1]) + " - " + str(row[3]))
        else:
            length = filter(unicode.isdigit, row[3])
            row[17] = length
            arcpy.AddWarning("C. Row: " + str(row[1]) + " - " + str(row[3]))


texts: Error when writing SHAPE to WKT with pyqgis and QGIS 2.14

texts:
It seems the answer linked here doesn't work in latest version of QGIS 2.14

layer = QgsVectorLayer("/home/reyman/Projets/cours/coursPython/pyQgis/shape/passages/passages.shp", "passages", "ogr")
if not layer.isValid():
  print "Layer failed to load!"

QgsVectorFileWriter.writeAsVectorFormat(layer, "wkt.csv", "utf-8", None, "CSV", layerOptions='GEOMETRY=AS_WKT')

return the warning :

Warning 6: layer creation option 'G' is not formatted with the key=value format
Warning 6: layer creation option 'E' is not formatted with the key=value format
Warning 6: layer creation option 'O' is not formatted with the key=value format
Warning 6: layer creation option 'M' is not formatted with the key=value format
Warning 6: layer creation option 'E' is not formatted with the key=value format
Warning 6: layer creation option 'T' is not formatted with the key=value format
Warning 6: layer creation option 'R' is not formatted with the key=value format
Warning 6: layer creation option 'Y' is not formatted with the key=value format
Warning 6: '' is an unexpected value for  layer creation option of type string-select.
Warning 6: layer creation option 'A' is not formatted with the key=value format
Warning 6: layer creation option 'S' is not formatted with the key=value format
Warning 6: layer creation option '_' is not formatted with the key=value format
Warning 6: layer creation option 'W' is not formatted with the key=value format
Warning 6: layer creation option 'K' is not formatted with the key=value format
Warning 6: layer creation option 'T' is not formatted with the key=value format


texts: osm2pgsql: modifing the default.style - newly added columns are empty

texts:
I would like to import a osm-db (from planet) with the "name_de" and "source" column. So I added/changed those lines in the default.style:

node,way   source       text         linear
node,way   name_de      text         linear

After import, both columns are created but empty. What I am doing wrong?


texts:
I am creating a PyQT and QGIS application and im trying to add a combobox to the UI. I have defined the layout of my UI in a separate module to the main module. I know it is an issue with the combobox as this error only occurs when I add it. The application runs but when I close it this error is displayed in the terminal - QPixmap: Must construct a QApplication before a QPaintDevice

Here is the code for the UI

import sys
import os

from qgis.gui import *
from qgis.core import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *

from ui_mainWindow import Ui_MainWindow

class DeforestationMapWindow(QMainWindow, Ui_MainWindow):
    def __init__(self):
        QMainWindow.__init__(self)
    self.setupUi(self)

    self.connect(self.actionZoomIn, SIGNAL("triggered()"), self.zoomIn)
    self.connect(self.actionZoomOut, SIGNAL("triggered()"), self.zoomOut)
    self.connect(self.actionPan, SIGNAL("triggered()"), self.panMode)
    self.connect(self.actionAddArea, SIGNAL("triggered()"), self.addArea)
    self.connect(self.actionRemoveArea, SIGNAL("triggered()"), self.removeArea)
    self.connect(self.actionGetInfo, SIGNAL("triggered()"), self.getInfo)
    self.cb.currentIndexChanged.connect(self.selectionchange)

#handlers for the events
def zoomIn(self):
    pass
def zoomOut(self):
    pass
def panMode(self):
    pass
def addArea(self):
    pass
def removeArea(self):
    pass
def getInfo(self):
    pass
def selectionchange(self,i):
    print "Items in the list are :"

    for count in range(self.cb.count()):
        print self.cb.itemText(count)
    print "Current index",i,"selection changed ",self.cb.currentText()

def main():

    app = QApplication(sys.argv)

    QgsApplication.setPrefixPath(os.environ['QGIS_PREFIX'], True)

    QgsApplication.initQgis()

    deforestMapWindow = DeforestationMapWindow()
    deforestMapWindow.show()
    deforestMapWindow.raise_()

    app.exec_()
    app.deleteLater()
    deforestMapWindow.close()

    QgsApplication.exitQgis()

if __name__ == "__main__":
    main()

Here is the code for the UI

from PyQt4.QtCore import *
from PyQt4.QtGui import *
import resources
class Ui_MainWindow(object):
    def setupUi(self, window):
        window.setWindowTitle("Deforestation")

        self.centralWidget = QWidget(window)
        self.centralWidget.setMinimumSize(900, 400)
        window.setCentralWidget(self.centralWidget)

        self.toolbar = QToolBar()

        window.addToolBar(Qt.TopToolBarArea, self.toolbar)

        buttonIcon = QIcon(":/resources/zoomIn.png")
        self.actionQuit = QAction(buttonIcon, "Quit", window)

        buttonIcon = QIcon(":/resources/zoomIn.png")
        self.actionZoomIn = QAction(buttonIcon, "Zoom In", window)

        buttonIcon = QIcon(":/resources/zoomOut.png")
        self.actionZoomOut = QAction(buttonIcon, "Zoom Out", window)

        buttonIcon = QIcon(":/resources/pan.png")
        self.actionPan = QAction(buttonIcon, "Pan", window)
        self.actionPan.setCheckable(True)

        buttonIcon = QIcon(":/resources/getInfo.png")
        self.actionGetInfo = QAction(buttonIcon, "Get Info", window)
        self.actionGetInfo.setCheckable(True)

        buttonIcon = QIcon(":/resources/draw.png")
        self.actionAddArea = QAction(buttonIcon, "Draw Area", window)

        buttonIcon = QIcon(":/resources/remove.png")
        self.actionRemoveArea = QAction(buttonIcon, "Remove Area", window)

        self.cb = QComboBox()
        self.cb.insertItems(1,["One","Two","Three"])

        self.toolbar.addAction(self.actionZoomIn)
        self.toolbar.addAction(self.actionZoomOut)
        self.toolbar.addAction(self.actionPan)
        self.toolbar.addSeparator()
        self.toolbar.addAction(self.actionGetInfo)
        self.toolbar.addAction(self.actionAddArea)
        self.toolbar.addAction(self.actionRemoveArea)
        self.toolbar.addSeparator()
        self.toolbar.addWidget(self.cb)

        window.resize(window.sizeHint())

I am running this on OS X
I would love for someone to explain what I am doing wrong and help me fix this.


texts: Selecting lowest n values in from attribute table in ModelBuilder?

texts:
I have a shapefile with a "Random Number" Column. I'd like to use ModelBuilder to Select and export the lowest N values. I've seen this post about selecting the top 10 records which is pretty similar.

Selecting limited TOP n from specific column in ArcGIS query builder?

But the problem I'm having is this post uses a Query Builder, but I'm not sure how to adapt this for model builder into the 'Select tool'. I'm using ArcGIS desktop


texts: Download Google Earth Engine Collection to Local Drive

texts:
I am fairly new to both python and Google Earth Engine, so please bear with me if the following question is kind of stupid because it ignores some basic knowledge.

I would like to process some Sentinel-2 hyperspectral images available from Google Earth Engine, but, as I said I am not very familiar with python. My idea was, therefore, to use python to get the images to my drive (in some kind of standardized format such as netcdf, BSQ, etc.) and then process them locally with MATLAB for which I already have an analysis pipeline.

Nevertheless I can't seem to find any way to do that: all the info I find clearly state how to export the image to Google Drive (I can do it and it works), but there's no clear info on how to get it on the local hard drive. Is that possible?

Thanks for the answers!


texts:
I've added a point on map like this but it is static in size.

How to resize that point on zooming?

Map on zoom out


texts: Google Basemap doesn't align with data after export

texts:
Help!

I put together a map in QGIS and everything looks good in the composer view, however when I export to PDF, the projection/data seems to have shifted from the Google Open Layers basemap.

I followed the instructions for adding the Google basemap from here:
https://maps.cga.harvard.edu/qgis/wkshop/basemap.php

Part of this guide had me change the project CRS to WGS 84 Pseudo Mercator.

Since I was having trouble, I attempted to fix the issue by right-clicking on each layer and changing the CRS to WGS 84 Pseudo Mercator as well. This made it worse.

The attached image shows what my PDF looks like at this point.


texts: Database Server Node vs. Database Connection Node

texts:
What are the main differences between the Server Node and the Connection Node?

The only difference I have been able to find is that the Server Node is intended for administrative tasks while the Connection Node is used when accessing data on one geodatabase.


texts: How To Add Values to an Attribute Field Based on a Selection - Python

texts:
Alright another problem that I have.
I am trying to give a priority ranking to lots based on various criteria.  Basically I have 7 areas that are criteria, if lots overlap each other the criteria is added up and given a sum value (so the highest number could be 7).  Now I am trying to create a python script that basically goes through my criteria and adds a value to a field added to an attribute table

For example, I would do a select by location for lots that are within 30m of wetlands - if YES then the field "CRITERIA" would be given a 1, if no it would be a 0
this same selection idea would be done for the rest of my criteria - then what I need is to have it total up the 1's that are in each row into a separate "RANKING" field.


texts: Using input layer as output in PyQGIS for QGIS Modeler?

texts:
I have written a simple Processing script to rename a given field. It works successfully, changing the input layer's field name. However, I have not been able to produce an Output. This is a problem because I would like to use this script as an interim step in a Graphical Modeler tool, which requires an Output to pass on to the next step.

def processAlgorithm(self, parameters, context, feedback):

    source = self.parameterAsVectorLayer(
        parameters,
        self.INPUT,
        context
    )
    input_field = self.parameterAsString(
            parameters,
            self.INPUT_FIELD,
            context)
    new_name = self.parameterAsString(
            parameters,
            self.NEW_FIELD_NAME,
            context)

    if source is None:
        raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT))

    # Rename the field
    source.startEditing()
    source.renameAttribute(source.fields().indexFromName(input_field), new_name)
    source.commitChanges()

    return {'MODIFIED_LAYER': source}

The last line does appear to return a Layer reference:

Execution completed in 0.07 seconds
Results:
{'MODIFIED_LAYER': <qgis._core.QgsVectorLayer object at 0x0000022D82850E58>}

However, when I use the script in a Model, it allows me to pipe the output into another algorithm, but at runtime, it does not work.

Input Parameters: # for 'Basic statistics for fields'
{ FIELD_NAME: 'placetype', INPUT_LAYER: None }

I am not sure what approach I should use to produce output which can then be used in the Modeler


texts:
Processing modeler: GDAL "clip raster by mask layer" tool

I have a raster layer that I wish to be clipped by a polygon layer (buffer zone). I intend to repeat to process many times in a more elaborate sequence for my project thus I aim to build an automated process in QGIS's processing modeler. When running the process, I get the following error message:

"Error executing algorithm Clip raster by mask layer 'NoneType' object has no attribute 'GetGeoTransform'"

NOTE: this is only the clipping sequence. You can see the log here:

Any ideas on how to solve this? The GDAL tool "raster clipper" works fine when I use it without the processing modeler.


texts: Summarise duration of points within polygon

texts:
I have table with points that are gps coordinates for multiple vehicles.

I have another table of multiple polygons that the vehicles travel through.

I’ve managed to write a python script that filters the vehicles by returning only the points inside the polygons. Now I need to figure out the duration of each vehicles inside of each polygon.

Unfortunately, each vehicle enters and exits each polygon multiple times so I can’t simply work out the difference between the max and min timestamps as that would also include time spent outside the polygons.

Is anyone able to point me in the right direction for how I might achieve this, either using SQL on its own or in combination with Python?

One way I’ve thought of is, instead of only returning the points in polygons is to add a Boolean column based on whether the row is an entry inside/outside. I’m not sure how to write this SQL query though.

I’d then sort the table on timestamp and iterate through it to look for when that value changes.

I’ve also thought of converting to lines and finding the intersects but then I’d lose the timestamp value, which is essential to the whole operation.

Any ideas gratefully received.


texts: Using python to access the attributes of features within the current extent?

texts: The script I'm writing requires a specific attribute of the building polygon that the data frame extent is currently focused on. Is there a way to grab this data?


texts: I want to add a circle which stays the same size

texts:
As stated in my title, I want to draw a circle with a fixed diameter. It should not change when I am zooming in or out. How do I do that? (I am using the CartoDB Service).

I am trying to visualize how many points on the map are within a 10 mile range for example.


texts:
So, basically I've created a QgsRubberBand which copies feature from layer. And I need to be able visually move this rubberband on mouse event 'move'. Actually same visualisation like during  self.iface.actionMoveFeature()

class mMapTool(QgsMapTool):

def __init__(self, mainLayer, canvas):

    self.mlayer = mainLayer
    self.canvas = canvas
    self.active = False
    QgsMapTool.__init__(self, self.canvas)
    self.rubberBand = QgsRubberBand(self.canvas, QGis.Polygon)
    mFillColor = QColor(0, 0, 0, 88);
    self.rubberBand.setColor(mFillColor)
    self.rubberBand.setWidth(2)
    self.reset()

def reset(self):

    self.startPoint = self.endPoint = None
    self.isEmittingPoint = False
    self.rubberBand.reset(QGis.Polygon)

def canvasPressEvent(self, e):

    global mdict
    mdict = {}

    if e.button() == 1:
        self.isEmittingPoint = True
        for feature in self.mlayer.selectedFeatures():
            centroid_tuple = (feature.geometry().centroid().asPoint().x(),
                              feature.geometry().centroid().asPoint().y())
            mdict[feature.id()] = [centroid_tuple]
            geo = feature.geometry()
            break

        self.rubberBand.addGeometry(geo, None)

def get_canvas_point(self, position):

    x = position.x()
    y = position.y()
    trans = self.canvas.getCoordinateTransform()
    return trans.toMapCoordinates(x, y)

def canvasMoveEvent(self, e):
    if not self.isEmittingPoint:
        return
    pos = e.pos()
    point = self.get_canvas_point(pos)
    self.rubberBand.movePoint(point)
    self.rubberBand.show()

But my current code doesn't work. I tried different options, searched for solution everywhere but it seems like I'm only one with this problem.


texts:
What does "OFT" stand for?

I'm getting the message "Invalid OFT field name given" in a Python program that uses GDAL, and before I debug it, I'd like to know what it means.


texts: How to import OSM relations to PostgreSQL databases?

texts:
I have created PostgreSQL databases using Osm2pgsql. I have all the lines and points data.

How can I import the relations stored in OSM file? for example turn restriction?
Is there an easy way to visualize it in QGIS? similar to JOSM symbols?


texts: How to connect / 'join' lines selectively but automatically in QGIS

texts:
I had a similar problem a while back with How to join (draw new) lines that are close but not touching - large dataset
This is the same dataset, however I found a problem with what I did before. The line sections where joined with no discretion and some that were joined still don't join when it comes to path finding algorithms.
The main dataset is a vector layer with lines (10k features). Some cross and intentionally 'join', others cross and intentionally don't join and some don't touch but should.
I have since been sent an additional layer with points of where lines should connect/intersect/join.
Part of my problem solving this has been the correct use of language as 'Join' and 'Intersect' is especially poor because it's typically used for 'joining' attributes and intersecting features.
In addition to this I'm not familiar enough with path finding and if there's a difference between line sections crossing and 'actually' intersecting... for my purposes I need to ensure that the lines intersect so a path finding algorithm will work (Network analysis -> 'Shortest path' seems to do the job)

Point where lines should join/connect/intersect - 1
Approx 25mm gap

Point where lines should join/connect/intersect - 2
Approx 0.03mm overshoot

Both issues seem to block path finding algor.

Is there a simple way to do this?
Methods I'm thinking about are:
Stopping short:

Snapping but somehow only for a subsect of the layer (e.g. union with the buffered points?)
Points to line/path (but I assume this draws lines between points which won't solve my problem

Overshoot

I think this one should be simple with a split tool, but again I need help to figure out how to only do it in the proximity of a junction point.


texts: progress bar/loader/spinner for map.locate

texts:
Im using Leaflet.Control.GPS (https://github.com/stefanocudini/leaflet-gps) to locate the users location.

How to add a progress indicator, or a spinner (or some other UX hint) which should appear only when "map.locate" is trying to find the current position of the user?

I already added via Leaflet.Loading (https://github.com/ebrelsford/Leaflet.loading) a spinner to my map, but only works when loading tiles.

Maybe it is possible to combine them?


texts: Different-sized areas bias in thematic maps and possible solutions

texts:
In my organization, we are looking for a possible solution to deal with the different-sized bias in thematic maps at census tract level/block-group level. We need to present our data as raw counts/totals. Therefore, normalization does not make sense. Other solutions such as cartograms don't make sense either.
We have two problems with large census tracts/block groups. First, larger areas naturally draw the map users' attention, especially if they are styled with darker colors. Second, the population in large census tracts is not evenly distributed like in the small census tracts. In large census tracts, the population is spatially clustered in a few small areas within the tract.
Here is my solution so far. I took advantage of the Microsoft building footprint data with all building footprints in the United States as a proxy of where the population is within large census tracts. Then, I created a hexagonal tessellation of 310 meters of side and 0.25 squared Km of area. On the left map, you can see the Orlando (FL) area's result if I keep all hexagons with one or more building footprints inside. The right-hand side map contains the original census tracts.

In that area of Orlando, there are one million building footprints and 40,000 hexagons. The hexagons contain between 0 to 350 buildings. If I zoom in on the left and side map, you can see that the hexagons only cover areas with buildings inside.

My next step will be getting rid of non-residential buildings by using land-use cover data. Finally, I am thinking of using a density-based cluster algorithm to remove the noise: isolated buildings that are not relevant for our business context:

In the end, we are going to use hot spot analysis. Interestingly, the output of the hot spot analysis is quite different if we use hexagons instead of census tracts:

Do you have any idea to share with me to improve my approach or suggest a better approach?


texts: How to centralize and change the center of the leaflet map dynamically

texts:
I have a map that covers the entire screen and I have a layer that overlays the map. That layer goes from bottom to top.
the user interacts with this layer by clicking on different options and these options draw markers on the map.

When I centralize the map, the markers are hidden behind the layer that superimposes the map, thus eliminating the visibility of the marker.

How could I centralize the map and then move it one third above the layer?
This question is a remake of
How can I change the location center of a map using Leaflet API?

Try doing it this way

var offset = map.getSize().y * 0.15;
map.panBy(new L.Point(0, -offset), {animate: false});

But every time I click on an option in the layer to create a marker, the map moves and I lose the center.

jsfiddle

Excuse the example image
The idea is that after the map makes the setView place its center above the blue layer that superimposes the map so as not to hide the markers.


texts: Updating attribute table of base shapefile from .csv file?

texts:
I am new to python programming. I am trying to create a script that can update attributes in a base shapefile based on a .csv file when changes are made to the .csv file. I also want the script to create new point features from the .csv file if they do not already exist in the shapefile.

Does anyone have a similar script already?

Here is what I have so far:

import arcpy, os, sys
arcpy.env.overwriteOutput = True
outFC = r"c:\Users\abneyr\Documents\Spatial_Data\Test.shp"
inputFile = r"c:\Users\abneyr\Documents\Python_Class\Study_Excel.csv"
coordsys = r"c:\Users\abneyr\Documents\Python_Class\lab08data\lab08\ShapeFiles\WGS_1984.prj"

arcpy.env.workspace = os.path.dirname(outFC)
featClass = os.path.basename(outFC)
if arcpy.Exists(featClass):
    arcpy.Delete_management(featClass)

arcpy.CreateFeatureclass_management(arcpy.env.workspace, featClass, "MULTIPOINT", "", "", "", coordsys)
fileRead = open(inputFile, 'r')
headerLine = fileRead.readline()
firstValue = headerLine.split(",")
nameIndex = firstValue.index("Name")
lines = fileRead.readlines()
fileRead.close()

boxDictionary = {}
with arcpy.da.SearchCursor(featClass, "*") as seCursor:
    for line in lines:
        Value = line.split(",")
        key = Value[nameIndex]
        if key in seCursor:
            continue

        with arcpy.da.InsertCursor(featClass, "SHAPE@") as isCursor:
            if key not in boxDictionary:
                boxDictionary[key] = arcpy.Array()
                print("Add {0} into the record".format(key))
                if "Latitude" in line:
                    continue

                xCoord = float(Value[1])
                yCoord = float(Value[0])
                boxDictionary[key].add(arcpy.Point(xCoord, yCoord))
                print("Add {0} into the shapefile".format(key))
            isCursor.insertRow([arcpy.Multipoint(boxDictionary[key])])


texts: How to write a loop selection by location and export in QGIS?

texts:
I would like to be able to make selections in sequence between buffers and surrounding areas. This would lead to areas (polygons) being selected repeatedly by 2 or more buffers and I want to be able to make selections and exports, listing every export to belong to a number (or numbers) of polygons

For instance:

I start by selecting the first line of shapefile A and I want to select all the polygons of shapefile B that are inside of shapefile A (B polygons contained in A), to then export them as shapefile_1.

And this process to be repeated for x records in A to produce more shapefiles (shapefile_1, shapefile_2, ... shapefile_x).

How can this be done in QGIS?


texts:
I'm configuring a new server with GeoServer and when I try to do a GetFeatureInfo query gets no feature found. I already check extension, bounds, names... and everything seems to be ok.

Any idea?

http://localhost:8080/geoserver/wms?bbox=-4.01855,36.75382,-4.01833,36.75390&styles=mapama&format=jpeg&info_format=text/plain&request=GetFeatureInfo&layers=sigpac2019:pg_malaga_2019_muni_94&query_layers=sigpac2019%3Apg_malaga_2019_muni_94&width=550&height=250&x=170&y=160

"no features were found"

The query is correct, the bbox are into layer limits.


texts: Extracting channel banks from DEM and channel network

texts: Currently I am working in project that require from me extracting channel banks. I have been thinking about the best method that I can use to come up with good function that extract the location of channel banks. I don't want to use borrowed software such as Geo-ras or HEC-Ras. My main focus is to come up with a function which I can apply to many DEMs based on the the shape of the cross-section and extract the width and height from it. It is similar to the this picture. At the current stage, my thought is create multiple points at one point from the channel network and then create multiple consecutive points across the channel with the same interval and use their elevation to create slopes and their is a change in the slope (threshold that will create my first bank point until the it finds the second bank then it stops. I believe that what they call slope breaking method. From that I can extract the depth and width in the channel correctly.


texts:
I have many nc files and for each one I would like to create a tif file with the average of the bands using r.series of grass processing in QGIS.
I wrote the code below but it doesn't work but no tif file is created and it gives me no error messages.
Using the processing tool, with one file it works.
import glob, os
path="C:/Users/Daniele/Documents/"
os.chdir(path)
for file in glob.glob("*.nc"):
        input_layer = (path+file)
        output_layer = input_layer.replace(".nc",".tif")

    param = { '-n' : False,
                'GRASS_RASTER_FORMAT_META' : '',
                'GRASS_RASTER_FORMAT_OPT' : '',
                'GRASS_REGION_CELLSIZE_PARAMETER' : 0,
            'GRASS_REGION_PARAMETER' : None,
            'input' : input_layer,
            'method' : [0],
            'output' : output_layer,
            'quantile' : '',
            #'range' : [nan,nan],
            'weights' : ''
            }

    processing.runAndLoadResults('grass7:r.series', param)
    ###########r.series average of bands of file nc

    ####### Reload your output layer ########

    #output_layer = QgsRasterLayer(path, file_tif,'ogr')


texts: DEM raster data to another raster with different resolution and extent

texts:
I have two raster datasets:

One fine resolution DGM and another one, with different resolution, in which I want to have the mean elevation from the DEM (if possible the mean of the cells from the DGM which intersects it)

The raster have different resolution and different extend.

How can I get the data from my DEM to my other raster


texts: True Curves change to densified lines upon Topology validation. Fixing 500,000+ slivers in parcel layer

texts:
We are experiencing an issue as detailed by this technical article: Bug: Curve features are converted to densified line segments from an arc in a SQL Server geodatabase after compression.

We are using ArcMap 10.7.1 Advanced. We have talked to ESRI technical Support and, they feel, the problem should not be in our database now. We recently went through two upgrades, originally editing in 10.5.1 then upgraded quickly to 10.6 and finally upgrading to 10.7.1. The Problem could be due to the ArcMap upgrades which do not have patches for the problem that did have a software patch in 10.5

We have 500,000+ topology errors in our parcel layer, due to slivers and overlaps of the true curves being generalized to densified lines. Does anyone have a recommendation of automating all or part of the process of cleaning this up short of manually aligning to edge on every sliver?

See the example screenshot below, where the blue line represents the parcel line after being densified, and the red line is the line the polygons were built from. I need to realign the polygons(blue line) to the polygon line boundaries (red lines):

Unfortunately, rebuilding the polygons from the lines is not an option as our cartography staff has made thousands of edits to attributes and geometry of the polygons in the last 6 months since this issue appeared. We would have to go back 7 months and attempt to track what polygons have been split, merged, reconfigured, resized, annexed, as well as trying to figure out what attributes changed and redo all of that.  We are not the primary editors to the tax parcel data.

We seek any ideas, resources, personal experience attempting a task like this, as this is a major issue for us. Fortunately, we believe we have stopped the process of creating these errors to begin with so no more slivers are being created, I just need to fix the 500,000+ that have already been created.


texts: Output shows symbols too thick. (QGIS 3.8.2)

texts:
I've loaded data from a CSV-file without problem and assigned rule-based symbols to it. These are SVG of my own design with a line thickness set at 0,2 mm (yes, I know that's thin).
Everything looks fine until I export the map. What should be 0,2 mm lines become more like 2 mm, making the symbols quite useless. I've had the problem before and managed to work my way around it, but I'm afraid I don't remember how. I feel like I'm overlooking something rather obvious but have no idea what…. Any ideas?

Right layout

Wrong output


texts:
I am trying to create a GIS app using Java, Spring and Hibernate.

I send a JSON from Postman like this:

{
    "pointName":"any location",
    "geom":"POINT(-71.060316 48.432044), 3857"
}

The JSON is received by the Controller Class and sent to a Converter Class where the information that arrives as String is converted to geometry with the following method:

public static Geometry wktToGeometry(String wellKnownText)
              throws ParseException {

                return new WKTReader().read(wellKnownText);
            }

But once the information is processed by this method, 3857 dissapears and I have the following error when I try to save in the DB:

org.postgresql.util.PSQLException: ERROR: Geometry SRID (0) does not match column SRID (3857)

This is because the reference system has been removed in the wktToGeometry method, and it is interpretated as reference system 0.

How should I format the JSON?


texts: How to select attribute over another through model builder?

texts:
I have a gas main feature that I was asked to add MEASUREDLENGTH field to. I currently have a model that uses Shape_Length to sum overall pipe, but was asked to alter the model to use the new MEASUREDLENGTH field if populated, else use the Shape_Length field to sum overall pipe.

Can this be done with the Calculate Tool?


texts: What data structure to use for raster data when using rasterio?

texts:
I want to write a method in Python that reads a raster (GTiff) using rasterio and returns its contents. Then, some processing is done. Finally, I want to write the data to a new GeoTiff. As far as I understand, with rasterio, I use the read() method to store the bands in a pandas ND-array and the source.profile to store it's metadata. This seems rather cumbersome to me compared to vector data.

When dealing with vector data, I can use a geopandas geodataframe to store the entire contents of a shapefile or a geopackage like this:

everything_from_a_shapefile = geopandas.read_file(input_shapefile)

Is there some similar data structure for storing raster data conveniently?

everything_from_a_gtiff = ???

I am new to this so maybe I am missing something obvious.


texts:
For ArcMap 10.3x, is there is a Python or ArcPy code or model or toolbox that will allow me to choose the layers that I want to include in a map and export as PDF maps?

i.e. let's say I want to make 10 county-level maps of schools, stores, hospitals, etc. and always include base layers of place names, roads, street names, etc. So in this case, I'll specify to make 10 PDF maps for county 1, county 2, county 3,..county 10 but always include base layers of place names, roads, street names, etc.?


texts:
I'm trying to install and run the Boundless SDK according to:
http://suite.opengeo.org/opengeo-docs/webapps/install.html. I've tried running:

  $ suite-sdk
  bash: /usr/bin/suite-sdk: Permission denied
  $ sudo suite-sdk // with and without the help and create commands
  sudo: suite-sdk: command not found

So I know that I am calling suite-sdk from the right place. I have jdk6 and Ant both installed.

  java version "1.6.0_33"
  OpenJDK Runtime Environment (IcedTea6 1.13.5) (6b33-1.13.5-1ubuntu0.12.04)
  OpenJDK 64-Bit Server VM (build 23.25-b01, mixed mode)
  Apache Ant(TM) version 1.8.2 compiled on December 3 2011

My /etc/environment files is:

  PATH="/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games"
  JAVA_HOME="/usr/lib/jvm/java-6-openjdk-amd64/jre/bin"

(I've also tried the JAVA_HOME path with and without /bin and /jre)

Does anybody know what I am doing wrong?


texts: Is it possible to dock the attribute table using pyqgis? For instance, docking it to the bottom of the screen?


texts: Error reading data from geomesa(bigtable)

texts:
I'm having a problem to read data from Geomesa(BigTable).
I can connect to the source of data, but when I try to get the data returns the following error:
java.lang.UnsupportedOperationException: Pleaes use addRange(byte[], byte[]) instead. Pleaes use addRange(byte[], byte[]) instead.
What could be the cause? Anyone already have this error before?
I tried to change the geoserver dependencies for the ones on the geomesa documentation and tried to use the updated version of almost all jar files.


texts:
I'm using the leaflet plugin photon to request adresses like here :
https://adresse.data.gouv.fr/map/

It works fine and the geocoder was added to my map.
Now I would like to put the geocoder outside of my map in a searchbox.
Is there a way to do it ?

API_URL = '//api-adresse.data.gouv.fr';

var searchPoints = L.geoJson(null, {
    onEachFeature: function (feature, layer) {
        layer.bindPopup(feature.properties.name);
    }
});

function showSearchPoints (geojson) {
searchPoints.clearLayers();
searchPoints.addData(geojson);
};

var photonControlOptions = {
resultsHandler: showSearchPoints,
placeholder:' Ex : Rue Diderot Paris',
position: 'topright',
url: API_URL + '/search/?',
noResultLabel: 'Aucun résultat',
};

var map = new L.Map("map", {
center: new L.LatLng(49.047, 1.605),
layers: [layerstreets],
zoom: 13,
zoomControl: false,
photonControl: true,
photonControlOptions: photonControlOptions
});

searchPoints.addTo(map);


texts: Placing an image on a lat/lon point in QGIS

texts:
Can we place an image as a marker on a lat/lon point in QGIS using PyQGIS? Currently, I use QgsVertexMarker()... but I'd like an image on a point rather than a marker.

registry = QgsMapLayerRegistry.instance()

canvas = qgis.utils.iface.mapCanvas()

m = QgsVertexMarker(canvas)

 m.setCenter(QgsPoint(float(lon),float(lat))

 m.setColor(QColor(255, 255, 0))

 m.setIconSize(10)

 m.setPenWidth(10)

m.setIconType(QgsVertexMarker.ICON_X)


texts: Separating date and time from the same column into two columns

texts:
If I have date and time in the same column how can I separate them in two columns in ArcMap?
And how can I put the date in this format (year,month,day)?


texts: Would a nice GPU on the Geoserver server (let's pretend it's old school housing) help in any way, or make any difference at all with the rendering tasks ?


texts: Georefencing angled geological cross-section maps in QGIS?

texts:
Is it possible to georeference this type of angled cross section, preferably in QGIS (Section-II in the image) ?

I am very much new to this type of georeferencing.

Googled & couldn't find a proper lead


texts: What is the visualization parameter of Water Occurance Layer transiting from 1% opacity red to 100% opacity blue?

texts:
I want to display my water occurance layer whose value ranges from 0% to 100% as the sample in GEE, but I did not find the corresponding visualization parameters in its Data User Guide.

Dose anyone have ideas? Thanks in advance.


texts: Excluding multi geometry types into a single Polygon type

texts:
After i did an intersection with two shape using geopandas through the following topic (Intersecting two shape problem using geopandas). But it output a multi geometry type shape,which in futher step will impossible export to polygon.
Results show as below:

for i in data:
   geom = i['geometry']
   print geom

Output:

POLYGON ((139.9125 40.592043, 139.925 40.592043, 139.925 40.58371, 139.9125 40.58371, 139.9125 40.592043))
POINT (139.9125 40.58371)
LINESTRING (139.9125 40.58371, 139.925 40.58371)

How can i exclude Point and LineString from multi-types into a single geometry type as polygon?


texts: Any way to deal with multiple entries for one field in QGIS?

texts:
I'm mapping out parcel numbers from an excel spreadsheet, matching them to the county APN field in a shapefile , which is in the format ###-###-##. My spreadsheet has four fields that I've concatenated to create the whole APN. I've run into a problem:

My final field sometimes contains multiple numbers (when a particular map covers more than one parcel). Is there a way to select all of the numbers in this field and have them selected out i.e.  the concat field would appear as such "100-184-12,13,17" and I want all three parcels picked "100-184-12",""100-184-13", and "100-184-17".


texts: Error 999999 calculating Euclidean Distance?

texts:
I am trying to calculate Euclidean Distance in ArcMap and I keep getting:

  ERROR 999999: Error executing function.

  Failed to execute (EucDistance).

I gave a maximum distance of 5000m, output cell size of 10m, set the extent and snap raster env var to my DEM. I tried giving different names to my output distance raster, disabling the background geoprocessing, setting different env var values, but I still getting the same error.


texts: Displaying Google Earth Historical Imagery in QGIS

texts:
I am currently carrying out research aimed at showing the evolution of plant cover in urban areas. I opted for the calculation of the vegetation index from Sentinel2 type images that I would like to refine the results using the historical images from Google Earth.
With the new version of QGIS, is it now possible to display historical images from Google Earth or is there another method knowing that the area of the study area is about 500km²?


texts: Is there a way to display more than 10 categories in cartodb?

texts:
I'm creating a map in CartoDB. I have roughly 18 categories, but each of those has a sub-category, which is 'yes' or 'no'.

Currently I am displaying the 18 categories in two different layers to get around the 10 category max limit.

How can I display all of these categories in as few layers as possible?


texts:
I want to programmatically add a header and a footer as I use atlas with multipages for one feature but I can't predict the number of pages for the table depending on the feature.
In order to add the header and footer, I want to copy a QgsLayoutItemGroup and all the associated items to position them on other pages. I can easily through the GUI, copy and paste the QgsLayoutItemGroup but I don't see how I can do it with PyQGIS and make change automatically. Pointers welcome!

Edit

I've changed the original function code to only deal with group to remove references between objects and rename them as the answer provided was keeping some unwanted references between group. It's clearly not enough generic but it fits my use case: a group with items but no other group within it.

from qgis.PyQt.QtCore import QUuid
from qgis.PyQt.QtXml import QDomDocument
from qgis.core import QgsReadWriteContext

def duplicate_item_group(item, layout, position=None):
    doc = QDomDocument()
    element = doc.createElement('Copied items')
    context = QgsReadWriteContext()

    group_uuid = QUuid.createUuid().toString()
    item.writeXml(element, doc, context)
    if type(item) == QgsLayoutItemGroup:
        print(len(item.items()))
        for subitem in item.items():
            subitem.writeXml(element, doc, context)
    doc.appendChild(element)
    layout_items = doc.elementsByTagName('LayoutItem')
    uuids = []
    for i in range(layout_items.count()):
        item_node = layout_items.at(i + 1)
        if item_node.isElement():
            item_node.toElement().removeAttribute('uuid')
            uuid = QUuid.createUuid().toString()
            uuids.append(uuid)
            item_node.toElement().setAttribute('groupUuid', group_uuid)
            item_node.toElement().setAttribute('templateUuid', uuid)

    group_elements = doc.elementsByTagName('ComposerItemGroupElement')
    for i in range(group_elements.count()):
        item_node = group_elements.at(i)
        if item_node.isElement():
            item_node.toElement().setAttribute('uuid', uuids[i])
    # print(doc.toString()) # To debug the XML content
    if position is None:
        return layout.addItemsFromXml(element, doc, context)
    else:
        return layout.addItemsFromXml(element, doc, context, QPoint(*position))


texts: Qgis2threejs plugin results images with bleeding edges and wider background canvas

texts:
I have started using 3D Qgis2threejs plugin in an attempt to make some beautiful 3D images of LiDAR analysis I am currently working.

I have two questions regarding that:

I am getting that bleeding edge (see it on screenshot) - why and what I can do to remove it?
Why do I get a bigger background white surface than the DEM? Is there a way to fix it?

Screenshot:
 https://mega.nz/#!AJ8DEJ7J!1iCbR4MHYkI1-R3bYpx0Rb1aLATB-4Ee9KJVU-g2veQ

Settings:
https://mega.nz/#!cQNnjQwA!zVD0j1ZW-nM9h81kjLQqQAkVwXlXZpkZ1QN7o6_vWWA


texts: pgr_dijkstraCost() Many-to-Many exceeding RAM disponibility

texts:
I am working with PostGreSQL 9.6, pgRouting 2.4.1 and PostGis 2.3.1 on a Linux with 3Gb of RAM and 4Gb of SWAP.

I am trying to launch the following query :

SELECT start_vid as building, end_vid as plant, agg_cost as junction_length from pgr_dijkstraCost
(
'SELECT gid as id, i as source , j as target , length as cost from   network_edges'
, array(SELECT id from "network_nodes" where buildingID!='NULL')
, array(SELECT id from "network_nodes" where plant_name!='NULL')
, false
)

The network_nodes as 135 000 row and the network_edges as 270 000 rows. Every time I launch the query, it exceed the RAM in 30s and I get disconnected from the database. I found a some informations on the web : https://www.postgresql.org/docs/current/static/kernel-resources.html#LINUX-MEMORY-OVERCOMMIT. I tried to minimize the value of work_mem, shared_buffer and max_connection_number but it didn't change anything.

Is there any way that I can limit the memory usage of this query, even if it takes longer or does it need to load everything in memory ?

The only solution that I see is to reduce the number of building and then launch this request multiple times.


texts: Features total  Filtered total, geometry problems

texts:
I have a user who's been adding or editing polygons on a vector layer; he noticed that when he saved the layer that sometimes the polygon would lose its geometry and disappear from the canvas - but the attribute data remained in the table.  He's also received a few 'unable to package..' errors when saving the layer.

I've noticed that in the attribute table, the total number of features is different than the filtered number - but there are no filters applied (see attached screenshot)   I took a look at the records in the .dbf file and the number matches the filtered total, but also includes a completely blank row (no blanks in the attribute table)  We're completely stumped why adding these simple polygon shapes have caused so many issues.

QGIS version 3.6


texts:
My address is not on Google Maps.

If I search Google Maps for my address it does not recognize it.

How can I add my missing address to Google Maps?


texts: How to export all singleband raster pixel values into excel or .csv in ArcGIS Pro?

texts: I have NDVI rasters and I would like to export ALL pixel values into an Excel compatible format.  I would ideally like to have a column for pixel values and another column for frequency.


texts: Openlayers3 - reproject a point from EPSG:4326 to EPSG:2193

texts:
I need to reproject a point from EPSG:4326 to EPSG:2193.I tried to achieve that doing the following:

Updated: Here is the little Angular2 service that I wrote: It takes a coordinate and and tries to reproject it. Afterwards it returns a filter string.

import {Injectable} from '@angular/core';
import {OlService} from './ol.service';

@Injectable()
export class OlFilterService {

constructor(private olService: OlService) {

}

pip = (coords: [number,number]) => {
    let ol = this.olService.get();

    let reprojectedCoords = ol.proj.transform(coords, 'EPSG:4326', 'EPSG:2193');

    let filter = 'Intersects(the_geom,POINT(' + reprojectedCoords + '))';

    return filter;

}

}

This code silently fails on me. I just don't get anything or 'false' back. Is it simply not possible or am I missing something?


texts: Using ESRI character font ("Environmental & Icons") with formatting tag in layout text box

texts:
I am attempting to use a character from the ESRI Environmental & Icons font collection (A bird in particular). I went into my character map tool and found the one I want (symbol right parenthesis). With the formatting tags, the input looks like this:

<BOL>Country Y</BOL><FNT name="ESRI Environmental & Icons">  )</FNT>
Disease X
Total Avian Cases: 12,345
<FNT size="8">Since:  01 JAN 2016 </FNT>

I had no problems using a different character set (ESRI US Forestry 1, for example). I suspect it has something to do with the "&" in the font title.
I also tried using "ESRI Environmental & Icons Regular," which is how it's listed in the font folder on my system.


texts:
I tried it this way but I think using DWithin might be better for this situation. What would be the correct order to find duplicates which exist within 10km from each other?

SELECT n, ST_ClusterDBSCAN(geog::geometry, eps := .08, minpoints := 2) over () AS cid
from (
select *
from cities as t
inner join
    (select n dn
     from cities as t
     group by n
     having count(*) >= 2) dups
  on dups.dn = t.n
order by t.n
) d

EDIT: I guess I would have to do something like this:

select * from (SELECT *
    from (
    select *, ROW_NUMBER() OVER(order by t.n) as rownum
    from cities as t
    inner join
        (select n dn
         from cities as t
         group by n
         having count(*) = 2
        ) dups
      on dups.dn = t.n
    order by t.n
    ) d
    where mod(rownum,2) = 0) even,
    (SELECT *
    from (
    select *, ROW_NUMBER() OVER(order by t.n) as rownum
    from cities as t
    inner join
        (select n dn
         from cities as t
         group by n
         having count(*) = 2
        ) dups
      on dups.dn = t.n
    order by t.n
    ) d
    where mod(rownum,1) = 0
    ) odd
    where st_dwithin(even.geog,odd.geog,5000)

but this is confusing... maybe it's better to just to do DWithin first but I'm not sure how to do that.


texts: Splitting polygon by sections obtained from IDW interpolation

texts:
I used IDW to show the rainfall pattern over a region ... the attached picture show the result ... my area split into three ranges that show 3 different colors...

How can I split one original polygon to three based on the boundary of colors or can the boundary be shown by lines?

My goal is to apply the new boundary one different heatmap.


texts: R Raster package: write NetCDF with time dimension

texts:
I am writing a rasterbrick with multiple layers to a netcdf file using the raster package in R:

writeRaster(gridfile, "Data", "CDF", overwrite=TRUE,
     varname="Data",varunit="mm.d-1",longname="Monthly data",
     xname="lon",yname="lat", zname="time", zunit="month",bylayer=FALSE,NAflag=-9999)

Is there a way to have the data with time coordinates in the written ncdf file?
E.g. in the NetCDF Viewer Panoply I can plot the data in [lon][lat] but not in [time][lat].

Also how can specify the timesteps as specific months in a year instead of the timesteps starting from 1 and increasing?


texts: Using arcpy.mp to edit TEXT_ELEMENTS with no success

texts:
I'm attempting to take a user input in a script tool to modify the text within text elements on a layout. The following is the code I am using:

aprx = arcpy.mp.ArcGISProject(r"S:/Initiatives/Mailings/MailOuts.aprx")
lyt = aprx.listLayouts("Buffer Map")[0]

mapType = lyt.listElements("TEXT_ELEMENT", "Map Type")[0]
mapType.text = arcpy.GetParameterAsText(0)

location = lyt.listElements("TEXT_ELEMENT", "Location")[0]
location.text = arcpy.GetParameterAsText(1)

schedule = lyt.listElements("TEXT_ELEMENT", "Schedule")[0]
schedule.text = arcpy.GetParameterAsText(2)

I don't receive an error when I run the script tool but the text elements remain empty when given user input.


texts:
I would like to plot on QGIS these coordinates using a different color for the individuals and a different size of the symbol when the same individual is detected in the same zone (same coordinates).
I can give different colors to the individuals using the command: double click on the layer-->style-->"categorized". But I can't find out how to give different size to the symbols.

For example, the individual Gen26 is detected in the same position 4 times (Ua8546, Ua8547, Ua8552, Ua8559), and in another position twice (Ua8542, Ua8543). I would like to give bigger dimensions to the symbol when is detected 4 times than when is detected twice; of course, maintaining the same color for each individual.

@she_weeds and @kazuhito, I've found a problem: the method seemed to work but saw one thing.
The points are different in colors depending on individuals, also the points seem to have different sizes depending on frequency of occurrence BUT one individual detected only once in a point has the same size (big) of another individual detected different times in the same point.
I've simplified (and modified) the dataset I posted one week ago

Gen20 has the same dimensions of Gen16, but Gen16 is detected 5 times in that position, whereas Gen20 once. I used this function: count_distinct("ID", group_by:=geom_to_wkt($geometry)).
I'm sorry, I thought it worked but I can't understand how to solve this problem (I tried in these days to solve, but nothing).
Hoping in your another magic answer


texts:
I have a Leaflet map with a control to switch between various basemaps or toggle on/off some overlay data.

  //Basemaps are all a-okay
  //Start with no polylines in my overlayData
  items = [];
  flightpaths = L.layerGroup(items);

  overlayData = {
  "Videos": flightpaths
  };

  L.control.layers(baseMaps, overlayData, {position:"bottomright"}).addTo(map);

But my overlayMaps objects are coming in at random times. I'm using Leaflet.Draw to allow people to create polylines, so each time there is a draw:created a new object gets created that I want the controller to be aware of. There is some other analysis I do in here but the basic flow is-

map.on('draw:created', function (e) {
  item = e.layer;
  items.push(item);
  flightpaths = L.layerGroup(items);

  overlayData = {
      "Videos": flightpaths
  };

  L.control.layers(baseMaps, overlayData, {position:"bottomright"}).addTo(map);
});

But this (unsurprisingly) add just a second control to my map. How do I get the original control to recognize that overlayData has been changed?


texts:
How to filter data on attributes with a slider?

By example, i would like to filter weights (CHARGE in my layer) attribute for polygons without to change the filter formula.

The feature is available for time window. I don't know other solutions for non date/time type attributes.

How can i do?


texts:
How can I use the name of the data-column in a print-collection? (e.g. to use as headline)

I mean the column I choose at "layer properties" => Style => Graduated => column. I want to use it in legend or as headline.


texts: Which projection: Map from South Tyrol (Italy)

texts:
I need to find out what projection this base-map has for digitizing the information that was drawn on it. It's from South Tyrol (South Tirol) in Italy. I have the following information:

Title: Provincia di Bolzano / Provinz Bozen
Year: It was published 1980, so it could be anything before that

It has this overview:

Next to it is has this code: 4 III NO-B

And the page number: Foglio 74 / Blatt 74

The x-coordinates are given like this: E=1650000

And the y-coordinates: N=5183000

At the bottom of the map it says (some sort of copyright notice):

  Questo supporto consente solo copie a luce fredda

  Cartogr. Autorizz. A norma di legge

And then it says:

  E.I.R.A.

  Ente Italiano rilievi aerofotogrammetrici

  Rilievo aerofotogrammetrico eseguito con strumenti "Galileo-Santoni"

Edit:

I did some research and it might be a Gauß-Boaga (Gauss-Boaga) or Roma 40 projection, but I'm not sure yet.
Link: https://it.wikipedia.org/wiki/Proiezione_di_Gauss-Boaga


texts:
I'm new to geopandas and want to analyse a Geodatabase.

The information I received contains a zip archive of 45 files, containing

17 ATX files - each 5KB
2 files without extension - each 1KB
7 GDBINDEXES files - each 1KB
8 GDBTABLE files - each 1-18KB, one being 54.5MB
8 GDBTABLX files - each 6KB, one being 891KB
1 LOCK file - 0KB
2 SPX files - 5KB and 2.04MB respectively.

could someone tell me which ones I need to use - and how I import them?


texts:
I'm trying to make a slope raster from a lot of DEM`s.

The DEM rasters are in this file Geodatabase
C:\Bachelorarbeit\Daten\GIs\Originale\dtk25_dgm1.gdb
The slope raster should be saved in this file Geodatabase
C:\Users\Hans im Glück\Documents\ArcGIS\Default.gdb

Is it possible to use spaces in my pathname?

I have ArcGIS 10.3 and found this basic script on ESRI's website:

import arcpy
from arcpy import env
from arcpy.sa import *
env.workspace = "C:/sapyexamples/data"
outSlope = Slope("elevation", "DEGREE", 1)
outSlope.save("C:/sapyexamples/output/outslope01")

How to I execute this script on multiple  DEM's contained within a folder?


texts: How to put landmarks per floor of the building?

texts:
If I want to route from RoomA of floor level 5 to RoomB of floor level 1 inside the same building. How to add direction confirmation landmarks with tolerance buffer 1.5 that will not conflict the landmark points per level?

Below with the sample direction, has a "Pass by PA" which is not a landmark from level 5.

Sample direction:

1: Start at DIRECTOR

2: Go southeast on I3 L05 Corridor

2.1: Pass by EXEC DIRECTOR. It will be on the left
2.2: Pass by DIRECTOR. It will be on the right
2.3: Pass by DEPUTY DIRECTOR. It will be on the right

3: Turn left to stay on I3 L05 Corridor

3.1: Pass by **PA**. It will be on the right
3.2: Pass by SENIOR MANAGER. It will be on the right
3.3: Pass by SENIOR MANAGER. It will be on the left
3.4: Pass by SENIOR MANAGER. It will be on the right
3.5: Pass by DEPUTY DIRECTOR. It will be on the left

4: Turn right to stay on I3 L05 Corridor

...

end


texts: Convert xy coordinates into GPS coordinates

texts:
I have quite a large dataset with only a x and y coordinate which came from a program that is no longer usable. However I have the GPS coordinates of a few xy coordinates:

| place |   X  |  Y   |  latitude   |  longitude |
| H13   | 773  | 1214 | 53,36007006 | 6,46534807 |
| N13   | 1351 | 1236 | 53,36118361 | 6,46499224 |
| H10   | 760  | 942  | 53,36151658 | 6,46638025 |

So my question is how can I calculate all the other points that don't have a latitude and longitude (GPS coordinates)?

X represents the latitude and Y the longitude.

I quess the lat and lng of the other points will be less specific as there is no Z coordinate however all the points are quite close to each other so I guess the difference will be quite small.

The place column represent a location on a map based on a 10x10 meter grid square so (A-B = 10 meter A-C = 20 meter etc.) Where the characters are from A-O and are horizontal and the numbers from 1 - 19 vertical.


texts: Text that remains without any active layers

texts: I think this might be a simple question to answer, but I'm not sure how to do it. I turned off all the layers in my .mxd file, and I still have text that shows up on the screen. I don't know if I accidentally added this or something, but I don't know how to get rid of it.


texts: Is there any built-in way of limiting the request width/height of ImageWMS sources and have the map just scale the image up to fill the extent? I'm building an application that includes a bunch of WMS providers that don't support resolutions higher than 2048.


texts: Numerical Vertex Edit Error: latitude or longitude exceeded limits

texts:
I am trying to shift a point to a specific location.
I have followed the instructions here: How to manually edit point coordinates in QGIS?
but got this error code:

  An error has occurred while executing Python code:

  QgsCsException: forward transform of (9675.040722, 120072.456618)
  PROJ.4:  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
  +towgs84=0,0,0 +to  +proj=utm +zone=56 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs Error: latitude or longitude exceeded limits

I entered 0554339, 6879645.
OTF is turned on.
The layer CRS is GDA94/MGA 56


texts: How can i make it so the legend on the left bottom don't cover the map?


texts:
I have quite a problem, hope you'll be able to help.
I've got a very large shp layer that I'm loading via wms to my openlayers project. Due the fact that the layer is very big (2 millions polygons) It causing my geoserver to crash.
I have 2 solutions to this situation:

I wrote a js code to load the layer only in certain zoom levels. This solution works but i don't want it to be executed via the js code, but via the sld.
write a MaxScaleDenominator statement. when I try this, I know it executes (I see the polygons only at the selected views), but geoserver still crashes. As I can see it, It still loading all the polygons and just after it filtering it, causing it to crash.

Any idea how can I make the zoom filter without loading all the layer ?

Thanks!


texts: Does longitude/latitude info come directly from GNSS chip on Android devices?

texts:
I'm trying to directly communicate with the GNSS chip on my device and I am not clear whether the RAW GNSS information is coming directly from the GNSS chip or not. I have used the application provided in the link and it seems to read longitude/latitude information.

Is this information directly coming from the GNSS chip? Further, is the GNSS driver enough to extract this information? Or is longitude/latitude info calculated based on GNSS raw data? If it's the latter, how do I calculate longitude/latitude from raw GNSS (a reference would be sufficient)?

Note: I am not interested in very high accuracy location nor am I interested in indoors location retrieval (using Wi-Fi and Carrier services). I want to record the GPS (or Glonass or whatever) and save it as longitude/latitude format.


texts: How to calculate difference in altitude along lines using GRASS?

texts:
Using SRTM altitude information and OpenStreetmap roads, I want to calculate altitude differences (sum of meters uphill and downhill separated = two values) per road element.

What would be the fastest way to achieve this?


texts: How to get latitude /longitude in Degree Minute Second (DMS) using QGIS field calculator?

texts:
I  have come up with the following expression to convert 5.1234   to   5° 7`24.24"

toint ($x) || '° ' || toint((($x) - toint ($x)) * 60) ||'`' || substr( (tostring(((($x) - toint ($x)) * 60) - toint((($x) - toint ($x)) * 60)) * 60),1,5) || '"'

Here’s the problem :
for some  points such as 5.1234, it works. but for others,  it doesn’t.  I suspect the problem is the integer conversion which round up the decimal figures instead of truncating it.

Is there any other option? thanks.


texts: Draw a perpendicular line between point and line layer

texts:
I have to iterate over a point and line layer to draw the perpendicular line from each point to the nearest line. The QGIS geoalgorithm "DISTANCE TO NEAREST HUB" isn't helpful because the nearest hub isn't perpendicular.

I found this code, but how can I call the points in the Point Layer and the lines in the Line layer generally? How can I can iterate over them with a loop to find/create the perpendicular line in between?

# basic example with PyQGIS
# the end points of the line
line_start = QgsPoint(50,50)
line_end = QgsPoint(100,150)
# the line
line = QgsGeometry.fromPolyline([line_start,line_end])
# the point
point = QgsPoint(30,120)

def intersect_point_to_line(point, line_start, line_end):
     ''' Calc minimum distance from a point and a line segment and intersection'''
      # sqrDist of the line (PyQGIS function = magnitude (length) of a line **2)
      magnitude2 = line_start.sqrDist(line_end)
      # minimum distance
      u = ((point.x() - line_start.x()) * (line_end.x() - line_start.x()) + (point.y() - line_start.y()) * (line_end.y() - line_start.y()))/(magnitude2)
      # intersection point on the line
      ix = line_start.x() + u * (line_end.x() - line_start.x())
      iy = line_start.y() + u * (line_end.y() - line_start.y())
      return QgsPoint(ix,iy)

line = QgsGeometry.fromPolyline([point,intersect_point_to_line(point, line_start, line_end)])


texts: How to display international boundaries in a OpenLayers map?

texts:
I'm interested in displaying only the boundaries, coastlines and between countries, of the world through OSM. I'm a bit lost and I don't know what I should do.

I'm trying to display the OSM through OpenLayers on a website.

This is how I displayed the map through OpenLayers

 var view = new ol.View({
    center: istanbul,
    zoom: 6
  });

  var map = new ol.Map({
    target: 'map',
    layers: [
      new ol.layer.Tile({
        preload: 4,
        source: new ol.source.OSM() //OpenStreetMap here
      })
    ],
    // Improve user experience by loading tiles while animating. Will make
    // animations stutter on mobile or slow devices.
    loadTilesWhileAnimating: true,
    view: view
  });


texts: Erdas 2014 process time signature editor supervised classification

texts:
I am currently using erdas imagine 2014 with 64bit OS. Basic configurations of my computer are; i7 3820 3.60 GHz CPU, 16 GB memory (RAM), gforce GTX 680 SLI 256 bit Graphic card ( 2GPU ), 250 GB ssd. It takes 3.50 minutes to complete a process called signature selection from aoi with 1.2 gb layerstack images and less than 500 points.

I wonder if this speed is normal for such a application for a computer like this one, or this is a long time for the computer to finish.

It takes 7.30 minutes on a notebook with i3 processor and 6 gb ram with onboard graphic card.


texts: How to do operation after using over() for a shapefile in R?

texts:
I am copying the answer of Counts the number of points in a polygon in R

library("raster")
x <- getData('GADM', country='ITA', level=1)
class(x)
# [1] "SpatialPolygonsDataFrame"
# attr(,"package")
# [1] "sp"
set.seed(1)
# sample random points
p <- spsample(x, n=300, type="random")
p <- SpatialPointsDataFrame(p, data.frame(id=1:300))
proj4string(x)
# [1] " +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
proj4string(p)
# [1] " +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
plot(x)
plot(p, col="red" , add=TRUE)
res <- over(p, x)
table(res$NAME_1)

You see a nice overlapping of the + on the map (refer to that answer). I have similar data sets but in addition to id in the CSV file other two columns bh and dk.  What I need is just to compute the mean for all the values of id, bh, dk for each polygon. Any idea on this. In other words, in that answer he returned the number of points for each polygon using table(res$NAME_1), while here I need to calculate the mean for all values within a polygon.


texts: Shifting entire river channel in HEC-RAS?

texts:
This question is specific to the HEC-RAS software.

I'm trying to use the software to shift a river channel. I've used the Geometric Data Editor for the Cross Sections to edit the positions of the bank stations, and it seems to have worked:

But when I save and export this geometry file and map it, either using the RAS mapper or in GIS, there doesn't appear to be any change. Here's a part of the channel, with blue showing the original river channel, cross-sections & bank stations and red showing the edited ones:

I know that the geometry data is altered because if I view the data in a text editor, I see the new values there; but it doesn't seem to be present otherwise. There seems to be an HDF file associated with this layer, but I don't know what that is or how to edit it.

If I cannot fix it quickly I will just try a solution in Python.


texts: QGIS tools change special characters to unreadables

texts:
QGIS tools change special characters (õ, à) of attributes into unreadable characters.

Example:

# load 1 wfs feature
uri = "http://geo.vliz.be/geoserver/MarineRegions/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=MarineRegions:eez&Filter=%3Cogc:Filter%20xmlns:ogc=%22http://www.opengis.net/ogc%22%3E%3Cogc:PropertyIsLike%20matchCase=%22false%22%20wildCard=%22*%22%20singleChar=%22.%22%20escapeChar=%22!%22%3E%3Cogc:PropertyName%3Egeoname%3C/ogc:PropertyName%3E%3Cogc:Literal%3E*Gomera*%3C/ogc:Literal%3E%3C/ogc:PropertyIsLike%3E%3C/ogc:Filter%3E"
vlayer = QgsVectorLayer(uri, "test", "WFS")

# save as shapefile with UTF-8 encoding
QgsVectorFileWriter.writeAsVectorFormat(layer = vlayer, fileName = 'H:/Gomera.shp', fileEncoding = 'utf-8', driverName = 'ESRI Shapefile')

# run qgis tool
import processing
processing.run("native:fixgeometries", {'INPUT':'H:/Gomera.shp','OUTPUT':'H:/Gomera_fix.shp'})

In the output shapefile, the special characters have been changed into unreadables (compare left to right attribute tables):

The input shapefile and output shapefiles are both UTF-8 encoded.
This happens with several tools, for example with fixgeometries, buffer or grass7:v.buffer, both in the python console as well as in the GUI.

How can I solve this problem?
This seems something related to internal QGIS tools rather than problems with reading special characters as input as discussed here:

QGIS Processing algorithm cannot deal with special characters in input
ArcPy Field Calculation with special characters

I'm running QGIS 3.0 on Windows.


texts:
We have electrical conductivity maps of our fields created using a Veris sensor that are in a .dat file type. We want to import these maps into our Ag leader SMS database, but in order to do that, they need to be converted to shapefiles.

How would I go about converting them?


texts: How to Select an SRID? And how is it different from the SRID passed into ST_Transform?

texts:
Say I wanted to find geometries that intersect a circle at lat/long as outlined here How to Select all polygons or lines that lie at least partially within a circle with a radius?.

SELECT *
FROM buildings AS b
WHERE ST_DWithin(b.geom,
  ST_Transform(
    ST_SetSrid(
      ST_MakePoint({longitude},{latitude})
      ,27700)
    ,27700)
,{radius})

Now my data was the OS Open Local data: https://www.ordnancesurvey.co.uk/opendatadownload/products.html (top result).

I've imported the .shp files.

Running the .prj through here: http://prj2epsg.org/search produces:

    PROJCS["British_National_Grid",GEOGCS["GCS_OSGB_1936",DATUM["D_OSGB_1936",
SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],
PARAMETER["False_Easting",400000.0],PARAMETER["False_Northing",-100000.0],
PARAMETER["Central_Meridian",-2.0],PARAMETER["Scale_Factor",0.9996012717],
PARAMETER["Latitude_Of_Origin",49.0],UNIT["Meter",1.0],AUTHORITY["EPSG",27700]]

So in my postgis db I run:

SELECT srid FROM spatial_ref_sys WHERE srtext LIKE '%British_National_Grid%';

Which produces two results:
27700 and 7405

By inspection the former corresponds exactly to my .prj

The query above, however, finds no results for any buildings in London (51.5236, -0.1047) irrespectively of radius size (500, 5000, 500000).

Additionally, dropping the ST_SetSrid produces:

ERROR:  Input geometry has unknown (0) SRID

Why doesn't the ST_Transform set the SRID as documented here: http://postgis.org/docs/ST_Transform.html

Would I ever want the two SRID's to be distinct?

And why is the query not producing results?


texts: Why the order of imports matters in a standalone PyQGIS processing script?

texts:
I've come across a weird issue while running standalone PyQGIS processing scripts. The order of imports in the script affects its normal execution.

You can reproduce the problem by opening a Python console and entering the following script (I use GNU/Linux, QGIS 2.6.1, processing plugin v.2.2.0-2, and Python 2.7.3):

# Prepare the environment
import sys
from qgis.core import QgsApplication
from PyQt4.QtGui import QApplication
app = QApplication([])
QgsApplication.setPrefixPath("/usr", True)
QgsApplication.initQgis()

# Prepare processing framework
sys.path.append('/home/YOUR_USER/.qgis2/python/plugins')
from processing.core.Processing import Processing
Processing.initialize()

print Processing.getAlgorithm("qgis:creategrid")

# Exit applications
QgsApplication.exitQgis()
QApplication.exit()

You should obtain:

ALGORITHM: Create grid
    HSPACING <ParameterNumber>
    VSPACING <ParameterNumber>
    WIDTH <ParameterNumber>
    HEIGHT <ParameterNumber>
    CENTERX <ParameterNumber>
    CENTERY <ParameterNumber>
    GRIDTYPE <ParameterSelection>
    CRS <ParameterCrs>
    SAVENAME <OutputVector>

On the other hand, if you switch the order of imports (lines 3 and 4), this way:

from PyQt4.QtGui import QApplication
from qgis.core import QgsApplication

the script now returns... None, because the algorithm was not found.

This issue implies that you cannot run processing algorithms out of QGIS if you (by chance) write imports in the wrong order.

I've checked in StackOverflow but according to Does Python import order matter, the order shouldn't really matter. Moreover, the Style Guide for Python Code tells us to import standard (more generic) libraries first, then related third party libraries, and finally, local application specific imports. I think PyQt4 lies in the 2nd category of imports, whereas PyQGIS would be local application specific, so PyQt4 imports should come first (I'm not an expert on this, though).

Do you have an idea of why this could happen? Have you ever experienced something similar?

EDIT 1: Changed implicit imports (from abc import *) by explicit ones (e.g., from abc import xyz) as suggested by @mike-t.


texts: How to add a point in a empty vetorLayer creadted in memory

texts:
I am working in c++/qt 4.5 and qgis api c++ 2.8.2
I had to create and empty vectorLayer and add Point
I saw in the qgis cookbook that it is possible to create a vector layer in memory ,I did that, and export to a file I did that too and It Works.

Here is , in detailled step , what I want to do :

a- create a vectorlayer in memory

b- add attributes to it

c- add features to it

d- save it in a shp file to be imported in QGIS

I don't know how to add attribute and set value .

is there a specific procedure ?

do you have somme link or some sample ?


texts: Determining amount of forest in overlapping buffers using ArcMap?

texts: I'm using the multiple buffer tool around a series of points for buffer distances 250 m, 500 m and 1000 m around each point. I am then using the tabulate area tool to determine the amount of forest cover from a land cover map within each buffer ring around each point. What I wanted to ensure is that within each buffer zone, such as the 500 m buffer zone that also contains the 250 m buffer zone, is including the amount of forest cover both within the entire 500 m buffer zone area, and is including the forest cover within the embodied 250 m buffer zone as well? I want to make sure it's including the cover in the entire buffer zone for each distance, and isn't excluding what has already been measured in the smaller, overlapping rings.


texts: Calculating Length of Lines within polygon (buffered area) using QGIS?

texts:
I created a buffer around points and then merged these buffers to eliminate overlaps. Now I want to calculate the exact length of tracks within this area. It might be the case that at some points the lines to intersect with the trade area (I need to cut it there).
Any idea on how to succeed?
I already tried the "Sum Line Length" Tool but didn't succeed.


texts:
I wrote the following method as an onClick handler. First and second click, I got result from DB. By the third time, the code stopped in the "getFeatures(trgtFilter)" line and didn't return. In debug mode, I saw that it is waiting for DB connection. Can someone tell me what I did wrong? I'm using GeoTools 15 and Oracle 12.

private Geometry getNewGeometry(String refID) throws Exception {
    if (trgLayer != null) {
        Connection con = null;
        OracleConnection oraCon=null;
        FeatureIterator<SimpleFeature> itr = null;
        try {
            con = ((JDBCDataStore) srcLayer.getFeatureSource().getDataStore()).getConnection(Transaction.AUTO_COMMIT);
            oraCon = (OracleConnection) new DelegatingConnection(con).getInnermostDelegate();

            Filter trgtFilter = editTask.getConfiguration().getReferenceFilter(trgLayer, refID);
            FeatureCollection fc = trgLayer.getFeatureSource().getFeatures(trgtFilter);
            itr = fc.features();
            if (!itr.hasNext())
                return null;
            ...
            CallableStatement cs = con.prepareCall("{? = call MY_FUNC(?,?)}");
            cs.registerOutParameter(1, OracleTypes.STRUCT, "SDO_GEOMETRY");
            cs.setObject(2, trgStrct);
            cs.setObject(3, srcStrct);
            cs.execute();
            ...
        } catch (Exception e) {
            throw e;
        } finally {
            if (itr != null)
                itr.close();
            if (oraCon != null) {
                try {
                    oraCon.close();
                    if (con != null && !con.isClosed())
                        con.close();
                } catch (SQLException e) {
                    LOGGER.error("", e);
                }
            }
        }
    }
}


texts: How can I get ASTER or SRTM DEM for years like 2016, 2020,... in Google Earth Engine? I can't filter date for SRTM DEM in GEE. It seems that SRTM DEM doesn't have multi-year data, it was captured only in 2000 and available in GEE. So what about ASTER DEM? Can I get it for the year like 2020? How can I filter date for it in GEE? (I can filter ImageCollection like Landsat but for images like DEM I can't filter date.)


texts: How are costs returned by the pyQGIS implementation of Dijkstra calculated?

texts:
How are the costs returned by pyQGIS implementation of Dijkstra (QgsGraphAnalyzer.dijkstra) calculated?

As far as I know, these costs correspond to the distance from the starting node to each of the other nodes in the graph, however the cost values I get are way too large to be the distances on the underlying road network.

For example: For a CRS EPSG:4326, the length of an edge in my road network, as returned by feature.geometry().length() is 36.769 metres, where feature.geometry().asPolyline() = [(10573, 5837), (10547, 5863)].

On the other hand, the following is the cost array returned by Dijkstra:

(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, graph.findVertex((10547, 5863)), 0)

cost[graph.findVertex(QgsPoint(10573, 5837))] = 17176150.0812.

Again, from my understanding the cost value above should be 36.769 rather than 17176150.0812.

How can the above be explained?


texts: How to convert an XML to Shapefile using FME?

texts:
I am looking to convert XML to shapefiles using FME.

The XML files are user generated while they draw points and lines (mainly) of interest on a web map. We generate the XML in the background containing the spatial information of user clicks.

I would like to ask whether we need to follow a specific format during XML generation so as to proceed to converting them to shapefiles using FME.

Here is an XML snippet

<?xml version="1.0" encoding="UTF-8"?>
<Configuration>
 <Note>
  <NoteTypes>
   <Type Desc="Point:a" Id="0"/>
   <Type Desc="39.66195294811828" Id="1"/>
   <Type Desc="20.84663271903992" Id="2"/>
  </NoteTypes>
 </Note>
 <Note>
  <NoteTypes>
   <Type Desc="Point:b" Id="0"/>
   <Type Desc="39.6634561263267" Id="1"/>
   <Type Desc="20.84773778915405" Id="2"/>
  </NoteTypes>
 </Note>
 <Note>
  <NoteTypes>
   <Type Desc="Point:c" Id="0"/>
   <Type Desc="39.66210161555193" Id="1"/>
   <Type Desc="20.8490788936615" Id="2"/>
  </NoteTypes>
 </Note>
 <Note>
  <NoteTypes>
   <Type Desc="Line:l1" Id="0"/>
   <Type Desc="39.66206857837209" Id="1"/>
   <Type Desc="20.84651470184326" Id="2"/>
   <Type Desc="39.66291928072009" Id="3"/>
   <Type Desc="20.84694385528564" Id="4"/>
   <Type Desc="39.66319183362618" Id="5"/>
   <Type Desc="20.84708333015442" Id="6"/>
   <Type Desc="39.66352219934388" Id="7"/>
   <Type Desc="20.84763050079346" Id="8"/>
  </NoteTypes>
 </Note>
 <Note>
  <NoteTypes>
   <Type Desc="Line:l2" Id="0"/>
   <Type Desc="39.66354697670904" Id="1"/>
   <Type Desc="20.84759831428528" Id="2"/>
   <Type Desc="39.66295231749322" Id="3"/>
   <Type Desc="20.8476734161377" Id="4"/>
   <Type Desc="39.66248980123196" Id="5"/>
   <Type Desc="20.84869265556335" Id="6"/>
   <Type Desc="39.66220072699667" Id="7"/>
   <Type Desc="20.84896087646484" Id="8"/>
  </NoteTypes>
 </Note>
 <Note>
  <NoteTypes>
   <Type Desc="Line:l3" Id="0"/>
   <Type Desc="39.66203554117646" Id="1"/>
   <Type Desc="20.84653615951538" Id="2"/>
   <Type Desc="39.66141609083402" Id="3"/>
   <Type Desc="20.84785580635071" Id="4"/>
   <Type Desc="39.66219246771504" Id="5"/>
   <Type Desc="20.84893941879272" Id="6"/>
  </NoteTypes>
 </Note>
</Configuration>


texts: Problem importing geometry from CSV in QGIS

texts:
I have a polygon layer (although I have also converted it to points and tried it that way). I would like to copy this layer into excel and add it back in after editing some things. I've attempted to do this by copying the attribute table, pasting it into excel, making my edits, saving the sheet as a CSV, and using "add a delimited text layer" to put it back into my QGIS project. I have not altered the "well-known text" column, and have used that when QGIS asks for the geometry (although I have also tried using the X and Y fields when attempting to do this when the layer is converted to points).

My problem is that the layer appears in the wrong place when I add it back in (the original is in New England, and the version I add back in is off the northern coast of Quebec, when compared to a Google Maps layer imported from the QuickMapServices plugin). The original layer is in the right place, before and after I add what should be its identical twin back in.

Should this be theoretically possible? And/or am I missing a step?

Other things I've tried:

Switching project and/or layer CRSs. I have not found one that fixes the problem - and in any case, the original layer is in the right place, and according to its attribute table, has the same geometry as my copied layer.

Just using a join. The original layer has one point for each geometry - I'd like to have multiple points in some locations (with the other fields differing). I don't see a way to do that outside of excel, other than manually drawing new features, which would be prohibitively time-consuming.


texts: Polygonal shapefile: dissolving adjacent ones

texts:
I would like to solve a issue related to a land use shapefile.

I would like to "merge/dissolve" just the adjacent polygons that have the same attribute. How can I do this using QGIS?


texts:
I'm using ArcGIS API for JavaScript 4.8. My goal is to get the Census tract ID value of the clicked Feature in a FeatureLayer and display the ID value in an alert box. I got stuck when trying to use Query/QueryTask to fetch the value. My code goes here:

require([
  "esri/Map",
  "esri/views/MapView",
  "esri/layers/FeatureLayer",
  "esri/layers/GroupLayer",
  "esri/widgets/LayerList",
  "esri/tasks/support/Query",
  "esri/tasks/QueryTask",
  "dojo/on",
  "dojo/domReady!"
],
function(
  Map, MapView, FeatureLayer, GroupLayer, LayerList, Query, QueryTask, on
) {
    var map = new Map({
      basemap: "topo"
    });

    var view = new MapView({
      center: [-91.05, 30.5161109],
      zoom: 11,
      container: "viewDiv",
      map: map,
    });

    // Generate link to Census tract details
    view.when(function() {
      on(view, "click", displayTractID);
    });

    function displayTractID(event) {
      // TODO
    }

    // Define popup template and get feature layers
    var census2000Template = {
      title: "Census 2000",
      content: [{
        type: "fields",
        fieldInfos: [{
          fieldName: "TRACT",
          label: "Tract",
          visible: true
        }, {
          fieldName: "POPULATION",
          label: "Population",
          visible: true,
          format: {
            digitSeparator: true,
            places: 0
          }
        }]
      }]
    };

    var census2000Layer = new FeatureLayer({
      url: "http://services.arcgis.com/KYvXadMcgf0K1EzK/ArcGIS/rest/services/Census_2000/FeatureServer/0",
      outFields: ["*"],
      popupTemplate: census2000Template
    });

    // Add feature layers to a group layer
    var censusGroupLayer = new GroupLayer({
      title: "Historic Census",
      visible: true,
      visibilityMode: "exclusive",
      layers: [census2000Layer]
    });

    map.add(censusGroupLayer);
  });

I need to fill the TODO part with Query/QueryTask. Any help or sandbox samples?


texts:
I have to use the NLCD 2011 dataset and I am trying to understand the "no change" class. This class is not defined in any of the legend descriptions (Anderson classification).

What is the no change class in NLCD 2011?


texts: Problem with Leaflet Boundary Canvas plugin

texts:
I'm trying to identify the source of a problem that I'm having, I cannot show a tile layer(created with gdal2tiles) inside a polygon generated in Leaflet, using the plug-in BoundaryCanvas(https://github.com/aparshin/leaflet-boundary-canvas)

Some facts:

The tile layer created with Gdal2tiles.py, can be shown in the map with Leaflet if I just use the object tileLayer. From this I understand that the tile layer was generated correctly.
If I use the object tileLayer.BoundaryCanvas(with boundary a polygon), the tile is not shown and in the browser there are errors missing PNG files(it's trying to download png's that are not in the path structure of the tile layer).
If I use the same plug-in with for example the TMS from OSM http://{s}.tile.osm.org/{z}/{x}/{y}.png, it works properly I can show that tile
inside the polygon.

At this point I don't know if the errors comes from the Leaflet plug-in or from the tile layer generation.

I'm using:

gdal2tiles.py --profile=geodetic -z 0-15 -w none ./terra.tif

gdalinfo

Driver: GTiff/GeoTIFF

Files: TerraColor_SanFrancisco_US_15m.tif

Size is 7204 x 7204

Coordinate System is:

GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]

Origin = (-122.809804036593036,38.285436652047451)
Pixel Size = (0.000138800000000,-0.000138800000000)

Metadata:

AREA_OR_POINT=Area

Image Structure Metadata:

INTERLEAVE=PIXEL

Corner Coordinates:

Upper Left  (-122.8098040, 38.2854367) (122d48'35.29"W, 38d17' 7.57"N)
Lower Left  (-122.8098040, 37.2855215) (122d48'35.29"W, 37d17' 7.88"N)
Upper Right (-121.8098888, 38.2854367) (121d48'35.60"W, 38d17' 7.57"N)
Lower Right (-121.8098888, 37.2855215) (121d48'35.60"W, 37d17' 7.88"N)
Center      (-122.3098464,  37.7854791) (122d18'35.45"W, 37d47' 7.72"N)
Band 1 Block=7204x1 Type=Byte, ColorInterp=Red
Band 2 Block=7204x1 Type=Byte, ColorInterp=Green
Band 3 Block=7204x1 Type=Byte, ColorInterp=Blue


texts: Problem with GDAL_WMS QGIS connecting to an ArcGIS Server?

texts:
I've been working, trying to get QGIS (v2.0 32bit) to connect to an ArcGIS Server service using the GDAL_WMS xml file. This is the service I'm trying to connect to:

http://geonb.snb.ca/ArcGIS/rest/services/GeoNB_Basemap_Imagery/MapServer/?SERVICE=WMS&REQUEST=GetCapabilities

This is the xml file I currently have:

<GDAL_WMS>
   <Service name="TMS">
        <Version>1.1.1</Version>
      <ServerUrl>http://geonb.snb.ca/ArcGIS/rest/services/GeoNB_Basemap_Imagery/MapServer/tile/${z}/${y}/${x}.jpg</ServerUrl>
        <Layers>Bathurst_ortho_final_graphic,SDEDO.FTON_ORTHOS_Clip2,MiramichiOrthos,SDEDO.snbortho,Year of Photography</Layers>
        <ImageFormat>image/jpeg</ImageFormat>
   </Service>
   <DataWindow>
        <UpperLeftX>2282341.95</UpperLeftX>
        <UpperLeftY>7698181.3</UpperLeftY>
        <LowerRightX>2737347.05</LowerRightX>
        <LowerRightY>7252058.7</LowerRightY>
        <TileLevel>13</TileLevel>
        <TileCountX>1</TileCountX>
        <TileCountY>1</TileCountY>
        <YOrigin>top</YOrigin>
    </DataWindow>
    <Projection>EPSG:2953</Projection>
    <BlockSizeX>512</BlockSizeX>
    <BlockSizeY>512</BlockSizeY>
    <BandsCount>3</BandsCount>
    <MaxConnections>10</MaxConnections>
    <Cache />
</GDAL_WMS>

And this is what I get:

This is not a WMS issue. The service doesn't support WMS connections and that's the reason for using the XML with GDAL. GDAL then creates a "WMS" like service on top of the images using the XML. It basically fetches images and puts them in the right spot. What I think is happening is it misinterpreting what tiles to fetch and how to stretch them. (Maybe)

Here's an example of a tile fetched directly from the service:
http://geonb.snb.ca/ArcGIS/rest/services/GeoNB_Basemap_Imagery/MapServer/tile/0/64/64.jpg
http://geonb.snb.ca/ArcGIS/rest/services/GeoNB_Basemap_Imagery/MapServer/tile/0/64/64.jpg http://geonb.snb.ca/ArcGIS/rest/services/GeoNB_Basemap_Imagery/MapServer/tile/0/64/64.jpg

This all works so there's something wrong with my XML. The documentation for the XML is here(www.gdal.org/frmt_wms.html), but it's not very descriptive. The specific section that relates to what I want to accomplish is the TMS service.

Any help would be greatly appreciated!


texts: Torque map: dates by year; year in slider text

texts:
I would like to put the year in the slider text and for the data to be scrolled through via year.

Apparently it was working until three weeks ago, and now doesn't. In particular the slider now shows wonky month/day combos rather than just the year. Just the year is what is desired.

How might I fix this?

Initially I fixed the really wonky dates by fixing the -torque-frame-count in the css from default 256 to 227 (for years 1787-2014), which is the number of years. But now, instead of getting random dates, it tends to flip between Jan 1 and Dec 31 on both sides of the date range. It seems disproportionately difficult to find information to get this done - surely this is a common request?


texts: Merge only overlapping shapefiles using FME

texts:
I am very new to FME.

I have 10 shapefiles which overlap but not from beginning to the end! Therefore I want FME to detect "segments" where there is overlapping and splits into such segments the requestor shapefile in order to be able to merge attributes of all overlapping shapefiles into that particular segment.
To simplify I have prepared a very basic sketch to hopefully make clear what I want.

the red shape file down is actually shape file A which is splitted based on overlapping other shapefiles B,C and D


texts: Eliminating large plateaus in coarse resolution DEMs?

texts:
I have a 15-meter, integer DEM I converted into a float raster. However, the original DEM contains large lowland areas that are nearly completely flat. Due to the original resolution and integer format, these flat areas have a step-like appearance with only one meter of difference in elevation despite the large area covered.

How can I correct these large plateaus by smoothing them with a very gentle slope but without touching too much the other features in the landscape?

I am using ArcGIS Desktop but any other free software would be fine.


texts: Is there a common standard for how GIS information is displayed?

texts:
New to GIS programming so this might be a silly question.  Is there a standard way to store/transmit how a GIS feature is to be displayed?  So I can see there's loads of standard ways of storing tha actual GIS data, WKT, WKB, GeoJSON etc, but how should I store display information?
e.g. I've got a polygon, it has a line width, line colour, colour, transparency, etc. How should this information best be stored that's technology neutral and preferably that's well understood by different systems?


texts: geodatabase/workspace content viewer control

texts:
I'm in the process of developing a windows form (desktop) app that will have the ability to display content from user-provided personal geodatabases (.mbd files). I want to create a custom control that will display in some part (a tree view?) all available databases, and in another part (list view?) will display the contents of the selected database. And by content I mean, to display a selectable icon for each dataset in the db.

My question is, is there something similar already implemented?

PS Im using ArcObjects 10, C# 4


texts: Is there a way to hide certain values in a layer on ArcPad 10.0?

texts: My organization uses Arcpad 10 on a trimble nomad. We have many data points in an environmental survey layer, but I would like to only display points from the last two years on Arcpad to help add clarity to the field map. Is there a way to hide the older data on arcpad, while not deleting it from the associated ArcMap document? I have tried creating a new layer from selected points of the desired years on arcmap, then uploading that to the Nomad, but all of the old points are still visible.


texts: Google Earth Engine: ImageCollection (Error) Parameter 'start' is required

texts:
I try to calculate the monthly mean of the Soil Moisture of a point from year 2016 to 2017 from Google Earth. I use a function to calculate the mean by iterating the months since 2016, which is from 1 to 24, but Google earth engine indicates the error:

  ImageCollection (Error) Parameter 'start' is required.

var point = ee.Geometry.Point([30.50,30.56]);

// Create a chart for SoilMoi100_200cm_inst trend in point 1.
var dataset = ee.ImageCollection('NASA/GLDAS/V021/NOAH/G025/T3H')
          .filterBounds(point)
          .filterDate('2016-01-01', '2017-12-01')
          .select('SoilMoi100_200cm_inst');

var months = ee.List.sequence(1, 24);
var start_year = 2016;
var start_date = '2016-01-01';

var byMonthYear = ee.ImageCollection.fromImages(
     months.map(function (m) {
     var start = ee.Date(start_date).advance(m, 'month'); // Starting date

      var y = start_year + Math.floor(m/13);

      var cal_m = (m % 12);
      if (cal_m === 0){
          cal_m=12;
      }
      var d = dataset
        .filter(ee.Filter.calendarRange(y, y, 'year'))
        .filter(ee.Filter.calendarRange(cal_m,cal_m, 'month'))
        .mean()
        .set('month', cal_m).set('year', y)
        .set('system:time_start', start.millis())

      return d;
}).flatten());

print(byMonthYear);