texts: Insert Javascript code in Popup.FramedCloud?

texts:
I would like to have a live stream from a web camera in a Popup.FramedCloud.
All I need to do is to insert this code in the HTML code parameter for popup:

<script src="http://www.live-events.ro/cam-acces-privat-delay.js.php?id=13"
type="text/javascript">

From my tests this is not working.
Please tell me if you tried this with success!

The script is correct, you can see it in action here.

I tried something more simple like "eval('a=1;alert(a)')" but with no success. I see the code in Firebug, but alert window does not appear.

The alternative I think of is to use other object for the popup, an Ext JS object maybe or some kind of reach tooltip library.

Thank you!


texts: Can you use CartoDB.js with maps hosted on the CartoDB website

texts: Is there any way to use the cartoDB.js library with maps hosted on the CartoDB website?  From the documentation it looks like the library can only be used for maps you make on your own website -- is that correct?


texts: Converting raster to vector by generating center lines?

texts:
I have a binary raster containing boundary pixels with the value 1 (white) and no-boundary pixels with the value 0 (black):

I would like to vectorize the raster by keeping one center line trough the boundary pixels. However, most polygonize methods I tried (e.g., gdalogr:polygonize, result in two lines around the outline of the boundary pixels (red lines) or no output (saga:gridskeletonization):

Did I miss any method with which this can be easily done? Or do I have to vectorize the raster first to a vector layer containing two parallel lines and find a way from there to skeletonize/centerline the two lines?


texts: X-Y coordinates unknown coordinate system?

texts:
I received a .dbf file with coordinates in a numerical format I have not been able to identify (for example (X,Y): 6151578.03, 2090875.69). These coordinates represent bus stop locations in the Orange County area (California) and I need to transform them into Lat-Long format (once I figure out what the original format is).

I am using ArcMap 10.2.2.


texts: Finding zero-length LINESTRINGs the wrong way gives funny answer

texts:
I tried

select st_astext(geom) from table where st_pointn(geom, 1) = st_pointn(geom, 2)

But I get some results which don't match the criteria. For example: LINESTRING(-117.1642086 34.2779308,-117.1642145 34.2779305)

Why do I get the wrong answer for that query?

(I have since changed my query to use st_length).


texts: What do the heatmap layer settings numbers in mapbox stands for?

texts:
I have layer settings of mapbox heatmap as follows from this official site:

What do these numbers stand for (the first, the second, the third, the fourth number) and what are their limitations please? (min, max)

In general there is written what the settings do, but how to customize them? I haven't found any detailed documentation link so far, just saw some basic results.


texts: Regroup multiple polygon with same attributes in one

texts:
After looking for an answer on some forums (without success :/) i post here in hope on of you can help me.

I have a polygon shapefile with thousands of entities (=row) and >10 attribute column.

A lot of these polygons share the exact same attributes, in fact there is only 12 different types of polygons.

I want to regroup all polygons with sames attributes in only one (so I just want to have 12 entities at the end and keep all the attribute data).

I have tried lots of options "regroup", "dissolve", "dissolve by attribute", etc... but each time one of these issues occurs :

I lost lot of columns (just keeping 1 or 2/3);
I got just one entity;
Some polygons "disappear"

The closest I found was SAGA "Dissolve", but it don't process on all polygons.

Any idea ?

I'm using QGIS 2.14.


texts: How to append the .csv data into GeoJSON using Python

texts:
The .CSV file has 6 indents, we need to append all the indent with GeoJSON file and updated file ready for open in QGIS

    import json
    import copy
    import csv
    import sys

    #print(sys.argv)

    def read_in_csv(csvfile):
    with open (csvfile.csv, 'rb') as f:
        reader = csv.reader(f)
        header = (f.next()).strip('\n').strip('\t')
        properties = list(reader)
        flatlist_properties = [x for y in properties for x in y]
    return flatlist_properties, header

    def append_properties_to_geojson(csvfile, geojsonfile):
        data = json.loads(open(geojsonfile).read())
        data = copy.deepcopy(data)
        values_to_add = read_in_csv(csvfile)[0]
        try:
            for entry in range(len(data["features"])):
            data["features"][entry]["properties"][read_in_csv(csvfile)[1]] = values_to_add[entry]
    except:
        print ("Note: Your CSV and geojson features do not contain the same number of values")
    return data

    if __name__ == "__main__":
    csvfile = sys.argv[1]
    geojsonfile = sys.argv[2]

    with open('updated_data.json', 'w') as outfile:
        json.dump(append_properties_to_geojson(csvfile, geojsonfile), outfile, indent=6)

    print ("Data dumped into 'updated_data.json'")

# Usage: python add_csv_to_geojson.py csvfile.csv geojsonfile.json

Error show

    ['add_csv_to_geojson.py', 'csvfile.csv', 'geojsonfile.geojson']
    Traceback (most recent call last):
    File "add_csv_to_geojson.py", line 37, in <module>
    json.dump(append_properties_to_geojson(csvfile, geojsonfile), outfile,  indent=4)
    File "add_csv_to_geojson.py", line 22, in append_properties_to_geojson
    data = json.loads(open(geojsonfile).read())
    File "C:\Program Files\Python37\lib\json\__init__.py", line 348, in loads
    return _default_decoder.decode(s)
    File "C:\Program Files\Python37\lib\json\decoder.py", line 337, in decode
    obj, end = self.raw_decode(s, idx=_w(s, 0).end())
    File "C:\Program Files\Python37\lib\json\decoder.py", line 355, in  raw_decode
    raise JSONDecodeError("Expecting value", s, err.value) from None
    json.decoder.JSONDecodeError: Expecting value: line 1 column 1 (char 0)


texts:
Does anyone familiar with the QGIS Zonal Statistics tool know why stats calculations might return empty?

Inputs

DEM for South Carolina: http://dnr.sc.gov/GIS/descdem.html
A vector shapefile: https://s3.amazonaws.com/www.improvz.com/gis/Archive.zip

Desired result: Zonal Statistics calculates elevation on all shapes.

Actual result: Zonal Statistics plugin returns empty values for all fields

QGIS 2.8 on MacOSX.


texts: How to use latitude/longitude with projection in R

texts:
I am using the SABS dataset:
http://nces.ed.gov/surveys/sdds/sabs/

In particular, let's say I am using the "SABS_1314_High" layer. (This provides the school attendance boundaries for many US school districts.) The identifying variable is "ncessch".

I also have a list of latitude and longitude coordinates that looks like this ("sampledata.csv"):

num lat long
1   37.378309   -122.052739
2   37.320722   -122.020981
3   37.436397   -122.149727

For each row of this, I want to know the "ncessch" value of the polygon that the point belongs in.

After much trial and error (I am completely new to GIS, though not to R), I am able to read in the shapefile properly, together with its projection:

library(rgdal)
poly3 <- readOGR (".", "SABS_1314_High")
proj4string(poly3)
[1] "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

I can also create a SpatialPointsDataFrame for my list of points:

library(sp)
locations <- read.csv("sampledata_2.csv", header = TRUE)
locations <- SpatialPointsDataFrame(locations[ , 2:3], locations)

At this point, it would appear that I need to use the over function in the sp package to check, for each point, which polygon it lies in, but that throws an error:

over(locations, poly3)
Error: identicalCRS(x, y) is not TRUE

It would therefore appear that the projections are not the same. Okay, so I tried a couple of different things:

library(mapproj)
locations.proj <- mapproject(locations[,2], locations[,3], projection = proj4string(poly3))
Error in min(x, na.rm = na.rm) : invalid 'type' (list) of argument
In addition: Warning message:
In is.na(x) : is.na() applied to non-(list or vector) of type 'S4'

Using a different package:

library(proj4)
locations.proj <- project(locations[,2:3], proj = proj4string(poly3))
Error in as.numeric(x) :
  cannot coerce type 'S4' to vector of type 'double'

What am I missing?


texts: How to read all raster values in an array?

texts:
I want to read the every cell value in a raster file. When I run the script, it should give me the cell location and it's value.

I can reach center point of the raster and read it's value with this script:

rasterLyr = QgsRasterLayer("C:\qgis_data\SatImage\SatImage.tif",
                                "Sat Image")
rasterLyr.isValid()

#extent means QgsRectangle
c = rasterLyr.extent()

qry = rasterLyr.dataProvider().identify(c,QgsRaster.IdentifyFormatValue)

qry.isValid()
qry.results()

If there is a toolbox which solves this problem in QGIS 3.0, it also helps me.

I use QGIS 3.0 as I mentioned and python 3.6


texts:
I'm a longtime GIS user who's just recently started learning more about data structures, cursors, and search/sort algos. I've used arcpy.mapping's UpdateCursor and the like for years now and didn't think much about how they worked. Recently, one of my assignments was to load a dictionary into a linked list in C - that meant rolling nodes with pointers to the "next" node's address and manually indexing the cursor to search something.
Is the underlying data structure for a shapefile's anything like this - is it implemented using a linked-list-like structure? Or a different data structure entirely? Also, does each row in the shapefile database containing a shape string and some other fields (id, name...) include information (like an address) regarding how to find the "next" or "previous" row?


texts: How to join LINESTRINGS in a specific order in PostGIS?

texts:
I have a path from point A to point B, which I want to double back on itself to create a path from A to B back to A.

It seemed like using ST_LineMerge on the MULTILINESTRING created by ST_Collect from the two LINESTRINGS was the answer, however it doesn't appear to maintain the ordering of the paths supplied in the MULTILINESTRING.

For example in psudosql

ST_LineMerge(
    MULTILINESTRING(
        (Point A, Point B),
        (Point B, Point A)
    )
)

Creates a path B to A to B, as opposed to A to B to A as I would expect based on the LINESTRING order and directions.

How can I reliably join these paths together?

Here is a concrete example...

SELECT ST_AsText(ST_LineMerge(ST_GeomFromText('MULTILINESTRING Z ((151.2569403 -33.4921765 14, 151.2568375 -33.4922183 14),(151.2568375 -33.4922183 14, 151.2569403 -33.4921765 14))')));
SELECT ST_AsText(ST_LineMerge(ST_GeomFromText('MULTILINESTRING Z ((151.2568375 -33.4922183 14, 151.2569403 -33.4921765 14), (151.2569403 -33.4921765 14, 151.2568375 -33.4922183 14))')));

which both return,

LINESTRING Z (151.2568375 -33.4922183 14,151.2569403 -33.4921765 14,151.2568375 -33.4922183 14)


texts: How to stylize polygons sharing the same attributes to represent a common outline

texts:
In QGIS 3.14.16, I want to outline multiple polygons which share same attribute value.
Styling the layer with the same attribute value bring the following result:

The style I wish to have is:

I could use the dissolve function to get the expected result. However, this is not a solution for me, as I need to modify frequently the attribute table.
Is there a way to configure such a style?
My final goal is to style polygons layer with:

Outlined polygons in function of a field (as described above),
Colored polygons in function of another field.

Graphically speaking, I want to jump from that style:

to this:

Polygons are colored with category according to field A. They are also outlined with category style according to field B.
EDIT
Here are examples of what I tried so far:
Example 1:

Example 2:

Example 3:

Example 4 (this last example is the best way I found to circumvent that issue, see comments below):


texts: Is it possible to create Tessellation table in schema other than sysspatial

texts:
As per geoserver teradata documentation, tessellation table is required in sys spatial schema only.

My question is - Is it possible to configure geoserver to pick the tessellation table from a different schema? Has anyone tried it with success?

I tried changing the tessellation table schema name while creating the datastore, but still geoserver looks for this table in syssptial schema only.

Any inputs?


texts: Changing visibility of layers on project with OpenLayers and GeoServer?

texts:
I'm working on a proyect with openlayers and geoserver, the idea is to have a select menu where the user can choose which layer is visible, something like this

  <select id="format">
    <option value="rios_dobles">Rios Dobles </option>
    <option value="rios_simples">Rios Simples</option>
    <option value="lagos_lagunas">Lagos Lagunas</option>
    <option value="vias">Vias</option>
    <option value="limite_parroquial" selected>Limite Parroquial</option>
    <option value="basicos">Basicos</option>
    <option value="Centros_Poblados">Centros Poblados</option>
    <option value="Aptitud_final">Aptitud Final</option>
    <option value="uso_final">Uso Final</option>
  </select>

The layers are set like this

var limites = new ol.layer.Image({
                source: new ol.source.ImageWMS({
                    ratio: 1,

                    url: 'http://localhost:8080/geoserver/gis2017/wms',
                    params: {'FORMAT': format,
                        'VERSION': '1.1.1',
                        STYLES: '',
                        LAYERS: 'gis2017:limite_parroquial',
                    }
                })
            });

And are displayed on the map like this

var map = new ol.Map({
                controls: ol.control.defaults({
                    attributionOptions: /** @type {olx.control.AttributionOptions} */ ({
                        collapsible: false
                    })
                }).extend([mousePositionControl]).extend([scaleLineControl]),
                interactions: ol.interaction.defaults().extend([
                new ol.interaction.DragRotateAndZoom()
                ]),

                target: 'map',

                layers: [
                    limites,
                    //lagos,
                    //vias,
                    //riosDobles,
                    aptitud,
                    //centros,
                    //riosSimples,
                    uso,
                ],

                view: new ol.View({
                    projection: projection
                }),
            });

At least as a test I've seen if I can make the layer limites visible or not, but it won't work, I've tried using the methots setVisible and setVisibility but regarldess of weather or not I put true or false, the layer won't display (or any other for that matter)

I don't know what to change or do, I'm completely stuck.


texts: ArcMap hangs on "Initializing Application..." when starting from Visual Studio?

texts:
This has worked just fine for years; I don't know what changed.  ArcMap 10.0 starts just fine on its own, it's just when using the debugger to test my add-in that it won't start.  The add-in, once compiled and installed, runs just fine.

Possibly related is an issue where in Windows Explorer, if I delete a file, it will not disappear unless I hit 'refresh'. This started happening around the same time that ArcMap 10.0 wouldn't start from Visual Studio.


texts: How can you know the sequence/names of the bands of a GeoTIFF file in Java?

texts:
I want to know the sequence of the bands (4 and more) of a GeoTIFF file for example whether the sequence is blue, green, red, near infra red, temperature, etc or green, blue, ...

How can I know that in Java, for example in GeoTools library? And by the way, are there specific tags in the GeoTIFF that indicates that (like TIFFTAG_IMAGEDESCRIPTION) or it is up to the image reader software to interpret the bands (even if it assumes the sequence mistakenly and considers blue instead of red and so forth).


texts: Opengeo : Add new server dialog doesn't show up

texts:
While trying to follow the tutorial of Opengeo for creating a new application, I am getting stuck on one of the first steps: namely the addition of a new server dialog.
(http://suite.opengeo.org/docs/apps/clientsdkintro/viewer/addnewserver.html)

I added the following lines in my app.js and restarted my server application. Unfortunately I cannot get the "Add a New Server ..." line as described in the tutorial when adding new layers.

* @require overrides/override-ext-ajax.js
  ...
proxy="proxy/?url="

The full content of my app.js is available at https://gist.github.com/anonymous/6899291

Note that I am starting my application using the following command on a Linux machine where Opengeo was installed : suite-sdk debug -g http://localhost:8080/geoserver myapp

Thanks


texts:
I have two polygon shapefiles that I need to merge. I applied 'Merge vector layers' tool from Data management in QGIS 3.16, but it took too long to process and I realized a message indicating that one shp has integer field type and the other integer64.
Assuming this is the problem why the process could not be completed, how can I change the data type of one of these layers? And which one should it be?


texts: Creating file geodatabases in multiple folders using iterator?

texts:
Software Spec:  ArcMap 10.2.1  Advanced Licence.

I have a single folder, let's call it "Processing".
Inside this folder are 15 empty folders.  Let's call these "Output_1" through to "Output_15". There is nothing else in the "Processing" folder.

I am trying to use ModelBuilder to create an empty file geodatabase inside each Output folder with the same name as the folder.
i.e. in the "Output_1" folder, the file gdb will be called "Output_1.gdb".

I thought this would be very straight forward using an "Iterate Workspace" iterator, where the input is the "Processing" folder.

I could then use the output workspace from the iterator as the file gdb location into the "Create file gdb" tool, and the %Name% from the iterator as the file dgb name.

However, the "Create File GDB" tool won't accept the workspace output from the iterator as the file gdb location.

I have tried changing the iterator settings so that the workspace type it looks for is only folders, but this doesn't fix this issue.

I have also tried changing the input variable (called "Processing") between a "folder" and a "workspace".  This makes no difference.

Is there something fundamental I am missing?  Or is this a bug?

If I first run just the iterator part of the model, I am then able to connect the "wokspace" output as the input location for the "Create file gdb" tool.  I can then run just the "Create file gdb" tool, but obviously this only runs for the last folder that was read by the iterator.


texts:
I want to smooth Lines in postgis, to be more specific, I want to smooth road data from OpenStreetMap.

Basicly I thought of something like this scenario.

This would be a possible road segment.

I want to use ST_LineInterpolate_Point() to select the middle Points

After this, I want to create New Points between the original ones and the darkgreen helpline. (I have no Idea how to achieve this...)

I read about circularlinestrings. I hope it will be possible to lay one of these over the darkgreen and red points.

Finaly I want to use ST_LineInterpolate_Point() again on the circularlinestring (if this is possible) and create some new points(yellow ones), to make a new normal LineString.

To sum up, I am new to PostGIS, but since I haven't found anything helpful for LineSmoothing I need a function where I can input a LineString and the number of yellow points per segment. The function should return a new LineString.

First of all: Is this possible to achieve with just pl/pgsql?

Second, since I am new to all of this, can somebody give me some advice where to start?


texts: ogr2ogr write multiple layers to one geopackage

texts:
I'm using ogr2ogr to write a SQL query from PostGIS to a GeoPackage in a linux shell.
Query by Query to many GPKG:

ogr2ogr -f "GPKG" point.gpkg PG:"host=localhost dbname=data user=user password=*******" -sql "select * from point-features"
ogr2ogr -f "GPKG" line.gpkg PG:"host=localhost dbname=data user=user password=*******" -sql "select * from line-features"
ogr2ogr -f "GPKG" polygon.gpkg PG:"host=localhost dbname=data user=user password=*******" -sql "select * from polygon-features"

I would like to write multiple layers from queries into one GeoPackage, but I cannot find a documentation to this.

Like:

ogr2ogr -f "GPKG" pointlinepolygon.gpkg pointlayer linelayer polygonlayer PG:"host=localhost dbname=data user=user password=*******" -sql "select * from pointlinepolygon-features"

Would be happy about support.


texts: Change QGIS background colour using python

texts:
I would like to change the background colour of a QGIS project using python.

I can change the colour using the GUI via the Project Properties dialog:

However, I cannot find the equivalent setting in the Python API. I have tried setting the canvas colour using iface.mapCanvas().setCanvasColor(), but the colour change does not stay after the project file is saved, and reopened again layer.


texts: How to compose WFS GetFeature link from WMS capabilities definition?

texts:
This publication of geoinformation provides download links for GetCapabilities documents via WMS. The platform also runs a WFS service as can been seen here.

I wonder if it is possbile to compose a download link for ShapeFiles/GML or another format which can be converted ... from the given WMS information. I managed to display the street tree information in QGIS and could load details for individual features. However, this is pixel data and I would like to download vector data if possible.

The second GetCapabilities document defines layers named pflanzstandorte2011 and pflanzstandorte_aktuell. It would be nice if these were accessible via the platform's WFS API.


texts:
I am having a problem opening a shapefile using QGIS using 'add vector later'. This is the error message I am getting:

/Users/Conrad/AmerindianLands 2/Amerindian_Land_Karasabai.shp is not a valid or recognized data source

Can anyone help?

I have made sure that it is being extracted from a folder containing the other necessary files with the following extensions (.dbf,.prj,.shp, .sbx,.sbn).


texts:
I naively assumed that the rasterType() function would give me the data value type for the raster (i.e. integer or float), but it just tells you whether it's a singleband or multiband raster.

Is there a function in QGIS that I can use to check whether the input raster to a script is integer or float (similar to what you get out of the VALUETYPE in the GetRasterProperties function in ArcGIS)?

Basically need to run a check so that if the input isn't an integer raster it should throw an exception and exit.

ArcPy snippet I am trying to convert:

EvidenceLayer = parameters[0].valueAsText
    valuetype = gp.GetRasterProperties (EvidenceLayer, 'VALUETYPE')
    #if valuetype not integer:
    if valuetype > 8:
        gp.adderror('Not an integer-type raster')
        raise ErrorExit


texts:
I have a map where i'm loading several WMS layers from different urls.  I'm trying to figure out a way to turn each layer on/off by an identifier, but i can't figure out how to add a unique identifier to the layer, or source of the layer.  Basically i have:

for (const entry of entries) {
    const tileWMSOptions = {
              url: entry.URL,
              params: { LAYERS: layerNames, TILED: true },
            };
            const layer = new TileLayer({
              source: new TileWMS(tileWMSOptions),
            });
            mapObj.addLayer(layer);
}

my 'entry' object has an ID that i would like to associate with this layer, so that i can do some sort of

map.getLayers();

run through them, and do something specific if the id of that layer matches what i'm looking for.  I don't see any way of doing this from https://openlayers.org/en/latest/apidoc/module-ol_layer_Tile-TileLayer.html or https://openlayers.org/en/latest/apidoc/module-ol_source_TileWMS.html .  Any ideas?


texts: Get a list of all the images that match the conditions in Google Earth Engine (not only the first)?

texts:
Normally when I run this code I get the least cloudy image but I want to know what the other images are.

I use this code:

    var image = ee.ImageCollection('COPERNICUS/S2')
        .filterBounds(PuntoAbast2)
        .filterDate('2018-07-01','2018-07-02')
        .sort('CLOUD_COVER', false);
    var mejorImagen = ee.Image(image.sort('CLOUD_COVER').first());

With the set dates I get only one image but If I extend the date range I will get more, but this code only will return one image.


texts: MapBox import from TileMill gives only 1 layer

texts:
This is my first time using MapBox and I'm trying to import a map from TileMill.

The map has 11 layers, in TileMill, but when I import it directly to MapBox, it only seems to have a single layer.

It still shows all of the layers, but if I run the following, it returns a single layer, which seems to have all of the data on:

mapbox.auto('map', 'MAP_NAME', function(map){
        console.log(map.getLayers());
});

For the map functionality, I need to be able to switch the layers on/off via user interaction, so having everything in a single layer is a slight issue.

Am I doing something wrong? I really hope that I don't have to export each of my TileMill layers individually, as this will take so much time..

Any help/advice would be much appreciated.

EDIT

After hours of searching, I've found this question on the MapBox support pages: http://support.mapbox.com/discussions/tilemill/2122-creating-a-batch-of-layers

It looks like TileMill is unable to export the layers and keep them as separate layers.
Instead I will probably have to export each layer individually, which is a bit of a pain, but oh well.

If anyone has any better suggestion, please let me know :)


texts:
In QGIS, how can I output a map

84 in. x 72 in. at 300 dpi

with QGIS 2.12.3 and Windows 10?

In print composer, I can dimension the canvas and set dpi, but it seems that I cannot add my map at that size (FAR from it!).

Can it be done in print composer? If not in print composer, then how?


texts: Converting point to polygon in a given distance

texts:
I'm working in QGIS 2.18 but I can work with QGIS 3 or ArcMap 10.2.

I have a multipoint layer providing from a lasDataset where I've isolated the buildings.

I want to create a polygon layer where every building or group of buildings are individual polygons. The problem is I can't use any attribute to group them because I don't have any. Concave Hull only creates one polygon. So I have to group them by other way. by distance? How can I do it?

The points I have

How I want to group them


texts: List Map Services and Publish/Delete Map Services using Python in ArcGIS 10.1

texts:
I looked at the link
http://resources.arcgis.com/en/help/main/10.1/index.html#//00s30000004s000000

which describes listing of map service for 9.3, how can we do this in 10.1, moreover I am trying to start stop, publish services , I read on ArcGIS help that the new ArcGIS API Administration uses the RESTful Architecture wherein an http request is made to start/stop a service. Is there any other way? Can we delete or publish/republish a service? Any help on this is most welcome. Thanks


texts: Divide Line Into Sections Using Points Offset from It in ArcGIS Desktop?

texts:
I have track (line) and milemarkers (point). Each point is offset by a different distance. Is there a way to find where each point would intersect the line if going perpendicularly out from the line?

As an analogy, imagine riding a train, sitting perpendicular to the direction of travel, looking out the window. Each time you see a milepost directly in front of you, you make note of your position on the track.


texts: Automatically creating user CRS from postgis in QGIS

texts: I have my custom CRS in table "spatial_ref_sys" in database. When I load layer from schema "public", QGIS create this CRS. But if layer from other schema, it's not. What I need to solve this problem?


texts:
I tried QTILES plugin for QGIS and it works very well.
I can export tiles in .mbtiles or zip format.
Now I'm trying to know how call these tiles from leaflet in this way:

L.tileLayer('http://www.cityplanner.it/experiment_host/supply/tiles/car_tiles/{z}/{x}/{y}.png', {
    maxZoom: 18,
    attribution: '..'
}).addTo(map);

(http://www.cityplanner.it/tilemill-alla-mappa-online/)

I don't know if is ok .mbtiles format or folder structure from .zip.

How leaflet reads .mbtiles?


texts: How to Substract one image band from the subsequent in Google Earth Engine?

texts:
i am trying to find the day of year (doy) for each year with the highest change in NDVI compared to the day before, to map in the end a shift in start of season.

My code to seperate each year so far (returns a list with 18 ImageCollections with unkown elements):

var months = ee.List.sequence(0, 17).map(function(n) {
  var start = ee.Date('2001-01-01').advance(n, 'year');
  var end = start.advance(1, 'year');
    return ee.ImageCollection("MODIS/006/MOD13A1").select("NDVI").filterDate(start, end);
});

Subsequently i would like to substract each ndvi-image from the previous one and compare get the doy of the biggest result. Anyone any Ideas how to do that?


texts: Deleting overlap between two polyline layers in ArcMap?

texts:
I have two stream layers and they overlap, I want to delete this overlap on one layer so I do not repeat any information.

How do I do this?

I can only seem to find information regarding polygons.


texts: Function access to surrounding object's variable

texts:
var CreateFeature = OpenLayers.Class(OpenLayers.Control, {

    clazz : "xyz",

    set_item : function(e) {

        var position = map.getLonLatFromPixel(e.xy);
        var latitude = position.lat;
        var longitude = position.lon;

        var marker = new OpenLayers.Feature.Vector(new OpenLayers.Geometry.Point(longitude, latitude), {
            draggable : true,
            clazz : this.clazz, /** this one should be set to 'abc' **/
        }, items_symbolizer);

        layer_items.addFeatures([marker]);
        marker.state = OpenLayers.State.INSERT;

    });

var br = new CreateFeature(layer_items, {
        clazz : "abc",
        title : "Create Feature",
    });

I want to create an Object with a function parameter that has access to a variable of the object's scope. However, this.clazz remains undefined. Any ideas?

I am looking for a way to either set this.clazz to 'abc' or to 'xyz'.


texts: PostGIS treating unique geometries as equal

texts:
I have two different geometries (albeit very closeby) that Postgres is treating as equal.

db=> select '0101000020E6100000C5BCD8C5FEF45DC0AFE29EF87B584740'::geometry = '0101000020E61000005A7CAFC6FEF45DC01A23C8F77B584740'::geometry;

 ?column?
----------
 t
(1 row)

Can I dial-up PostGIS equality resolution so that these are unequal?


texts:
Okay, I've been banging my head against the wall on this and I'm pretty sure I'm missing something silly. So I'm hoping you guys can help.

The idea here is to allow the user to type in the label they want. It appears on the map in the center and then they can move the label to the location they want it to be. (Then I'll add functionality for styles, but basic functionality first, frills later.) The problem is that I can't move the labels. I can't select them either. In fact, I can't select anything on the layer, but it is on top (or at least it should be). And yes indeed, 'labelSelect' is set to True.

Anyway, here is the code to initialize the layer and the controls:

editableLayer = new OpenLayers.Layer.Vector("User Layer");
editableLayer.setZIndex(4000);
layerArray.push(editableLayer);

map = new OpenLayers.Map({
div: "map",
layers: layerArray,
controls: [
    new OpenLayers.Control.LayerSwitcher(),
    new OpenLayers.Control.PanZoomBar(),
    new OpenLayers.Control.Navigation(),
    new OpenLayers.Control.ScaleLine({geodesic: true}),
    new OpenLayers.Control.KeyboardDefaults(),
    new OpenLayers.Control.MousePosition({
        displayProjection: new OpenLayers.Projection("EPSG:4326"),
        numDigits: 5,
        emptyString: '',
        prefix: "x: ",
        separator: " y: "
    }),
    new OpenLayers.Control.EditingToolbar(editableLayer),
    new OpenLayers.Control.SelectFeature(editableLayer, {toggle: true, clickout: true}),
    new OpenLayers.Control.DragFeature(editableLayer),
    overviewMapControl
]
});

Here is the code for creating the label:

var labelElement = document.getElementById("LabelText");
var labeltext = labelElement.innerHTML;

var labelFeature = new OpenLayers.Feature.Vector(
            new OpenLayers.Geometry.Point(
                map.getCenter().lon,
                map.getCenter().lat
            ),
            {
                text: labeltext
            },
            {
                label: labeltext,
                labelSelect: true
            }
        );
editableLayer.addFeatures([labelFeature]);

So, what did I mess up? I tried ModifyFeature with the mode set to DRAG as well, that didn't work either. I suspect the problem is somehow in my layer ordering, but even when I add code for map.setLayerIndex it doesn't help and when I inspect the elements in the browser it seems to be on top.

Thanks.


texts: Passing an ICursor to ArcObjects ITableSort

texts:
I'm trying to pass a cursor (from a querydef with a prefixclause of DISTINCT) to itablesort.cursor. I've read another post on GIS SE on the same subject that warned against doing this. That post is 4 years old -- has anything changed since then? My code works without using the cursor but it is slow when working with larger datasets. Hence why I want to use the querydef.

Here's excerpt of the code. The GFLayer at the beginning is an iGeoFeatureLayer pointing to a file GDB featureclass.  If I comment out the line, .Cursor = pTblCursor, everything runs without error.

Dim pDataset As IDataset = pGFLayer.FeatureClass
Dim pPosCursor As ICursor = Nothing
Dim pTblCursor As ICursor

Dim pWS As IFeatureWorkspace = pDataset.Workspace
Dim pQueryDef As IQueryDef2 = pWS.CreateQueryDef
pTblSort = New TableSort

If pGFLayer.FeatureClass.FindField("SYM") = -1 Then
    MsgBox(pGFLayer.Name + " does not have a SYM field.  Layer will not be processed.")
    pLayer = pEnumLayer.Next
    Continue While
End If
With pQueryDef
    .PrefixClause = "DISTINCT"
    .Tables = pDataset.Name
    .SubFields = "POS,SYM"
    .WhereClause = ""
    pTblCursor = .Evaluate
End With

If Not pTblCursor Is Nothing Then
    With pTblSort
        .Cursor = pTblCursor
        .Table = pGFLayer.FeatureClass
        .Fields = "POS"
        .Ascending("POS") = True
    End With
    pTblSort.Sort(Nothing)
    pPosCursor = pTblSort.Rows
End If


texts: feature.get(0).getinfo()[0] does not initialize inside function GEE

texts:
I have a list of features which I am going to use map.function to extract specific value from features but inside function when I am using .getinfo()[0] to return first value give me an error: "Internal error: function argument not initialized."
Let me know if anyone is familiar with this problem.

Here is the sample code:
var map_m = function(i) {
  i=ee.Number(i)
  var feature = ALL.get(i).getInfo()[0]
  var ff=feature

  return ff
}


texts: Avoid internal patches when clipping a raster in R based on a polygon

texts:
I want to clip a raster in R based on a polygon. The analog function in arcgis that does the same is Data Management Tools -> Raster -> Raster Processing -> Clip.

What I have until now is:

library(raster)

set.seed(123)

r <- raster()
res(r) <- c( 0.016666, 0.016666 )

extent <- shapefile("colombia_2000.shp")

r <- crop(r, extent, snap='near')
r <- setValues(r, sample(x=c(0,1),size=ncell(r), replace = T))

undata_mask <- mask(r, extent)

You can find the "extent" file in the following link:

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

The polygon (extent) I am using represents the area between 1000 an 2000 meters above see level in the Andean Mountain range in Colombia, and as you can see there are some internal patches that are not covered with the polygon, which belongs to the mountain peaks that exceed this altitud range.

When I perform the crop and mask functions in R, those internal patches remain untouched as if they belonged to the polygon. When I perform the calculation in ArcGis, the internal patches are indeed clipped, which is what I want to achieve in R.

How Can I get rid of this internal patches?

Here are some images of What I just said.

The first one is with R and the second one is with ArcGIS


texts: Calling individual bands from multiband.tif in python perform action and save output in mutilband.tif?

texts:
I am trying to write a script that would go through my directory and open each *.tif file, evaluate each of their 120 bands and write the results out to another 120 multiband *.tif file. The only way I could figure out how to do this was to iterate through the input directory,

open each raster,

iterate through its bands,

evaluate each band for the threshold value (2),

if condition is met assign a number if condition is not met assign a value of zero,

save each band to an output.tif file

then use the composite function to collect all 120 output.tifs

But I'm stuck. When I try to iterate through the bands I get an error message (see below).

First is there an easier way to do this without creating 120 files for each of my 350 multiband input.tifs?
If not, can you help me figure out how to iterate through each band?

This is the error message:

  ========= RESTART: C:/Python27/ArcGIS10.5/Scripts/StartJD.py ========= about to start cycling through rasters in the raster cycle at count 0

  Traceback (most recent call last):   File
  "C:/Python27/ArcGIS10.5/Scripts/sapflow_StartJD.py", line 36, in

      sapRas = Raster(inRas + "\Band_" + count) #get individual bands for input raster   File "C:\Program Files
  (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Functions.py", line 4329, in
  Plus
      in_raster_or_constant2)   File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Utils.py", line 53, in swapper
      result = wrapper(*args, **kwargs)   File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Functions.py", line 4326, in
  Wrapper
      ["Plus", in_raster_or_constant1, in_raster_or_constant2]) RuntimeError: ERROR 000732: Input Raster: Dataset \Band_ does not
  exist or is not supported

This is an example of the input file names within the directory, 26_sapsum_2011_converted (note no extension is listed but the accompanying file does have the extension 26_sapsum_2011_converted.tif.aux).

Here's my script:

# objective of this code is to assign the Julian day when sapsum = 2 and
write a multiband output file for each year
# Import arcpy module
import arcpy, os, sys
#import the environment settings and spatial analyst extension
from arcpy import env
arcpy.env.extent = "MAXOF"
# check out spatial analyst extension
if arcpy.CheckExtension('Spatial') == 'Available':
    arcpy.CheckOutExtension('Spatial')
else:
    arcpy.AddMessage("Error: Couldn't get Spatial Analyst extenstion,
exiting")
sys.exit(1)

from arcpy.sa import *

# Set the input/output directories
env.workspace = r"C:\Netcfd\BCSD\sapsum\sapsumtiff\Year\26_Y_20s"

# set number of bands
number_of_bands = 120
# Get a list of Rasters
rasters = arcpy.ListRasters('*.tif')

print 'about to start cycling through rasters'

for raster in rasters:
  # cycle through rasters in directory
 inRas = Raster(raster) #open and rename raster input
 count = 0
 print "in the raster cycle at count", count
 for band_number in range(0, number_of_bands):
      count = count + 1
      sapRas = Raster(inRas + "\Band_" + count) #get individual bands for
      # input raster
      print 'in the band number cycle and the current band number is', count
      output = con((sapRas == 2),count,0) #if sapRas is equal to 2 assign
      #count if SapRas not equal to 2 assign a value of 0
      #save output in temporary file
      Rasname, Rasext = os.path.splitext(inRas)# sparate the name from the
        #extension
      save.output(env.workspace + "/JD" + Rasname + "_comp.tif")
      for i in range(0, number_of_bands):
           listcomp = [] #create a list of output tifs
           for ras in arcpy.ListRasters('_comp.tif'):
                listcomp.append(ras)
           arcpy.CompositeBands_management(listcomp, env.workspace +
          "Composite" + inRas + '.tif') #collect output tifs into a
                                        #composite
print "Done! and look it's beer o'clock"

raw_input("Press <enter>") # this keeps the window open until you press
                            "enter"


texts: Connect and modify map through localhost in ubuntu (content management system)

texts: How can i set cartaro in linux,facing problem with installation instruction given on site.


texts: Can QGIS remember the state of open tables?

texts:
Is there an option to choose while saving a QGIS project,to save the state of the table views and their location on the screen?

I am on 3.6.3-Noosa


texts:
I'm selecting all records and then removing certain values from this selection in a Python script. I've succeeded using a for loop in a SearchCursor to remove records from the selection but it is very slow. Is there anyway to remove the records all at once using "IN" like how this would work in the field calculator

SELECT * FROM FeatureClass WHERE Field IN ('a', 'f', 'l')

This is what I have so far:

gRoutes =  ('R10T', 'R10TBP', 'R10TWP', 'R20T', 'R20TBP', 'R20TWP', 'R32T', 'R32TBP', 'R32TWP', 'R45T', 'R45TBY', 'R45TBP', 'R45TWP', 'R64T', 'R64TBP', 'R64TWP', 'R96T', 'R96TBP', 'R96TWP', 'R10TBY', 'R20TBY', 'R32TBY', 'R64TBY', 'R96TBY')

# field name
QF1 = "Product_Code"

WC1 = '"' + QF1 + '" IN ' + "'" + gRoutes + "'"

arcpy.SelectLayerByAttribute_management("layerFC","REMOVE_FROM_SELECTION", WC1)

This is what the WC1 expression prints as:

'"Product_Code" IN \'[\'R10T\', \'R10TBP\', \'R10TWP\', \'R20T\', \'R20TBP\', \'R20TWP\', \'R32T\', \'R32TBP\', \'R32TWP\', \'R45T\', \'R45TBY\', \'R45TBP\', \'R45TWP\', \'R64T\', \'R64TBP\', \'R64TWP\', \'R96T\', \'R96TBP\', \'R96TWP\', \'R10TBY\', \'R20TBY\', \'R32TBY\', \'R64TBY\', \'R96TBY\']\''

And this is the error I'm getting:

Start Time: Thu Oct 20 14:59:48 2016
ERROR 000358: Invalid expression
Failed to execute (SelectLayerByAttribute).
Failed at Thu Oct 20 14:59:48 2016 (Elapsed Time: 0.02 seconds)


texts: QGIS 3.4 Stream delineation problem, lines running too far into ocean

texts:
I am using QGIS 3.4.4 Madeira version trying to delineate the stream and catchment areas of a small Bay in Honduras based on the technique of Hanz van der kwast tutorial "Stream and Catchment delineation in QGIS 3".
I imported DEM > reprojected > Used SAGAs Fill Sinks (Wang and Liu) > USED SAGAs Strahler Order to delineate streams.

What happens is that the tool cannot decipher what is land at sea level and what is the actual sea and creates channels running straight off into the ocean. How can I correct the DEM before using stream delineation tool so I don't have my streams running far off into the ocean? I want them to stop when they hit their drainage point (i.e. Bay/Ocean).

I've already tried clipping the raster to a polygon of the land area to simply cut out the ocean from the raster, but then the results don't work as I am guessing I cut out important information that the tool needs to process (streams and basins are completely missing).


texts: How to resample a batch of rasters using OGR/GDAL?

texts: I'm looking for a function similar to GRASS r.resample. I have a folder full of TIF files and want to reduce their resolution ("scale them down") to reduce file size.


texts:
I am new to QGIS and trying to import two thousand geo pngs to create a raster map. each png is indexed 2-4 colours.

I have tried several things:

importing one by one works, as does to drag-and-drop all of them at once, but i have to confirm the CRS for every file which is not practical. is there a way to say 'yes to all' ?
i've created a VRT, but when that is imported, the files all render black, except the first one that renders red. i've tried playing around with the vrt styles, but it seems that is meant to use the colors from only first 1-3 bands, and in my VRT every file is a separate band?
I also have separate layers as b&w geo-tiffs, that were used to create the coloured pngs. however importing them has exactly the same issues as A+B above, either I have to manually create each one, or the map ends up black.

I guess i don't understand something basic about raster rendering in gdal.


texts: Mapping Commuting Zones from Counties for US 1990 Data in ArcMap?

texts:
This question follows Seeking 1990 US Commuting Zone shapefiles?

I am attempting to map US commuting zones as defined in 1990. I have successfully mapped counties from a county shapefile obtained here. After completing the following advice, I get the error message below:

  In the xls, there is a column called CZ90 which I assume is the 1990
  commuting zone code. Once you've joined XLS to FIPS shapefile, use the
  dissolve tool with the CZ90 field as the dissolve field to output a
  shapefile/feature class with one polygon per commuting zone.

I used the dissolve tool and checked CZ90 as my only field; however, I get the following error message. It appears that I correctly joined the CZ90 excel file with the FIPS shapefile.

How can I successfully dissolve the counties into commuting zones?


texts:
I am trying to create an intersection between two different layers of two different geodatabases using ogr2ogr, but without success.

I have been using the following command line prompt:

ogr2ogr.exe intersection.shp .\Sea.gdb .\Roads.gdb -dialect sqlite -sql "SELECT ST_INTERSECTION(Sea_Depth.geometry,Roads_100m.geometry) as geometry FROM Sea a,Roads b"

But I get the error:

Warning 1: layer names ignored in combination with -sql.
ERROR 1: In ExecuteSQL(): sqlite3_prepare_v2(SELECT ST_INTERSECTION(Sea_Depths.geometry, Roads_100m.geometry) as geometry FROM Sea a,Roads b):
no such table: Sea

For ease of context, the sea layer is made of polygons and the road layer is lines. Intersection.shp will ideally be the resulting file.


texts: How to smooth large vector polygons from Raster

texts:
I have a landuse classification  raster file (6 classes) of a large area. But I need it in vector polygons so I vectorised it. No problem so far, but the result is - of course - that every pixel is exactly vertorised with all its egdes :) But I need it smooth like contours. So, okay, I could generalize it with Boyle's Algorithm for instance, but its just too much for my computer to calculate that...

So what could I do? Thanks!
Martin


texts: How to Export Raster with same classification colors seen in ArcMap?

texts:
I have a single band raster, but i like to export it with the same colors as in Arcmap.

Any help please.


texts: Map-matching GPS points to the road network!

texts:
I am quite new to Postgres and PostGIS.

I am trying to do a simple map-matching (or path-inference if you wish) of some GPS points. I am using postgres(postgis).

I have a table in my database containing my GPS points:
gpspoints (Lat, Long, time, speed, azimuth, geometry, ...)

I have imported a shapefile of the road network into my database as well:
routes (gid, idrte,version, nomrte, norte, clsrte, geom, ...) -- I DO NOT HAVE the azimuth of the links, and I am not sure how can I calculate it using the postgis function ST_azimuth. All I have is a shapefile of the road network containing columns enumerated above.

I want to associate each GPS point to the nearest link (in a buffer of 20 meters around the point) only if the direction of the link agrees with the azimuth of the GPS point (+ or -15 degrees) and retrieve the projected position. Otherwise I want it to search for the NEXT nearest link, within the 20 meter buffer, which has an acceptable azimuth! (Just like the picture!)

I want the new coordinates of the projected GPS points to be added in the table "gpspoints" as "projectedLat" and "projectedLong".

(In the picture below, actual points are demonstrated using a direction while projected points do not have any direction mark)


texts:
I have been successfully adding text delimited layers until today. Now the program freezes if I try this, and I have to press escape on my keyboard. Everything else works fine. I have just downloaded the new version (2.12.3), but have the same problem.

Does anyone know if there is a simple solution?


texts: How to generate a Global Voronoi diagram?

texts:
Using arcMap, QGIS and similar software it is easy to create a non-global Voronoi diagram.

The resulting shapefile polygons however are showing weird behavior when crossing the -180 meridian or the poles. What tool can I use to create a globally consistent Voronoi diagram similar to this website? A code based solution (Python/Javascript) is also ok.

Data used to create Voronoi:

+----------+--------------+-------------+--+
| name     | x            | y           |  |
+----------+--------------+-------------+  +
| Miami    | -80.25551469 | 25.79555801 |  |
+----------+--------------+-------------+  +
| Boston   | -71.07082218 | 42.3585578  |  |
+----------+--------------+-------------+  +
| Honolulu | -157.8555857 | 21.30573449 |  |
+----------+--------------+-------------+--+


texts: Call external programmes from python plugin in QGIS on windows

texts:
I need to call an external programme from python plugin in qgis, I didn't use os.systeme because I want to hide the launched programme with 'startupinfo.dwFlags |= subprocess._subprocess.STARTF_USESHOWWINDOW'.

import subprocess
import os
startupinfo = None
if os.name == 'nt':
    startupinfo = subprocess.STARTUPINFO()
    startupinfo.dwFlags |= subprocess._subprocess.STARTF_USESHOWWINDOW
subprocess.Popen("TeighaFileConverter",shell=True, stdin=None, stdout=PIPE, stderr=PIPE,startupinfo=startupinfo)

it's workin on pyzo but not working in QGIS editor.
this code is working on linux but not on windows :/


texts:
I have a ArcGIS JavaScript API map. On map click I want to display a circle and highlight all the points within that circle.

However, when I create the circle and zoom in the map, the circle zooms in as well. How can I get the circle to stay the same pixel radius no matter what the zoom level of the map?

    <script>
        require([
            "esri/Map",
            "esri/layers/GeoJSONLayer",
            "esri/views/MapView",
            "esri/widgets/Zoom",
            "esri/core/watchUtils",
            "esri/geometry/geometryEngine"
        ], function (Map, GeoJSONLayer, MapView, Zoom, watchUtils, {buffer}) {
            const url =
                    "https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_month.geojson";

            const template = {
                title: "Earthquake Info",
                content: "Magnitude {mag} {type} hit {place} on {time}",
                fieldInfos: [
                    {
                        fieldName: "time",
                        format: {
                            dateFormat: "short-date-short-time"
                        }
                    }
                ]
            };

            const renderer = {
                type: "simple",
                field: "mag",
                symbol: {
                    type: "simple-marker",
                    color: "orange",
                    outline: {
                        color: "white"
                    }
                },
                visualVariables: [
                    {
                        type: "size",
                        field: "mag",
                        stops: [
                            {
                                value: 2.5,
                                size: "4px"
                            },
                            {
                                value: 8,
                                size: "40px"
                            }
                        ]
                    }
                ]
            };

            const geojsonLayer = new GeoJSONLayer({
                url: url,
                copyright: "USGS Earthquakes",
                popupTemplate: template,
                renderer: renderer //optional
            });

            const map = new Map({
                basemap: "gray",
                layers: [geojsonLayer]
            });

            const view = new MapView({
                container: "viewDiv",
                center: [-168, 46],
                zoom: 3,
                map: map
            });

            let csvLayerView;
            let highlight;
            let resultFeatures;

            view.whenLayerView(geojsonLayer).then(layerView => {
                csvLayerView = layerView;
            });

            view.on("click", ({mapPoint}) => {
                const geom = buffer(mapPoint, 500, "miles");

                view.graphics.removeAll();
                view.graphics.add({
                    geometry: geom,
                    symbol: {
                        type: "simple-fill",
                        color: [51, 51, 204, 0.2],
                        style: "solid",
                        outline: {
                            color: "white",
                            width: 1
                        }
                    }
                });

                const query = geojsonLayer.createQuery();
                query.geometry = geom;
                csvLayerView.queryObjectIds(query).then(oids => {
                    if (highlight) {
                        highlight.remove();
                    }
                    highlight = csvLayerView.highlight(oids);
                    return csvLayerView.queryFeatures(query);
                }).then(({features}) => {
                    resultFeatures = features;
                });
            });

        });
    </script>


texts:
I have multiple .asc grid files, I need to merge these into one raster grid file.

I would like to know the procedure to merge .asc files.


texts: Reproject raster without projection info

texts:
I'm pretty new to gdal and got a presumably easy question for experienced gdal users. My problem: I have a raster ascii file of which I don't know the projection or epsg code. However, I want this raster to be projected in epsg 32663 (WGS84).

The information I get from gdalinfo contains the coordinate origin and the corner coordinates:

I tried directly assigning a projection with gdal_translate, but ofc it didn't work as I don't know what else (and where) I have to provide additional information (by it didn't work I mean the spatial position is entirely wrong when comparing with a correctly projected raster in epsg 32663).

I also tried gdal_translate with providing gsp's (my corner coordinates; don't know if this is even legitimate), or gdalwarp - no success.


texts:
I am using some legacy code which requires my data in PostgreSQL box type e.g. this format (37357,50073),(34242,357).

However, I am using PostGIS with box2df i.e. I get: BOX2DF(588893 4515997, 600342 4517442.5). Is there a PostGIS/PostgreSQL function which can transform my data from the BOX2DF format to the BOX format?

I cannot simply do a type cast:

  ERROR:  cannot cast type box2df to box


texts: Vector Tiles generating taking time in Geoserver

texts:
I have displayed map box vector tiles for EPSG:3857 projection customized by me with 21 zoom levels.below is code using in my application

<!DOCTYPE html -->

  <title>Vector tiles</title>

  <link rel="stylesheet" href="https://cdn.rawgit.com/openlayers/openlayers.github.io/master/en/v5.3.0/css/ol.css">

vector tiles

`

but tiles displaying like as show in below image

Those empty space taking time to generate.I tried to seed the layer but getting error Chunk [format=application/x-protobuf;type=mapbox-vector] is not a valid entry .
How to speed up layer to generate vector tiles and what is issue with seeding


texts: FME Server is not able to connect to ESRI enterprise geodatabase and shows no error

texts:
I have an environment with a FME Server 2018.1.2 Build 18590 and ArcGIS Server 10.6.1.
It has been running 2 pretty similar workspaces with no problem (A and B).
The ArcGIS server was upgraded to 10.7.1. After that the following behavior is observed:
Workspace A Fails with this log:
 Creating reader for format: Esri Geodatabase (ArcSDE Geodb)
 Trying to find a DYNAMIC plugin for reader named `GEODATABASE_SDE'
 Loaded module 'GEODATABASE_SDE' from file 'C:\Program Files\FMEServer\Server\fme\GEODATABASE9.dll'
 FME API version of module 'GEODATABASE_SDE' matches current internal version (3.8 20180604)
 Performing query against GEODATABASE_SDE dataset `<Unused>'
 Opening the Geodatabase reader

Here the log stops.
Then, Workspace B is triggered, and shows the following in the same part as the Workspace A
 Creating reader for format: Esri Geodatabase (ArcSDE Geodb)
 Trying to find a DYNAMIC plugin for reader named `GEODATABASE_SDE'
 Loaded module 'GEODATABASE_SDE' from file 'C:\Program Files\FMEServer\Server\fme\GEODATABASE9.dll'
 FME API version of module 'GEODATABASE_SDE' matches current internal version (3.8 20180604)
 Performing query against GEODATABASE_SDE dataset `<Unused>'
 Opening the Geodatabase reader
 FME has checked out an Esri license. The product checked out is 'ArcServer'
 Installed ArcGIS version is '10.7'
 Connection made to server 'MYSERVER' for dataset 'UPDM' (transactional version 'sde.DEFAULT') using instance 'sde:sqlserver:MYSERVER' and user ID 'myuser'
 ArcSDE release: '10.0'. Underlying database: 'SQL Server'
 R_3 Reader: Using Enhanced geometry ......
.....

Log continues and finish successfully.
Here is the tricky :
If Workspace A is triggered again, it works fine.
After several times (3-4) start failing again.
Then I run the Workspace B and the Workspace A works again.
I noticed that it breaks just prior to the activation of the license. always is in the same location (First Feature reader with format Esri Geodatabase ArcSDE Geodb).
I know that this article states that I should be using FME 2019.1 but the thing that caught my attention is the Workspace A works, and fail once in while. moreover the Workspace B does not fail.


texts:
How to do after a search, zoom in on the value found on the map? My code is:

$("#txtSearch").autocomplete({
    source: setoresComerciais.features.map(function(d){
        return d.properties.sco_num_sc + " - " + d.properties.sco_dsc_loc
    }),
    select: function(event, ui){
        map.fitBounds(stComerciaisLayer.getBounds(ui.item.value));
    }
});

First, I get return d.properties.sco_num_sc + " - " + d.properties.sco_dsc_loc, this return, for example: 1 - FORTALEZA
After, I use select, where have the param ui, when do ui.item.value return 1 - FORTALEZA, so this:

map.fitBounds(stComerciaisLayer.getBounds(ui.item.value));

Where:

stComerciaisLayer = L.geoJSON(setoresComerciais, {
                style: function (feature) {
                    return feature.properties && feature.properties.style;
                },
                onEachFeature: onEachFeature,
(...)

It was to return exactly what I wanted, the result of the search and the zoom in it, but don't work. What did I do wrong?


texts: How to highlight multiple cities/towns in Google Maps?

texts:
I'm trying to create a service area map in Google Maps. When one searches for a city or town, Google outlines the town in a light red. If possible, I'd like to have multiple towns selected in this method.

I'm aware its possible to draw it out using the line tool, but there are about 60 towns and it would take forever. Any ideas?


texts:
I'm currently writing a python plugin where I would like user with 0 experience of GIS to use QGIS to view data geographically anyhow,

I'm running two processing functions:

import processing
processing.runandload('qgis:reprojectlayer', str(self.input_file_path) + "shapefiles/temp.shp",'EPSG:4326', file_name + '.shp')
processing.runandload("grass:v.voronoi", file_name+".shp",False, True, extent, -1, 0, 3, file_name_poly)

Both of them are adding the new layer in qgis "Layers Panel", can I turn this off? (and just saving the files in the selected path). Now I got one layer called Reprojected and one called Voronoi Polygons in the Layer Panel.


texts: Assigning QGIS Attribute Table commands, like zoom to selected feature, to shortcut key?

texts:
Is there a way to assign QGIS Attribute Table commands, like zoom to selected feature, etc a shortkey?

If I try to assign them a shortkey using setting >> shortkey. Attribute table commands doesn't show up.


texts: I have an ArcPy script tool that exports a map with appropriate labels (10.1 SP1) and my coworker who made the original mxd that I'm basing it on is using Myriad Pro as the font for the text elements. How might I get this font in my mxd? Or could I simply get her to create the mxd which the tool will reference?


texts: Extracting Transit Layer polylines (e.g. subway lines) and using them in Google My Maps?

texts:
I am trying to create a Google My Map of Tokyo Metro and Toei Subway lines. So far, I've been painstakingly drawing the subway lines manually, only using a mouse. Here's what I've got so far (WIP, of course). It is extremely labor intensive to do it this way, and I am convinced there must be a way to automate this process, even if it requires utilizing the Google Maps API or even Google Earth Pro.

Here's another example of the kind of thing I want to do, though obviously not at this scale: https://japonyol.net/editor/ajax.html

Does Google Maps allow you to take a polyline object from the Transit Layer (e.g. the polyline of a subway line) and add the object to a Google My Map?


texts: How to zoom exactly to a desired scale in Flex?

texts:
I am using ArcGIS API for flex 2.4 with ArcGIS server 10. I am trying to write the widget where user can define the map scale (type it in text box) and it will zoom to that scale. I successfully created a widget. But When I type scale in the text box it still zooms to the nearest LOD scale level. e.g. if i zoom to 5000 scale, it will take it as 4,514.

Could anyone suggest what change I need to make? Is there anything that needs to be defined in MapManager?

<?xml version="1.0" encoding="utf-8"?>
<!--
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2010 ESRI
//
// All rights reserved under the copyright laws of the United States.
// You may freely redistribute and use this software, with or
// without modification, provided you include the original copyright
// and use restrictions.  See use restrictions in the file:
// <install location>/License.txt
//
////////////////////////////////////////////////////////////////////////////////
-->
<viewer:BaseWidget xmlns:fx="http://ns.adobe.com/mxml/2009"
                   xmlns:s="library://ns.adobe.com/flex/spark"
                   xmlns:mx="library://ns.adobe.com/flex/mx"
                   xmlns:esri="http://www.esri.com/2008/ags"
                   xmlns:viewer="com.esri.viewer.*"
                   widgetConfigLoaded="init()">
    <fx:Script>
        <![CDATA[
            //this function called when the widget's configuration is loaded
            private function init():void
            {
                if (configXML) // checking for valid content in the configuration file
                {

                }

                mapScaleObj.text = String(map.scale);
            }

            private function doScale():void
            {
                map.scale = Number(mapScaleObj.text);
                mapScaleObj.text = String(mapScaleObj.text);
            }
        ]]>
    </fx:Script>
    <viewer:WidgetTemplate id="helloWorld"
                           width="250" height="110">
        <viewer:layout>
            <s:VerticalLayout horizontalAlign="center" verticalAlign="middle"/>
        </viewer:layout>
        <s:SkinnableContainer >
            <s:HGroup>
                <s:Label height="100%" color="#FFFFFF" fontSize="14" fontWeight="bold" text="Scale 1:"
                         verticalAlign="middle"/>
                <s:TextInput id="mapScaleObj"   text="{map.scale}"
                             restrict="0-9"
                             enter="doScale();"
                             />
            </s:HGroup>
        </s:SkinnableContainer>
    </viewer:WidgetTemplate>
</viewer:BaseWidget>


texts: How to batch process (100s) of interpolations in QGIS or other open source software programs?

texts:
Is there an efficient way to run 100s of interpolations in QGIS or other open source software programs?

I have groundwater contamination data collected from monitoring wells at about 800 different sites.  I would like to interpolate a set of groundwater plume contours for each site.  I'd like to avoid running 800 separate interpolations.  Any ideas?


texts: Can I use GeoWebCache with GeoServer Layers which have different native and declared SRS?

texts:
We are using GeoServer to serve raster tiles with a native SRS of 27700, but are forcing use of a declared SRS of 900913.

When I try to seed these tiles using GeoWebCache through the GeoServer web admin interface no tiles are cached. However when I use the native SRS of 27700 seeding works fine.

Is what I'm trying to do possible? And if so, how do I do it?

Thanks,


texts:
A few weeks ago we crashed our server by trying to tile a large image service (600gbs) (multiple times without clearing the cache) which caused us to exceed the servers free disk space.

We fixed the issue by deleting the service and clearing the server cache.  From this point forward we have not been able to publish and tile an image service, however we are able to publish a regular map service.

Does anyone know what might be causing us to not be able to publish/tile image services?

Is there a potential break in code somewhere in a config file on the server?

We are running ArcGIS 10.1 and we have substantial free disk space.


texts: If a person stands on the north, or south pole. How do they determine which way is east/west?


texts: Use and difference of "Routing and Analysis metric" and "Radial metric" in SDNA Integral Network

texts: I want to know what it means to select Euclidean Metrics and Angular Metrics in SDNA Integral Network between "Routing and Analysis Metrics" and "Radial Metrics ".


texts: How to create group boxes in a QGIS custom form using Qt Designer?

texts:
I am trying to create a custom form for QGIS using Qt Designer.
Let's say I have a field which uses relational values. I want these values to be checkboxes. I can do it easily with the included form designer in QGIS.
Now I am trying to get this with Qt Designer using checkboxes grouped together in a GroupBox widget. But I can't manage to get the values back to QGIS.

The form designer in QGIS.

The custom form made from the form designer using relational values.

The custom form I would like to create in Qt Designer.


texts:
The following OverpassTurbo query for train routes in Sicily, Italy, returns ten's of tiny tracks:

[out:json][timeout:25];

//Sicily: 39152
rel(39152);map_to_area -> .searchArea;

(
  node["railway"="rail"](area.searchArea);
  way["railway"="rail"](area.searchArea);
);

out body;
>;
out skel qt;

I was wondering if there were a way to somehow merge them so that I end up with much fewer?

The goal is to get a GPX file that contains all the train routes displayed in the Transport Map available on OSM.org:

Actually, it's a lot more than "ten's of tracks": The query above returns… 4,158 occurences of "<trk>".

--

Edit: Drawing the routes manually isn't very smart, but it does the job


texts:
I am writing a set of plugins doing complementary functionality. I want all the plugins under the same sub menu under Vector. This works fine, but when I am reloading the plugin during development, it will not remove the old menu item, it will just add on a new one, so that I end up with a lot of menu items for the same function, although only the last one is working.

So when I reload the plugin doing the "Load spectral data" part, I do not want to do anything on the menu for Altitude correction. (Which btw shows exactly the same when that is the one I am working on)

The relevant parts of my plugin:

metadata.txt:category=Vector

DataLoader.py (__init__) : self.menu = self.tr(u'&Spectral data')

DataLoader.py (initgui):
    self.add_action(
        icon_path,
        text=self.tr(u'Load spectral data'),
        callback=self.run,
        parent=self.iface.mainWindow())

DataLoader.py (unload):
    for action in self.actions:
        self.iface.removePluginVectorMenu(
            self.tr(u'Load spectral data'),
            action)

What should I add to tell it to remove "Load spectral data" under the "Spectral data" submenu?


texts: GROUP BY timestamp interval 10 minutes PostgreSQL


texts: OpenStreetMap uMap show classified features

texts:
For the first time ,I am using the OpenStreetMap uMap to share my data online (http://umap.openstreetmap.fr/)
I am trying to show classified features (different colors in the same layer) but I can't find how to do it. all my features are shown with the some color.
is there any solution? if not , is there any other web application that can enable me to do it ?


texts:
I am trying to add a SOC machine to SOM server (distributed setup). From the server, I can see the service running. I have all the accounts properly setup. I may be still missing a super critical settings somewhere. This is the error that I am getting:

  Machine is not a valid server container. Machine
   is not found.

My environment is Windows Server 2008 R2 (both SOM and SOC machines), ArcGIS Server 10.

I am curious if some had any idea what else to try. Is there a way to check if SOC is properly configured without using ArcCatalog or ArcGIS Server Manager.

Thanks
Jay


texts: Coordinate system problem when importing ESRI, MapInfo, kmz and CSV files into GRASS

texts:
I'm very new to GRASS.

I've imported some government files (from http://listdata.thelist.tas.gov.au/opendata/) that are available in ESRI and MapInfo format and they display successfully, and show coordinates in GDA94 / MGA zone 55, as expected.

I've also imported a CSV file with points in decimal degrees long and lat using:

v.in.ascii input="stops.txt" output=metro_stops format=point separator=comma skip=1 x=6 y=5 z=0 --overwrite

and I've imported a kmz file from Google Earth.

The points from the CSV file and the kmz file display with coordinates in decimal degrees, so they don't appear anywhere near the ESRI data.

I realize my problem is to do with telling GRASS the coordinate system of the data, or transforming it somehow, but I don't know where to look.


texts:
Working on being able to zoom to the extent of features in a search row to produce individual maps for row (like data driven pages but for individual rows on the dataset).

I just can't seem to get it to zoom to extent?  I have looked around but can't find what I need in the syntax help guide.

Script below

import arcpy

from arcpy import env
arcpy.env.overwriteOutput = True
#Set File workspace
file_workspace = 'N:\\GIS\\Projects\\AA_Leith_Hawkins_TestBed\\Search_Cursor\\Search_Cursor.gdb'
env.workspace = file_workspace

#buffer set up
Holdings = 'N:\GIS\Projects\AA_Leith_Hawkins_TestBed\Search_Cursor\Search_Cursor.gdb\Data\Holdings'
distances = [1000, 4000]
unit = "Meters"
#Make a feature layer
arcpy.MakeFeatureLayer_management(Holdings,'Holdings_Layer')

#searchcursor for evey row in dataset
with arcpy.da.SearchCursor(Holdings, ['Holding_Reference_Number'])as Holdings_Ref_cursor:
    for row in Holdings_Ref_cursor:
    print row[0]
    query = "Holding_Reference_Number = " + str(row[0])
    print query
    File_output = file_workspace+ '\\' 'Buffer_'+str(row[0])
    #print File_output

    #Select Feature using the reference number from holdings layer
    arcpy.SelectLayerByAttribute_management('Holdings_Layer', 'NEW_SELECTION',"Holding_Reference_Number = " + str(row[0]))

    # Export holding to geodatabase
    Holding_Boundary = file_workspace+ '\\' 'Holding_'+str(row[0])
    arcpy.management.CopyFeatures('Holdings_Layer', Holding_Boundary)

    #Mutliple ring Buffer using Selected Features
    arcpy.MultipleRingBuffer_analysis('Holdings_Layer', File_output, distances, unit, "", "ALL")
    arcpy.MakeFeatureLayer_management(File_output, 'Buffer_Layer')
   #arcpy.Buffer_analysis("Holdings_Layer", ofc, var_Buffer, "FULL", "ROUND", "ALL", "")
    print 'Buffer complete'

    #Intersect Features
    Intersect_out_feaatures = file_workspace+ '\\' 'Intersect_'+str(row[0])
    arcpy.Intersect_analysis([Holdings,File_output],Intersect_out_feaatures, "", "", "INPUT")
    print "intersect Complete"

    #add Layers to the Map
    mxd = arcpy.mapping.MapDocument('N:\\GIS\Projects\\AA_Leith_Hawkins_TestBed\\Search_Cursor\\Search_Cursor_mxd.mxd')
    df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
    addLayer = arcpy.mapping.Layer(File_output)
    addLayer2 = arcpy.mapping.Layer(Holding_Boundary)
    addLayer3 = arcpy.mapping.Layer(Intersect_out_feaatures)
    arcpy.mapping.AddLayer(df, addLayer, "TOP")
    print addLayer
    arcpy.mapping.AddLayer(df, addLayer2, "TOP")
    print addLayer2
    arcpy.mapping.AddLayer(df, addLayer3, "TOP")
    print addLayer3

    #zoom to layer

    print df
    lyr = arcpy.mapping.ListLayers(mxd, '', df)[2].name
    print lyr
    extent = lyr.getExtent()
    print extent
    df.extent = extent

    #Export Map to PNG File
    Png_output = "N:\\GIS\Projects\\AA_Leith_Hawkins_TestBed\\" + str(row[0]) + '.png'
    arcpy.mapping.ExportToPNG(mxd,Png_output)
    print 'Map Created'
    del mxd

The error message I get is (I have included the prints from the run as well so you can see what they produce to try and make it easier:

108050022
Holding_Reference_Number = 108050022
Buffer complete
intersect Complete
Buffer_108050022
Holding_108050022
Intersect_108050022
<geoprocessing Data Frame object object at 0x1262EC00>
Buffer_108050022
Traceback (most recent call last):
  File "N:\GIS\Projects\AA_Leith_Hawkins_TestBed\Search_Cursor\Script\SearchCUrsor.py", line 65, in <module>
extent = lyr.getExtent()
AttributeError: 'unicode' object has no attribute 'getExtent'


texts: Using CartoDB to change projection of maps - then exporting as SVG

texts:
This seems like a fairly simple question, but I've looked everywhere for the answer and can't find it.

I've imported geoJSON UK constituency data into CartoDB and it looks great. I now want to export this as an SVG to open in Adobe Illustrator, however when I do this my map looks like it has the wrong projection on it. (i.e The map is warped and I want it to appear flat.)

What can I do to fix this? Is there something I can do in CartoDB to get it to export flat?


texts: I am dynamically creating a join between a GeoFeatureLayer and an external ITable derived from a SQL query and opened with OpenQueryClass(). It works. Now if I could only delete the Join!? But even worse than that: when I re-enter the code that creates a new join (after I run some code that seems to remove the existing joined fields from the GeoFeatureLayer) it fails to create the ITable the 2nd time because "The FeatureClass already exists." So it appears that not only must I reset the GeoFeatureLayer I must also remove the named ITable.  Is this correct?  If so, how?


texts: How can I disperse points within polygons using ArcGIS?

texts: I have huge number of points. These points are within polygons. I want to disperse the points in the polygons randomly. I tried to use the disperse point extension in Arcview, but it was not useful for me. Disperse Markers tool in ArcGIS toolbox was also not useful and can’t create permanent disperse points, it just creates a view of point disperse. Are there any tools or scripts to disperse  points in polygons?


texts: geoserver support to proxy other tiling services?

texts:
I'd like to use a central server to provide a common location for all of our mapping / imagery resources.  Our data is primarily in the form of geotiffs with a variety of base layers from various sources on the web.

Geoserver, given a set of geotiffs, does a nice job of serving up tiles (via WMTS or TMS).  However, I don't know of a way to easily proxy other external sources of data like the sample base layers from http://services.arcgisonline.com/ArcGIS/rest/services.  Unfortunately, it doesn't seem like many services implement a WMS layer (or don't make it public).  Currently we have logic in our application to special case various web services but it is a bit of hack.  Ideally I'd like to present a common web services interface to all our tiling providers.  This would enable WMS clients like world wind applications to "just work".

My only thought right now is to build a small python web service that would expose a list of capabilities via WMS and then add this service into geoserver as a cascading WMS store.  This seems like a lot of work to simply serve up map tile URLs.

I am curious if this this problem has been solved before or if there is an easier solution.

Thanks!


texts:
It seems I am missing simple functionality in Pro. How to update Contents pane in ArcGIS Pro other than to remove a layer from it and add it again?

I am working on some script which adds a subtype for a specific layer and it works fine, however, Contents pane is not automatically updated and subtype is not visible in it. There is a workaround where the script will remove a layer from the map(which might remove styling) and add it again.

But is there a simple refresh button for Contents pane or something similar I am missing here?


texts: What is the best way to extract the geometry of a road from OSM data? The goal is to display segments of highway that are closed due to weather. The weather data consists of positions (in KM measured along the road) from major cities, and length of the affected stretch of road. The main problem is clipping the road based on these distances; afterwards rendering seems straightforward.


texts: Can Server Object Extensions filter map data and thus modify the returned map image? (ArcGIS Server 10)

texts: I was wondering if I could use a Server object Extension (SOE) to modify the image data that is returned by a Map Service. By that I mean, for each incoming http request (and based on the authenticated user), set a filter on the Server Object's map layer(s) that will be reflected in the map image. Is this possible and would there be performance issues? Also would there be performance issues of doing a database query inside an SOE?


texts:
Is it possible to download Google Maps data points?

It would be useful in order to gather all the restaurants in a dataset for example.


texts:
I've been trying to access a several WMS servers in order to visualize imagery on my QGIS but when I put the URL on the searcher no results are given, as if there were no WMS server with such URL . I've tried WMS URLs that have worked before, but I can't find them now. Nevertheless, those WMS servers with whom I had established connection prior this incident I haven't had any issues.

That's what I get after doing the search, basically nothing:

Any ideas?


texts: Does the origin of Sentinel-2 L1C images on AWS represent a point at the center of the pixel or a point at the upper left corner?

texts:
Does the origin of Sentinel-2 L1C images on AWS represent a point at the center of the pixel or the upper left corner of the pixel? i.e. AREA_OR_POINT=Point, whether AREA_OR_POINT=Area means that a pixel value should be assumed to represent a sampling over the region of the pixel and AREA_OR_POINT=Point means a point sample at the center of the pixel.
Here's an example of 'asoc' box from a Sentinel-2 image

<XMLContent>
  <?xml version="1.0" encoding="UTF-8" standalone="no"?>
  <gml:FeatureCollection xmlns:gml="http://www.opengis.net/gml" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengeospatial.net/gml http://schemas.opengis.net/gml/3.1.1/profiles/gmlJP2Profile/1.0.0/gmlJP2Profile.xsd">
    <gml:boundedBy>
      <gml:Null>withheld</gml:Null>
    </gml:boundedBy>
    <gml:featureMember>
      <gml:FeatureCollection>
        <gml:featureMember>
          <gml:RectifiedGridCoverage dimension="2" gml:id="RGC0001">
            <gml:rectifiedGridDomain>
              <gml:RectifiedGrid dimension="2">
                <gml:limits>
                  <gml:GridEnvelope>
                    <gml:low>1 1</gml:low>
                    <gml:high>10980 10980</gml:high>
                  </gml:GridEnvelope>
                </gml:limits>
                <gml:axisName>x</gml:axisName>
                <gml:axisName>y</gml:axisName>
                <gml:origin>
                  <gml:Point gml:id="P0001" srsName="urn:ogc:def:crs:EPSG::32723">
                    <gml:pos>399965 8700035</gml:pos>
                  </gml:Point>
                </gml:origin>
                <gml:offsetVector srsName="urn:ogc:def:crs:EPSG::32723">10 0</gml:offsetVector>
                <gml:offsetVector srsName="urn:ogc:def:crs:EPSG::32723">0 -10</gml:offsetVector>
              </gml:RectifiedGrid>
            </gml:rectifiedGridDomain>
            <gml:rangeSet>
              <gml:File>
                <gml:fileName>gmljp2://codestream/0</gml:fileName>
                <gml:fileStructure>Record Interleaved</gml:fileStructure>
              </gml:File>
            </gml:rangeSet>
          </gml:RectifiedGridCoverage>
        </gml:featureMember>
      </gml:FeatureCollection>
    </gml:featureMember>
  </gml:FeatureCollection>
</XMLContent>


texts: Using computer vision for linear feature detection

texts:
I am considering to use CV methods to detect linear features in very high resolution imagery (i've added an example below). These linear features are made by animals which transverse this section of reed vegetation.
I've seen many examples of CV being used to detect animals, objects, etc but not so much to actual map these type of features in the landscape. In classical remote sensing we would look at edge detection algorithms (e.g. sobel filters) or image segmentation but those algorithms always feel like they end up tailored to that specific case. Ideally, i could train a model in location X and have good results in location Y.
I'm still in the design phase. So, I would like to know, what algorithms you would recommend and if you can, what type of challenges you expect me to experience?


texts: Using British National Grid Transformation OSGB_1936_To_WGS_1984_NGA_7PAR in ArcGIS Desktop?

texts:
In ArcGIS 10 a new transformation has been added for converting between WGS84 and British National Grid.

What is OSGB_1936_To_WGS_1984_NGA_7PAR?

How should it be used when compared with the other available transformations?

Is it better?


texts: Multiple Markers on OpenLayers using Javascript

texts:
Using Openlayers 2.13, I am trying to set up multiple markers for New expatriates in our suburb.

As a new user of both JS and Openlayers, I am stuck - I have managed one marker but cannot add another - do I need to upgrade to Openlayers 4 or use an array?

    <title>OpenLayers Solution</title>


texts: How to load ALKIS (GML3) data in GeoTools

texts:
I wonder how I can load ALKIS GML (German gov geo data) in GeoTools.

Here is a small sample dataset from the government: http://www.schleswig-holstein.de/LVERMGEOSH/DE/Service/Download/Testdaten/alkis_bestandsdaten_etrs89zip__blob=publicationFile.zip

I put together some simple code to load it:

    GML gml = new GML(Version.GML3);
    CRSAuthorityFactory crsFac = ReferencingFactoryFinder.getCRSAuthorityFactory("EPSG", null);
    CoordinateReferenceSystem osgbCrs;

    try {
        InputStream in = new FileInputStream("/home/me/Bestandsdatenauszug_NAS_ETRS89_UTM_0348.xml");
        osgbCrs = crsFac.createCoordinateReferenceSystem("EPSG:25832");
        gml.setCoordinateReferenceSystem(osgbCrs);

        FeatureCollection fc = gml.decodeFeatureCollection(in);
        System.out.println("size: "+fc.size());
    } ...

But sadly I get the following errors:

WARNING  Error parsing: NAS-Operationen.xsd
WARNING  Could not find a schema

As I am sure the file itself is valid and well-tested, what is going wrong here? How can I successfully load this file?


texts: Are there any international water bodies codes (oceans, seas, etc.)?

texts:
As you may know, there are 2 letter international continent codes:

AF - Africa
AN - Antarctica
AS - Asia
EU - Europe
NA - North America
OC - Oceania
SA - South America

But what about water objects such as oceans, seas and so on? Are there any international codes for these types of objects (alphabetical, numeric)?

If no, are there any unique codes for them within a country?


texts: Counting building polygons in neighborhood polygons

texts: I have two layers, one of them is neighborhood data and the other is the buildings in the city. I want to count the buildings in a neighborhood (by using borders) and add this value to the neighborhood table as a field. How can I do this in ArcMap?


texts: Select feature one by one not working in pyQGIS?

texts:
I have used following code for selecting the features one by one. but when i run the script it is just selecting last feature and zoom it. I want to select each feature and zoom it.

def SelectFeatures(self):
    layer = self.iface.activeLayer();
    iter = layer.getFeatures();
    for feature in iter:
        geom = feature.geometry()
        self.iface.mapCanvas().setSelectionColor(QColor("yellow"));
        layer.setSelectedFeatures([feature.id()])
        selected_features = layer.selectedFeatures()
        self.iface.mapCanvas().zoomToSelected( layer )
        self.iface.mapCanvas().refresh()


texts:
I've one question or maybe two about Points to path in QGIS.
First, I present my case:
I have to do "travel paths" from one address to another address.
I have a database - example:

How I should arrange the data and prepare them to use the Points to path? Should I use some other tool?


texts: Copying some values from several columns to another column in attribute table?

texts:
I am trying to copy selected values from several columns to one single column in the attribute table.

Specifically, where columns (Grid_code, grid_code1, grid_code2.... etc) have values of 4 and 10, I want them to be copied to a new column, but to ignore all other values (which are all 0s).

This image illustrates roughly what I'm trying to do (in reality there are more columns and rows!)

I think I should do this in field calculator, but can't work out how- I'm inexperienced with VB and Python.

I'm using ArcGIS 10.3.1 for Desktop.


texts: Calculating areas and joining attributes between overlapping features

texts:
I am working with QGIS 3.4.8 and Python 3.6.8 on OSX El Capitan 10.11.6.

I have two shapefilea, one which is a grid of cells (each cell a feature) and the other one which is urban settlements as reported by Columbia's SEDAC project (green features in picture).

I wish to calculate the area that is urban in each cell, but I have the problem of urban settlements crossing between cell (like the yellow shown in picture) and also of cells with multiple urban settlements. I have not been able to find out a proper way to get around it.

Similarly, I have a map of Africa with country divisions and my cell grid. I wish to attribute a country to each cell, but some cells cross borders.

One way to get around it is to use centroids, which actually is an arbitrary measure and makes me lose some near-the-ocean cells, which could be negligible but I would like to make an effort to recover. I was thinking that maybe calculating the area of a country in a cell and attributing to the country with higher area the country of the cell would be somehow more rigorous. However, I have not been able to figure how can I possible do that.


texts: Simple band arithmetic on GDAL virtual rasters (VRT)

texts:
I want to calculate the difference of two bands. I have the individual bands stored somewhere, but I was thinking that maybe there's a way of producing a VRT file that is the result of this difference. I know that you can use pixel functions in GDAL, but this limits the use (C++, and only available when the function is linked in...).

It would be nice if simple "virtual" operations were available!


texts: Reading same hgt file downloaded from two sources

texts:
I have two tiles for SRTM1 for the same bounding box region, but downloaded from two sources.

I read a value for a given latitude and longitude with python.
I get always a difference of one meter for the same latitude and longitude, I changed random numbers within the boundary, the same problem remains.
Then I opened both the hgt files with binary viewer, indeed there are not the same in binary.
My question is, which one should I trust and from where comes this error.

My python code:

import numpy as np

#https://e4ftl01.cr.usgs.gov/provisional/MEaSUREs/NASADEM/Africa/hgt_merge/
file_1 = 'n34e036.hgt'   #Source of file is from the link above

#https://cloud.sdsc.edu/v1/AUTH_opentopography/Raster/SRTM_GL1/SRTM_GL1_srtm/North/North_30_60/
file_2 = 'N34E036.hgt'   #Source of file is from the link above

SAMPLES = 3601 # Change this to 1201 for SRTM3

def read_elevation_from_file(hgt_file, lon, lat):
    with open(hgt_file, 'rb') as hgt_data:
        # Each data is 16bit signed integer(i2) - big endian(>)
        elevations = np.fromfile(hgt_data, np.dtype('>i2'), SAMPLES*SAMPLES).reshape((SAMPLES, SAMPLES))
        lat_row = int(round((lat - int(lat)) * (SAMPLES - 1), 0))
        lon_row = int(round((lon - int(lon)) * (SAMPLES - 1), 0))
    return elevations[SAMPLES - 1 - lat_row, lon_row].astype(int)

#test example
lat = 34.0 + (3599.0/3600.0)
lon = 36.0 + (4.0/3600.0)

elv_1 = read_elevation_from_file(file_1,lon,lat)
elv_2 = read_elevation_from_file(file_2,lon,lat)

print(elv_1)  # this yields 423
print(elv_2)  # this yields 422

Edit: I tried also for tile N52E009.hgt, there is also difference but here for two meters.

In binary viewer:
The right side is for file_2 and the left side for file_1


texts: Correction of GPS-Track, snap to road (line, network)

texts:
Problem as shown in the image:

I like to correct the GPS track (orange) to the right road (black).
In a network of roads it is difficult to snap to the right path - for human it's obvious which one is the right - how can we compute?
Until now I tried to find a solution in QGIS:
A. in QGIS

buffer road
clip track
extract vertices
snap to road
points to path
snap to road
-> corners are often cropped

I also tried to smooth the line before. The problems remain...
B. in QGIS:

buffer road
clip track
extract vertices
snap to road
points to path
select by location (intersecting elements of road)
-> results with gaps

C. GRASS in QGIS v.edit
-> problems of snapping not to the right road perists
D. in POSTGIS:

split track in segments of 10m (roads are mostly segments of 10m)
calculate Hausdorff distance:

CREATE TABLE tracks2 AS
SELECT  a.id, a.geom
FROM roads a, track1 b
WHERE ST_HausdorffDistance(a.geom, b.geom)< 3;

-> to many wrong segments, to many missing segments
E. in POSTGIS:
I was thinking of a solution with the change of azimuth between points, but I got troubles with the points between the locations of change and corners.
(see my previous question: PostGIS: ST_Azimuth between Points in same Table)
And now I found this: Matching segments of different lengths
How can I implement this in POSTGIS? Or is there a completely different approach (Python, ...)?
I also researched the topics map matching / conflation without finding a solution that suits my problem and my possibilities.


texts: How to calculate and compare RMSE between DEM using different interpolation methods?

texts: I have data from aerial photo, it was built from stereo photo so I have 3D information that I am using for the DEM. I extracted it into shapefile. I want to get RMSE (Root Mean Square Error) from different interpolation data using ArcGIS. How would I go about solving this problem, for interpolation I plan to use TIN, Kriging and IDW, are these the correct tools? Also do I need more data?


texts: JavaScript remove Leaflet layer when checkbox is unchecked

texts:
I have an JS object where the 1st value of each key is the checkbox id, 2nd is the API where the data is and the 3rd is the color

var layers={childcares: ["#childcares",'http://127.0.0.1:8000/childcare_data/','orange'], childcare_buff: ["#childcaresbuff","http://127.0.0.1:8000/childcare_buff_data/","orange"],
            religious: ["#religious","http://127.0.0.1:8000/religious_data/","blue"], religious_buff: ["#religiousbuff","http://127.0.0.1:8000/religious_buff_data/","blue"],
            park: ["#parks","http://127.0.0.1:8000/park_data/","green"],park_buff: ["#parksbuff","http://127.0.0.1:8000/park_buff_data/","green"],
            school: ["#schools","http://127.0.0.1:8000/school_data/","red"],school_buff: ["#schoolsbuff","http://127.0.0.1:8000/school_buff_data/","red"]
            };

now I have these two functions

    function color(feature, color_chosen, opacity){
        return{
        color: color_chosen,
        fillOpacity: opacity}
       };

this color one above is fairly simple and just sets the color

the legend_click one is where I am getting problems

    function legend_click(id,layer_api,color_layer) {
        $(id).click(function () {
          var layer_add=new L.GeoJSON.AJAX(layer_api,{
            style: color(layer_add,color_layer,0.8),
            onEachFeature: function(feature,layer_add){
             layer_add.bindPopup('Buffer Type = ' +feature.properties.buff+ '
'
                            +'Address = '+feature.properties.address+ '
' +'Buffer Distance = '+feature.properties.buff_dist)
                    }
                   }
                );
            if ($(id).prop('checked')==true) {
                layer_add.addTo(map);
                //map.fitBounds(layer_add.getBounds())
            }
            else if ($(id).prop('checked')==false) {
                console.log('noury')
                map.removeLayer(layer_add);
            }
        })
    };

looping through the layers like this

for (var key in layers){
        console.log(layers[key])
        legend_click(layers[key][0],layers[key][1],layers[key][2])
    };

so in the map when I click on the checkbox the layer correctly gets added to the display but when I uncheck the checkbox is it throws NO error and prints out noury to the console but it does not remove the layer...


texts: Equal Time Maps for Transit Network (Isochrones) on QGIS

texts:
I am on a committee looking to redesign transit service in Sacramento and need to be able to make equal-time maps or "isochrones," like the attached.

A prominent difference between transit isochrones and road isochrones:

each transfer has a temporal penalty. That is, it takes time to wait for the connecting bus.
that penalty is lower where the connecting routes are more frequent and higher where the connecting routes are less frequent.
some routes only serve stations, while others are considered to have an even distribution of stops along their routes.

( I am interested in a hybrid model between frequency and coverage that still creates a high-frequency grid between major traffic generators, but with routes that take different branches between the generators. This offers both frequency and coverage. Currently, the debate is binary, frequency or coverage.)


texts:
I am trying to transform CRS to UTM using a GeoJSON file.
This is what I tried:

from sentinelhub import CRS
country_crs = CRS.UTM_32N
country = country.to_crs(crs={'init':CRS.ogc_string(country_crs)})
country.plot()
country.crs

But this code results in a completly wrong transformation followed by a warning message:

C:\Users\ASUS\Anaconda3\envs\ASV\lib\site-packages\pyproj\crs\crs.py:53: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  return _prepare_from_string(" ".join(pjargs))

trying this:

country = country.to_crs("EPSG:22391")
country.plot()
country.crs

resulted in the following: (negative values in the Y-axis):


texts: Sites which can be reached in a certain time

texts:
I'm trying to analyse the sourrounding of an archaeological site in QGIS. Therefore I've created a cost map with r.walk in GRASS (cumulative costs range from 10939500 to 14211500).
Is there anly possibility to show only areas which can be reached in 1 h? My Point of Destination is a archaeology site and I want to know which setllements in the neighborhood can be reached in a 1 h walking distance. I'm not interrested in an exact least-cost-path solution but there must be a possibility to create a simple map, or am I wrong?

The solution must be easy but I'm totally lost.

Can you help me?


texts: Changing date format for labels in ArcGis Pro with label expression?

texts:
I am trying to use an expression to convert my date field to a simple %d,%b,%Y format. Right now it is in "%Y%m%d %H:%M:%S %z"

How do I do it properly using Python to do this in ArcGIS Pro?

Q&As from ArcMap have given me errors.

My labels are coming from my [Date] field.

Here's one example of a code I tried -

def FindLabel ( [Date]  ):
  return str(datetime.datetime.strptime([Date],"%Y/%m/%d %H:%M:%S %z").date())

and I got this error:

My labels look this :

I also tried this code:

from datetime import datetime
def FindLabel ( [Date] ):
  d = datetime.strptime([Date], '%Y/%m/%d %H:%M:%S %z')
  return d.strftime('%d, %b %Y')


texts: Cropping a map into multiple pieces with QGIS

texts:
I have an .ecw map of north Portugal that has been made by joining 132 little maps at scale 1:25000 (I got the map from here: https://www.elrincondelmapa.com/Mapas_Nacionales.php). I want to crop it to get the original 132 pieces (no need to get the exact ones) and then I would like to export them as .jpg files.
I'm an absolute beginner to QGIS and digitalized maps in general. I found several explanations for extracting one region of a map with the clipper tool, but I would like to automate the task, so I could avoid repeating the same action 132 times. It would be possible to draw a grill on the map, crop it to get the 132 parts and export the images as .jpg files?
I have QGIS version 3.16.3.