texts:
I am trying to export a table from postgres to a geopackage. But I get the following errors with the following command.

ogr2ogr -f GPKG /opt/data/ways.gpkg "PG:host=postgres dbname=osm user=postgres port=5432 " -sql "select * from ways"
ERROR 1: COMMIT transaction failed: disk I/O error
ERROR 1: Terminating translation prematurely after failed translation from sql statement.
ERROR 1: Transaction not established
ERROR 1: COMMIT transaction failed: disk I/O error

If I do the same thing but export to shapefile I dont have any issues. It will successfully export. Also it may be worth noting that Postgres is running in a docker container and I am trying to run this command in a seperate container. If I run the same command to export to gpkg from my macbook it works fine.

I am wondering if it is a sqlite issue but cant find anything there that has helped either.


texts: Cannot load WFS with BBOX Parameters via Python script

texts:
I tried to load a WFS with BBox-Parameters using a python script. The script is working, but loads the whole layer, instead of just the part inside the set BBox.
The URL itself works on both the Browser and the QGis-Client and loads only the features inside the BBox.
Why isn't the BBox recognized, when using a python script?
Code:
url ="https://inspire.brandenburg.de/services/boartsubstr_wfs?version=2.0.0&srsName=EPSG:25833&service=WFS&Request=GetFeature&typename=app:bodenarten&namespace=xmlns(app%3Dhttp://www.deegree.org/app)&BBOX=367808.0001,5844987.0001,378111.0001,5858643.0001"

layer = QgsVectorLayer(url,'Example', 'WFS')

if not boden_layer.isValid():
    print ("Fail")

QgsProject.instance().addMapLayer(layer)


texts:
I am using several ESRI .NET DLLs in some custom Python scripts. For example ESRI.ArcGIS.Geodatabase.dll

On my development machine these DLLs are in the C:\Program Files (x86)\ArcGIS\DeveloperKit10.0\DotNet folder.

Now I wish to deploy the scripts to another machine. However unless the user has the ArcObjects SDK for .NET installed these DLLs are missing from their machine.

Even worse, requesting that the user install the SDK, means they also need to install Visual Studio, which is a 600MB download (for the free Express version). If they don't have Visual Studio the ESRI installer won't continue.

So should these DLLs be bundled up with the scripts (which may cause compatibility issues if service packs are added), or is there an easier deployment method?

Update:

The .NET DLLs are now installed by default in version 10 of ArcGIS. They are placed in the GAC (Global Assembly Cache). You can see them in Windows Explorer (in Windows 7) in C:\Windows\assembly (not really a folder, but you can view what is in the GAC). Looking at the assembly properties indicates the DLL should be in a folder such as C:\Windows\assembly\GAC_32\ESRI.ArcGIS.System\10.0.0.0__8fc3cc631e44ad86\ESRI.ArcGIS.System.dll but this file does not seem to exist.

Python for .NET seems to now require you use the full name when adding a reference to these DLLs. Looking at its sourcecode it looks like it previously used LoadWithPartialName so you before you could use the code below. This now returns a FileNotFound exception.

import clr
clr.AddReference("ESRI.ArcGIS.System")
from ESRI.ArcGIS.System import *

Now it seems you need to use the following:

import clr
clr.AddReference("ESRI.ArcGIS.System, Version=10.0.0.0, Culture=neutral, PublicKeyToken=8fc3cc631e44ad86")
from ESRI.ArcGIS.System import *


texts:
I have some Chinese labels to display on a map and I want them to be vertical, as the picture shows.

In another word, I want to rotate the label while keeping each character upright. How could I achieve this?


texts: Custom position of in-panel widget ArcGIS Web App builder

texts:
I need the opportunity to toggle between 'draggable' position of an in-panel widget (ArcGIS Web App Builder 2.19) and its fixed position in the footer of an application. Is there an option to do it for now?
I tried to solve it via forcibly changing a widget position by .setPosition() method and adding some css-styles (with a lot of !important rules) to a .jimu-panel class but it doesn't work correctly. What am I missing?


texts: Should I use OpenLayers + GeoNetwork + GeoServer together?

texts:
At the moment I'm doing an internship and my assignment is as following:

design (and build) a (prototype) web application
with a map on which the users can show selected data.
The server collects metadata and a direct download url from the source.
Users can download the datasets directly from the source.
A system like this that already exists and it should be able to use layers from this system through a WMS protocol.

Now there has already been research done, but it wasn't done internally and I'm not sure if she understands what GIS is. She refers to GIS as a GIS server. GIS is the whole package though, isn't it?

The research report came up with using OpenLayers, GeoServer and GeoNetwork together.
OpenLayers is referred to as the webviewer, GeoServer as the GIS server and GeoNetwork as a geoportal.
As I understand it though, GeoNetwork can display maps just fine. Should I disregard the research done and just use GeoServer + GeoNetwork? Or maybe some other combination? As long as it's free OpenSource (part of the assignment) the language doesn't matter.

I've tried (and am still researching) to find the best solution myself, but a little nudge to the correct tools would be awesome. Before this I hadn't even heard of GIS, though I do find it pretty interesting.


texts: Cross section profile distance not providing accurate distances?

texts:
I have created a DEM file, which I am using to draw a cross section along a line to see the variation of height from point A to point B along land seen below.

When it use the Profile tool on the DEM file, it produces a cross-sectional area as seen below.

I noticed from the profile that the elevation variation is correct, but my distance is out massively. A, because it doesn't give me units along the x-axis. B, because when I measure the same distance line from point A to point B on Google Maps, it gives me an accurate data which is massively out from the scale provided from using the profile tool.

Does anyone know why this is happening or how to fix it?

The steps I have taken to get to where I am now is..

Downloaded 'OSTerrain 50' ESRI shapefile from https://www.ordnancesurvey.co.uk/opendatadownload/products.html#TERR50
I converted the Contours from line.data to point.data as seen in the answer given in Creating DEM Contour files from shapefile line data in QGIS
this allowed me to create a DEM file using the SAGA Triangulation tool. from this I created a cross-sectional profile.

The line I am requiring a cross-section profile of is starting at coordinates --> 60.471758, -1.291595 and ending 60.466071, -1.272155

I am using QGIS V 3.4


texts: Creating polygon with one click in QGIS?

texts:
How can i create a polygon inside a hole made by a set of other polygons?

blue X space on the image. I want to create polygon with one click, like in autocad with option BOUNDARY.


texts:
This is a crosspost with stackowerflow and I hope for understanding and that  it will no be closed as it concerns both cartography and R.

I have already asked about plotting points dependent on region, but now my question is about overlaying points given by coordinations over gadm map. I want to show different meteostations, using data based on this site or here is  data used in code I've tried using such code:

require(ggplot2)
library(maptools)
library(rgdal)
library(RColorBrewer)
library(gpclib)
library(rgeos)
library(PBSmapping)
gpclibPermit()
rus<-load("C://RUS_adm1.RData")
proj4.str <- CRS("+init=epsg:3413 +lon_0=105")
gadm.prj <- spTransform(gadm, proj4.str)
rus<-gadm.prj

met <- read.csv2('C:\\meteo.txt')

cds <- data.frame(
                  Longitude=met$Long,
                  Latitude=met$Lat)

k<-as.matrix(cds)

popSP <- SpatialPointsDataFrame(coords=k,met["Elevation"], proj4string=proj4.str)

spplot(popSP, sp.layout=list('sp.polygons', gadm.prj))

, following advice from already mentioneed question, but it ends up with plotting points w/o gadm layer , I guess it's because of non-mentioning regions

Sorry, if question is dumb, but I will be grateful for any help


texts: Join returning some null values using ArcMap

texts:
I am a student in an ArcGIS specialization using ArcGIS 10.7.1 for Desktop, and have been struggling to complete part of an assignment.
My assignment is to interpolate ozone data from monitoring stations in California onto a raster, then join it to the map of census tracts dissolved to the approximate county boundaries, and symbolize based on the mean ozone concentration.

I carried it out by joining the ozone data to the air quality stations and exporting the joined data as a layer, then creating a TIN and subsequent raster. This appeared to be correct based on the screenshot provided of the raster in the assignment description.
I used the Zonal Statistics as Table tool to get the mean ozone concentration by census tracts and by counties and joined the tables to their respective layers; however, the tables showed some null values for census tracts and counties that should not have been null, and those areas were left blank on the map after applying symbology. I looked up the issue in Stack Exchange and it said to make sure everything was in the same projection (I used NAD_1893_California_Teale_Albers in meters), but it still turned up null values. I attached the counties layer with the selection by attributes MEAN is Null 1, as well as the ozone raster 2, which seems to be correct based on the screenshot provided in the sample.

Do you know what I may have done incorrectly and how to fix the error?


texts: How can I generate a grid that would have a quadrant like naming? Or is there a way just to rename quadrants through data driven pages. Sort of Newbish but basically creating a new firestation map book and the typical row names (as see attached), does no suit us too well. I would like to see (2nd attachement with 10) more of a quadrant. How would I generate a data driven template/grid with that? Thanks!


texts:
I have a requirement to trace path on the points present in the shapefile which is loaded in QGIS. This should enable the user to specify their own routes within the present network.
There will be existing points and lines present in the shapefile. Whenever the user clicks on the points there should be a line drawn between the points selected and the route should be saved as a shapefile.
I have searched and found in many places of options like pgRouting and shortest path routing and that is not what I require. My need is simple and due to these options being available I'm not able to find a simple solution. It's very much similar to How to route railcars along railroads?, which does not yield a clear solution.
I am new to GIS.


texts: Update Cursor to Polyline Creates Records with Missing Geometry

texts:
I am working on an arcpy script to create a polyline feature class with a single polyline feature between each 2 points in a sequence of points.  I use a search cursor to first create a table with all the necessary data, including the start and end point coordinates for each record.  Then I use an update cursor to load the data into a polyline feature class.  The output feature class has a row for every record, but is missing the geometry for most, but not all of the polylines I want to create.  I've tried rewriting it several ways, always with similar results. Any ideas on how to fix this?

rowsTbl = arcpy.da.SearchCursor(attributeTable, ['SURSEGUID', 'TRACKSEGMENT', 'START_TIME', 'END_TIME', 'START_X', 'START_Y', 'END_X', 'END_Y'])

#New polyline featureclass fields = ['SHAPE@', 'SURSEGUID', 'TRACKSEGMENT', 'START_TIME', 'END_TIME']

with arcpy.da.InsertCursor(tempOutFC, fields) as cursorTbl:
    for row in rowsTbl:
        array = arcpy.Array()
        point = arcpy.Point()
        point.X = (row[4])
        point.Y = (row[5])
        array.add(point)
        del point
        point = arcpy.Point()
        point.X = (row[6])
        point.Y = (row[7])
        array.add(point)
        del point

        #Create polyline from the array
        polyline = arcpy.Polyline(array)

        #Clear array for future use
        array.removeAll()
        #Append to the list of polyline objects

        cursorTbl.insertRow([polyline,row[0], row[1], row[2], row[3]])
        del polyline

del cursorTbl

Here's a screenshot of the attribute table I am trying to load:

As you can see, most output polylines are missing, although their records show up in the attribute table (with a shape length of zero).


texts: Zonal Statistics with very large mosaic in ArcGIS

texts:
I have the complete National Elevation Dataset from USGS in over 3000 ArcGrid files (coming to almost 190 GB), and a polygon shapefile of Colorado counties. I'm trying to use ArcMap Zonal Statistics to get the mean elevation for every county in Colorado, but after hours of processing, it repeatedly errors out. Is the mosaic just too big to work with?

I've tried Raster Calculator (to go from float to int) and Mosaic To New Raster, but both have failed.


texts: Select lines by poly using ModelBuilder?

texts:
I have a set of polys for cities and a set of lines for streets. What I'd like to do is anywhere a line segment (streets) is within a city (poly) have the value for the city applied to the line segment for example Wall St runs through City of NY but might also extend into the next city. I'd like the segments in NY to reflect NY the segments in the next city to reflect as such.

I've connected the poly and line features to ModelBuilder, created a new feature using "make feature" tool for both polys and lines. I've then fed the lines and polys into the Select Layer By Location tool. I run the model and get nothing? The model runs correctly The location is output is good (gdb I created) it's results are blank? I did refresh the gdb within ArcCatalog and still nothing.

Any ideas what I'm doing wrong.


texts:
I'm new to Mapbox GL and want to understand how i would go about the following;

I have 5 principle layers (from the same source, all vector tiles from Geoserver) which will only be drawn in isolation (ie if one is on the others are off).

Within any of those 5 layers are properties which i want to use as a filter, those properties can apply to any layer. As example below, one property is used to colour polygon (building) I would like to then do something like show only pink buildings with property x OR show any colour building where property x is present.

In GIS terms this is like running a filter within a layer and providing check boxes for the values within a property.

In total i will have 5 primary layers to which 31 different filters could be applied (filters could be multi - choice -any or only).

Something similar to this request on another platform.


texts: QGIS 3.6.2 with GRASS - raster calculator - Raster1 + Raster2 returns missing pixels

texts:
I have 5 large point datasets (shapefiles) and I'm interested in getting a heatmap on the combination of all 5 point layers.

I use the QGIS heatmap tool to create a heatmap on each layer. This is the tool:

Here's the parameters for the tool:

I run this tool for all 5 point layers to have 5 different heat maps, which have some pixels overlapping and other pixels on their own.

I.e. the rasters do not cover the same extent.

I need to add these values together (as in addition '+').

But, when I use Raster Calculator to 'add' them together, the pixels on their own disappear. I tried many different calculations, such as Raster1 + Raster2; Raster1 OR Raster2; (Raster1 + Raster1) /2  + (Raster1 + Raster1) /2 but still the same issue with data vanishing where the two layers don't overlap.

I also tried to make a large raster that covers the whole point layers area, with each pixel valued 0, but the size of the pixels were not the same! This larger raster was used to attempt to add all the other Heat Map rasters so I wouldn't be performing arithmetic on empty space. I'm not sure if the different sized raster pixels would be an issue. Making that huge raster with the same (tiny) size of the pixels in the heatmap layer is a bit impractical because it would be too big to process.

I tried nan_to_num() but still the same result.

Before (i.e. multiple Heat Maps). You can see 3 areas:

After (i.e. one combined Heat Map). You can see 2 slightly smaller areas:

How do I stop these pixels from disappearing, particularly where the layers don't overlap?


texts: Error accessing Server Status of GeoServer

texts:
I have installed GeoServer 2.17, using openJDK 11. When access the Server Status page, form aside menu on the main page after login, return a error below:
Link to full stack trace in text file.


texts: What is the unit in which value returned by OL3's geom.getArea() method?

texts:
If It is Going to return in some standard unit What is it ? so that i can convert to my desired unit.
Is it possible to set the Measuring Unit type like ESRI's API support?, That    i'm able to set in terms of Sq meter, Meter, Sq Miles, Miles, Sq Km... etc
If there is no such standard unit going to be returned, how do i can convert it to Miles, KM, Feet etc..


texts: GDAL.RasterizeLayer not ported to the C# bindings?

texts:
I am trying to use the OSGeo.GDAL.Gdal.RasterizeLayer(…) from GDAL 1.9.2 and it requires some SWIGTYPE parameters:

OSGeo.GDAL.Gdal.RasterizeLayer(
    OSGeo.GDAL.Dataset dataset,
    int bands,
    OSGeo.GDAL.SWIGTYPE_p_int band_list,
    OSGeo.OGR.Layer layer,
    OSGeo.GDAL.SWIGTYPE_p_void pfnTransformer,
    OSGeo.GDAL.SWIGTYPE_p_void pTransformArg,
    int burn_values,
    OSGeo.GDAL.SWIGTYPE_p_double burn_values_list,
    string[] options,
    OSGeo.GDAL.Gdal.GDALProgressFuncDelegate callback,
    string callback_data
);

Reading up on some old mails somewhere (I've lost the link) regarding ReadRaster and WriteRaster in the C# bindings, Tamas Szekeres states that methods requesting SWIGTYPE parameters, cannot be used (yet).

Question:
Does this mean that the RasterizeLayer method hasn’t been ported to the C# bindings or are there some clever way to utilise this method?


texts: Use ogr2ogr to filter GeoPackage by geometry type

texts: I have a GeoPackage that has features of different geometry types (Point & Polygon, see below) on the same layer. How can I filter all polygons from the file and create a new GeoPackage using ogr2ogr?


texts: Why sectors drawn by Shapetool become distorted and can it be avoided?

texts:
I try to draw sectors with 60...90 degrees azimuths by Shapetool (see picture 1). Every sector seems to get a bit different from each other. I guess this has something to do that sectors get distorted further away from central meridian of my projected coordinate system (RT90 in Sweden).

In second picture each sector should be from 60...90 azimuth so that southernmost side of sector should be horizontal. But clearly that is not the case. Sectors to the right have southernmost side pointing upwards while sectors to the left have them pointing downwards.

Is there a way to get "undistorted" sectors and be shown somehow in my projected coordinate system (RT90)?


texts:
I am working on a city-planning research project using QGIS. For a given area of interest, I have downloaded an open street map layer (the 'highway' keys) with all the roads in that area.
What I am trying to do is:

Determine how many kilometers of road in the OSM layer are within +/- x degrees of true north (ie. for all the road segments in the map, which of them are within some tolerance of true north).

Create a new layer containing just the road segments that were found within (1) above.

My preference is to do this programmatically in python.
I am relatively new to QGIS, what would be the best way to approach this problem?


texts: Insert large list of fields and values into feature class

texts:
I want to create a toolbox
I have a list of point coordinates, that define a polygon and a list of attributes, that each polygon has.
points_list = [(323917.71, 479699.45), (323923.36, 479699.74), (323923.89, 479709.52), (323918.24, 479709.23)],
attributes = [('ADDRESSID', '5'), ('BUILDINGID', '4'), ('LANDID', '35'), ... ]

I want to insert the information into a feature class. I also want to define the projection, that would be selected by the user, inside the Toolbox.

I have tried:

list_of_fields = [key_value[0] for key_value in attributes ]
insertCursor = arcpy.da.InsertCursor(fc,
                                     ['SHAPE@', list_of_fields], spatial_ref)
polyArray = arcpy.Array()
for pointsPair in points_list:
    newPoint = arcpy.Point(*pointsPair)
    polyArray.add(newPoint)
newPoly = arcpy.Polygon(polyArray)
list_of_values = [key_value[1] for key_value in attributes ]
insertData = newPoly, object_id
insertCursor.insertRow(insertData)

Nothing is recorded in the feature class.


texts: OpenLayers 3 Style - text label with size not depending on scale

texts:
I try to label polygons from a WFS-Source with an attribute. This works fine, but my problem is that the size of this label is scale depending. It has the same size on every scale. But i want it to become smaller/bigger if i zoom out/in the map. Like the polygons become smaller zooming out the map.
How can i get this with Openlayers?

This is the problem in pictures:

Here is the code i use:

var Flst_GbR = new ol.layer.Vector({
        source: new ol.source.Vector({
        loader: function(extent) {
            $.ajax('http://192.168.10.9:8096/geoserver/IIS_Becker/wfs', {
                type: 'GET',
                data: {
                    service: 'WFS',
                    version: '2.0.0',
                    request: 'GetFeature',
                    typename: 'IIS_Becker:Flurstücke_GbR',
                    srsname: 'EPSG:2398',
                    outputFormat: 'application/json'
                }
            }).done(function(response) {
                Flst_GbR.getSource().addFeatures(new ol.format.GeoJSON().readFeatures(response));
            });
        },
        strategy: ol.loadingstrategy.all,
        projection: customprojection
        }),
        style: function(feature, resolution) {
            var stylestößen = new ol.style.Style({
                stroke: new ol.style.Stroke({
                    color: 'rgba(255,0,0,1.0)',
                    width: 1
                }),
                fill: new ol.style.Fill({
                    color: 'rgba(253,184,184,0.5)'
                }),
                text: new ol.style.Text({
                    text: feature.get('flurstueck')
                })
            });
            var stylegörschen = new ol.style.Style({
                stroke: new ol.style.Stroke({
                    color: 'rgba(255,128,0,1.0)',
                    width: 1
                }),
                fill: new ol.style.Fill({
                    color: 'rgba(255,193,132,0.5)'
                })
            });
            if ( feature.get('gemarkung') == 'Stößen') {
                return [stylestößen];
            } else {
                return [stylegörschen];
            }
        }
    });

How do i have to change the code on the ol.style.Text-object to get it work?


texts: Shapefile from heading and distance directions for property boundary

texts:
I have a survey with heading and direction points describing the boundary for a rural property in Brazil. I am in correspondence with someone there who needs a shapefile (I believe he mentioned kml or kmz format file as well) of this boundary.

The format of the information is as follows:

Heading                 Distance            Alignment

NW 79º 41' 27"         170.08 meters         1-2
NW 58º 27' 42"         119.83 meters         2-3
...

The alignments I assume references points on some old maps I have of the property, beginning with a starting point described by some known area which has been marked in some way. Neither this map nor the survey has any lat/lon coordinates.

I know Python but do not have experience with GIS.

What terms can I search for, or what documentation would help me understand how to create the polygon he needs from the kind of information I have?

Is shapely or QGIS the tool I should familiarize myself with?

Is there a process for converting the kind of information into a polygon and what is the name of this process so I can familiarize myself with it?

There seems to be no way I could think of to ascertain the starting point considering it is described by surrounding landmarks but perhaps there is a way to create the general polygon and have a third party figure out how to orient it so to accurately overlay it on a map. The only other way is to contact the person living on the property and have them somehow use their smartphone to record the lat/lon of the starting point (maybe have them take a picture to send to me and I extract the coords from the exif data?).

Would this be helpful?

I attempted to create a polygon by accumulating points as suggested by @Spacedman (after converting the the bearing of deviation from N/S to degrees the python function I used for conversion are below) and the shape did not resemble the map boundary at all. As suggested by @MartinF perhaps I could invert some of the coords to try and see if I can make the shape more consistent with the map, but then I found a more recent map(2005).

Except now in this new map there seems to be some distances mixed in under the bearings column, I can't really make heads or tails of it.

(section)  (heading)       (distance)
Trecho     Rumo            Distancia (m)
1-2      R=1600.94m      D=138.36
2-3      56º02'14" SW    635.47
3-4      R=2900.00m      D=298.80
4-5      50º08'02" SW    187.45
5-6      R=10.00m        D=13.67
.....
17-18    38º28'17" NE    415.16
18-19    38º52'52" NE    305.91
19-20    40º11'22" NE    182.30
20-21    34º37'52" NW    122.53
...

At least I now know it is SAD-69 format and it contains zone and hemisphere info (although hemisphere is obvious) along with other UTM info.

Python functions to transform bearings to degrees:

def dms2dd(degrees, minutes, seconds):
    dd = degrees + (minutes/60) + (seconds/(60*60))
    return dd

def degree_from_bearing(bearing):
    # first transform from `NW 50º22'12"` to individual elements
    ns, we, deg, mins, secs = parse_dms(bearing)
    degree = dms2dd(deg, mins, secs)
    if ns == 'N':
        start = 0
        if we == 'W':
            start = 360
            degree *= -1
    elif ns == 'S':
        start = 180
        if we == 'E':
            degree *= -1
    else:
        raise ValueError(f'Heading must be North (N) or South (S), got "{ns}" instead.')
    result = start + degree
    assert result >= 0, f'Resulting degree is negative ({result}).'
    return result

def transform_data(data):
    # implementation of @Spacedman's algorithm
    x = y = 0
    coords = [(x, y)]
    for (heading, distance) in data:
        delta_x = distance * Decimal(math.cos(heading))
        delta_y = distance * Decimal(math.sin(heading))
        x += delta_x
        y += delta_y
        coords.append((x, y))
    return coords


texts: Creating a PostGIS database on an external hard drive?

texts:
Is it possible to create a PostGIS database on an external hard drive? If so, is there anything to be mindful of and what would be the host address to connect to it with pgAdmin and QGIS ?

Background: I am part of a UK wide team, all remote workers. We have a PostGIS db on a cloud server for our shared project work which with the offline editing plugin works fine. We all use locally stored OS Mastermap. We could all have a local PostGIS db and use OS Translator II to populate it with the OSMM we want and styled. But all our machines are 'locked down', so only the IT department can install software or change files on the C drive. So as an option I thought we could each use an external HDD for our OSMM.


texts: Mapbox GL JS style vector tile interaction in OpenLayers?

texts:
I am currently researching different approaches for integrating vector tiles (especially MVT and/or GeoJSON tiles, served from GeoServer) with OpenLayers (^4.0.0) in a manner that would retain access to original vector feature attribute data and interactions therewith. Since there's potentially a lot of vector data to interact with, vector tiles would be highly preferred over raw GeoJSON due to rendering performance (of course TopoJSON might be an alternative to all this).

Interaction with standard OpenLayers API MVT integration is possible pixelwise with forEachFeatureAtPixel. However, due to the nature of vector tiles, it seems not possible out-of-the-box to access the "original" feature data behind e.g. MVT. with OL methods.

Following approaches have been suggested for accessing data:

1) Using Mapbox's vector-tile library
https://github.com/mapbox/vector-tile-js

2) Loading features in GeoJSON and converting client-side to vector tiles for rendering with:
https://github.com/mapbox/geojson-vt

At least to my knowledge, OpenLayers already incorporates some parts of vector-tile-js, but it's stayed unclear to me to which extent and in which methods.

I've compared OpenLayers with Mapbox GL JS's, the latter of which seems to have figured out interaction with rendered vector tiles and original data, for instance in this example:
https://www.mapbox.com/mapbox-gl-js/example/query-similar-features/

Digging deeper into the Mapbox GL JS API, it seems that Mapbox uses two different data query methods, map.querySourceFeatures and map.queryRenderedFeatures. Looking into Mapbox GL JS source, this seems to have been somehow managed via maintaining both a cached original data and vector tiles. https://github.com/mapbox/mapbox-gl-js/blob/master/src/source/query_features.js

My real question here is, is there a similar (built-in) approach in OpenLayers, that would allow me to access both "raw" and "rendered" vector tile data, or perhaps using geojson-vt or vector-tile-js to switch between the two? For instance in the manner that would allow one to use methods akin OL's Vector method "getFeaturesInExtent"? i.e. figure out which vector tiles are in View, query corresponding "raw" data and pass on (all) accompanying feature attribute data.


texts: How to host the results of qgis2leaf on google drive?

texts:
I am wondering how to host the results of the plugin qgis2leaf on google drive to create a free hosted webmap?

Here is an example:
https://39b425ec138a1e68ecbac84487470661e6b8b943-www.googledrive.com/host/0B5ETLwh9h00iX1YyU01tNk1wU0U/

They said they created that map with qgis2leaf and it is hosted on google drive. I can get my data on google drive and the box for the map shows up but the contents of the map are not there.


texts:
I am trying to find a fast way to merge two shapefiles and i ran over pyshp which apparently could do the job. However looks like this is not working for me because i receive the following error

file "c:\net\shapefile.py" ine 443 in __dbfHeader
fieldDesc = list(unpack("<11sc4xBB14x", dbf.read(32))) error: unpack requires a string argument of length 32.

I have no idea what could be wrong and i've been trying to find a solution but no luck

My code looks like this:

import glob
import shapefile
import sys
sys.path.append('C:\NET')
files = glob.glob("C:\TestSHAPE\merge\*.shp")
w = shapefile.Writer()
for f in files:
    r = shapefile.Reader(f)
w._shapes.extend(r.shapes())
w.records.extend(r.records())
w.fields = list(r.fields)
w.save("merged")


texts: Where did checkbox option in Points to Path tool go?

texts:
I made considerable use of the "Points to Paths" tool. I forgot which QGIS version I was using, but it was most likely 2.x, probably 2.14.
According to my memory and notes the tool had a checkbox option titled 'Lines per Vertex'.  However, I don't see that option as of 2.18. or 3.8.2.
However, the according to the tool's github page, no changes have been made to the tool in the last 4 years.
The checkbox option was very useful for connecting points as their values changed in sequence, and I'd like to take advantage of its function again.
Accordingly, I'm confused. I know for sure that I used that checkbox, but there's no evidence that it existed. Am I dreaming?


texts:
I want to build a site using Django and have a user authentication system. In that site I want to serve maps and other data, but not to all users.

What I am interested in is: can I give the maps some viewing permissions? I could build maps with Geoserver and OpenLayers, I just don't know how to serve them for each user.

Expected behavior: Have some maps (WMS, GeoJson features etc.) and a user with some permissions. Logged in he should see a few of those maps (as for his permissions). If he logs out he shouldn't see the map (by visiting a url or something similar). If his permissions are revoked he can no longer see maps.

Can this be done using Django and Geoserver (or other map service solution)? What should I start learning?


texts: How to retrieve transport infrastructure data from google map?

texts: I imported into postgis osm data using osm2pgsql. My objective was to get a database with roads elements for a specific region of Italy (i.e. traffic lights, driving directions and turning direction at junction). However there are a lot of missing data, because almost all the cells are empty. I would like to know if I did something wrong and if there is an alternative way to obtain these data. For instance I would like to know if there exists a way to retrieve data from google map (transport infrastructure or even  real time traffic information)?


texts: How To Create Map View without XAML in ArcGIS Runtime for .NET?

texts: I have three projects here. There are Client, Server, and UserControl which contains with ArcGIS. Let me call the UserControl as UserControl map. I use the UserControl in the client. To synchronize between client and server, I created libraries which contain functions of calculating map, such as the function of GetMapPosition, GetMapDepth, etc. In the GetMapPosition, I need to use the function of ScreenToLocation. For access it, I must create the map view in my libraries class. That's why I must create the map view without XAML. Anybody have a solution for this?


texts: I'm developing an application, and I'm using arcgis javascript api, I'm trying to print the map using Print() object, but when pressing on the print button, it syas printing, and after amount of time, it says time out, and didn't print anything, this happens when I added the dynamicLayers to the map, what's wrong?


texts: Where to download shapefiles of interest for Honduras

texts:
Where can I download free Shapefile data for Honduras including:

roads
departments
municipalities
points of interest


texts: QGIS 3 wrong conversion from QGIS layer to KML or GeoJSON file

texts:
I want to save my work as the .kml file
Unfortunately, I've encountered a serious problem. Once loaded in Google Earth I have something like below or nothing.
A similar thing occurs with the geoJSON file. Does anyone know where the problem might be? I've never seen it before
Both the project and work are in EPSG:4376 WGS84 coordinate system


texts: UnicodeEncodeError in QGIS when setting a path

texts:
When I'm trying to batch reproject a bunch of files using GDAL warp (reproject) from processing module in QGIS (2.14) I'm getting UnicodeEncodeError:

UnicodeEncodeError: 'ascii' codec can't encode characters in position 4-8: ordinal not in range(128)
Traceback (most recent call last):
  File "C:/PROGRA~1/QGISES~1/apps/qgis/./python/plugins\processing\gui\BatchOutputSelectionPanel.py", line 83, in showSelectionDialog
    print filename, selectedFileFilter
UnicodeEncodeError: 'ascii' codec can't encode characters in position 4-8: ordinal not in range(128)

This is a very common error when working with Russian/Unicode paths (and/or Russian locale of QGIS) and is not specific to this particular module but rather to using any Unicode paths within a workflow.

Is there any way to overcome this error?


texts: Formula for transforming from EPSG:4326 to EPSG:3003 or EPSG:3004

texts:
This system is used in Italy.
I came across this:

It says that lambda zero is 9° if it's to epsg:3003 otherwise if it's to epsg:3004 it's 15°

I came up with this Java code, but the result is wrong:

public static Posizione fromGeo2MonteMario(Posizione from, String srid)
{
    double fi = from.getLatitude();
    double la = from.getLongitude();
    double b = la - (srid.equals("EPSG:3003") ? 9d : 15d);
    fi = convert2Rad(fi);
    la = convert2Rad(la);
    b = convert2Rad(b);
    double w = Math.sqrt(1 + (0.0067681702 * Math.pow(Math.cos(fi),2)) );
    double a = Math.atan( Math.tan(fi) / Math.cos(w * b));
    double v = Math.sqrt(1 + (0.0067681702 * Math.pow(Math.cos(a),2)) );
    double y = 6397376.633 * asinh( Math.cos(a) * Math.tan(b) / v);

    Posizione to = new Posizione();
    to.setNorth( (111092.08210 * a) - (16100.59187 * Math.sin(2*a)) + (16.96942 * Math.sin(4*a)) - (0.02226 * Math.sin(6*a)) );
    to.setEast( y + 500000);

    return to;
}

public static double asinh(double x)
{
    return Math.log(x + Math.sqrt(x*x + 1.0));
}

public static double convert2Rad(double degAngle)
{
    return degAngle * Math.PI / 180;
}

For example, latitude 42 and longitude 11, is EPSG:3003 so it should be around: East: 1665646.4 North: 4651793.5
You can test it here (or check the example in the pdf)


texts: GeoJSON from react-leaflet doesn't render the layer over the map

texts:
I have a GeoServer running on my localhost also I have a react app which consumes a WFS service to get a GeoJSON data to render as a layer over the map.

These are my dependencies on my package.json

"axios": "^0.19.0",
"leaflet": "^1.5.1",
"react": "^16.8.6",
"react-dom": "^16.8.6",
"react-leaflet": "^2.4.0",
"react-scripts": "3.0.1"

And this is my HomePage.js in my react app

// @packages
import React, { PureComponent } from 'react';
import axios from 'axios';

import {
    GeoJSON,
    Map,
    TileLayer
} from 'react-leaflet';

// @constants
const GEOSERVER = 'http://localhost:8080/geoserver/wfs';

const REQUEST_PARAMS = {
    outputFormat: 'application/json',
    maxFeatures: 250,
    request: 'GetFeature',
    service: 'WFS',
    typeName: 'gaia:mining_titles',
    version: '1.0.0'
};

class HomePage extends PureComponent {
    constructor() {
        super();
        this.map = null;
        this.state = {
            center: [6.217679, -75.570648],
            data: null,
            zoom: 8
        };
        this.handleOnMapMounted = this.handleOnMapMounted.bind(this);
    }

    componentDidMount() {
        axios.get(GEOSERVER, { params: REQUEST_PARAMS })
            .then(({ data }) => this.setState({ data }))
            .catch(error => Promise.reject(error));
    }

    handleOnMapMounted(evt) {
        this.map = evt ? evt.leafletElement : null;
    }

    render() {
        const {
            center,
            data,
            zoom
        } = this.state;

        return(
            <Map
                center={center}
                id="home-page-map"
                ref={this.handleOnMapMounted}
                zoom={zoom}
            >
                <TileLayer
                    url="https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png"
                    attribution="&copy; OpenStreetMap contributors"
                />
                <GeoJSON data={data} />
            </Map>
        )
    }
}

export default HomePage;

When the componentDidMount is called the data is set on the state but doesn't render anything and doesn't show errors it just nothing happens.

I already enabled the CORS configuration on my web.xml from my GeoServer

EDIT:

This is the response from the GeoServer

{"type":"FeatureCollection","features":[{"type":"Feature","id":"mining_titles.4133","geometry":{"type":"MultiPolygon","coordinates":[[[[689969.6239,566781.7359],[689969.6246,566381.7341],[690169.6234,566381.7345],[690169.6237,566231.7353],[689873.623,566283.7334],[689757.6244,566421.7352],[689757.6241,566621.734],[689835.6237,566781.7357],[689969.6239,566781.7359]]]]},"geometry_name":"geom","properties":{"id":4133,"objectid":1,"oid_":9918,"id_1":4133,"codigo_exp":"HCS-152","codigo_rmn":"HCS-152","codigo_ant":"HCS-152","fecha_cont":"2006-06-29Z","fecha_insc":"2007-05-17Z","grupo_trab":"PAR PASTO","estado_exp":"TITULO VIGENTE-EN EJECUCION","modalidade":"CONTRATO DE CONCESION (L 685)","minerales":"MATERIALES DE CONSTRUCCION","titulares":"(69021277) CLAUDIA CECILIA BETANCOURT YEPES\\ (12972275) GERMAN NAZARENO BETANCOURT PANTOJA","municipios":"ORITO-PUTUMAYO","fecha_term":"2037-05-16Z","shape_leng":1742.81141444,"shape_area":118628.3028,"grupo":null}}],"totalFeatures":8861,"numberMatched":8861,"numberReturned":1,"timeStamp":"2019-07-16T16:39:34.595Z","crs":{"type":"name","properties":{"name":"urn:ogc:def:crs:EPSG::3116"}}}


texts: Is it possible to set the layer which the trace tool traces?

texts:
Let's say I have at least 3 Layers. I want to reshape a polygon of layer 1, tracing a polygon from layer 2.

This is no problem if you have only those 2 layers.
But now I need layer 3 (or 3 to 10) turned on for information when tracing layer 2. What can I do that the trace tool only traces layer 2?
I know what to do when I create a new polygon (then I just select the polygon of layer 2 and set the trace tool to "trace selected features only"). But when I want to rehape the polygon of layer 1 I don't know what I can do.

I want trace only on the red lines.
Does anyone know an answer?


texts: Python for loop doesn't work for Spatial ETL Tool

texts:
I've created a Spatial ETL tool from ArcGIS 10.1 and would like to run this tool to all 200 featureclasses in a geodatabase.  Therefore, I've written a python script to do this automation. (btw, i'm new to arcpy) However, I get an error message stating that my featureclasses do not exist in the geodabase.  It does exist though.

"Failed to execute. Parameters are not valid.
'Export1' does not exist.
Failed to execute (SpatialETLTool152)"

I couldn't find any solution on the internet.

`import arcpy, os
from arcpy import env

arcpy.ImportToolbox("M:/GIS-Data/Poin2Line.tbx")

env.workspace = r"L:\\ExportLines\\Test.gdb"
FakeTide_All1_tif = r"L:\\ExportLines\\FakeTide.tif"
Test_gdb = r"L:\\ExportLines\\Test.gdb"
Destination_Comma_Separated_Value__CSV__Directory_ = r"L:\\ExportLines\\Output"
fClass = arcpy.ListFeatureClasses()

for fc in fClass:
    arcpy.gp.toolbox = "M:/00-GIS-Data-4-All/Poin2Line.tbx";
#arguments (Tiff, output, input)
    arcpy.gp.SpatialETLTool152("'L:\\ExportLines\\FakeTide.tif'", Destination_Comma_Separated_Value__CSV__Directory_, fc )`


texts: Graduated symbols not appearing for all polygons

texts:
I'm using ArcGIS 10.2.2.

Trying to map some city populations using graduated symbols. Should be simple, but while most of the polygons are being assigned symbols, some aren't.

For example, in the attached picture I have a city highlighted in blue - North Bay. A symbol isn't being assigned to this city at all. I'm displaying the POP_TOT_19 field, and obviously there's data there for North Bay looking at the attribute table. It should be assigning the largest symbol, yet shows nothing. What am I doing wrong?

There are no joins. All entries were unselected at the time of classification.

I think I have a very round-about way around it - but not a fix, and I still don't know why it's doing this. Best I can think to do is export the few cities individually and manually change the symbology to the appropriate symbol before grouping to the layer with the rest of the data. Anyway, thoughts?


texts: Can shp2sde on ArcSDE 9.3 server operate with SDE 10.3 database?

texts:
We have 2 servers running ArcSDE 9.3, and one Oracle database at SDE 9.3. Only one of the ArcSDE 9.3 servers is to be upgraded to ArcSDE 10.3, and the Oracle SDE database will be migrated from SDE 9.3 to SDE 10.3. The hope is that the sdp2sde commands on the non-upgraded ArcSDE 9.3 server will be able to operate on the upgraded SDE 10.3 database. I firmly believe that the upgrade to SDE 10.3 will update the database schema such that shp2sde on ArcSDE 9.3 will no longer work. Am I right?

I don’t believe we will be able to download a version of shp2sde for SDE 10.3. It was in was in ArcSDE Administration Command Line Utilities, aka. SDE Commands.

From what I can tell, if we already had 10.2.2 we would perhaps have some SDE commands which would work against SDE 10.3, but not shp2sde.

ESRI expect us to use Python rather than shp2sde.

This is the closest related information that I can find:
Can ArcSDE command line tools be used against 10.3 and newer geodatabases?

Confirmation that shp2sde is not in 10.3: shrug-l: SDE Commands on 64 Bit Servers
The commands not supported on 64-bit systems are:
cov2sde, sde2cov, sde2shp, shp2sde, shpinfo, sde2tbl, tbl2sde


texts: Counting unique occurrences during Spatial Join?

texts:
I have a point shapefile (points.shp) with fields [owner], [name], etc.. [owner] and [name] are both of type string. The owner contains a single string, the unique owner. Name contains a list of names, seperated by ";".

[owner]     [name]
   a     ";wood;grass;"
   a     ";grass;house;"
   a     "house"
   b     ";wood;tree;"
   b     ";grass;house;"
   c     ";tree;"

Now, what I am currently doing is aggregating all points of each unique name-string in [name] by a distance and Spatial Joining those back to count the amount of points containing this string within the Aggregated Polygon. The following script is running in Model Builder.

Select all [name] LIKE '%;wood;%'
Aggregate Points by distance z = agg_wood.shp
Spatial Join agg_wood.shp with points.shp (selection still alive)
= adds a field "Join_Count" to agg_wood.shp
Select [name] LIKE '%;house;%' etc.

The aggregated Polygon Shapes (about 7000) are merged into a final Shapefile. So, my final, merged shapefile will have all counts of [name] occurrences:

shape_wood.shp:
[FID]  [Join_Count]
  1         10
  2         11
  3         2
  4         1
  5         13
  6         4

Finally, I merge all aggregated polygon shapefiles, which works. The point file has something about 800,000 points, about 5,000 unique name strings, and about 20,000 unique owners.

My problem is I would like to add a field to agg_wood.shp during the process that sums the number of unique occurrences for the field [owner] for each aggregated polygon. Lets say there are 5,000 points in an area and 500 of those contain "wood" in the list of strings from the field [name]. These are aggregated into one or more polygons. However, there may be 1 up to 500 different owners who generated "wood". I want to add a field [occ_own] with the number of unique owner-strings who added "wood" in the field [name]. So my intermediate output file (agg_wood.shp) would look like this:

shape_wood.shp:
[FID]  [Join_Count]  [Occ_Own]
  1         10        1
  2         11        1
  3         2         2
  4         1         1
  5         13        13
  6         4         1

shape_house.shp:
[FID]  [Join_Count]  [Occ_Own]
  1         23        3
  2         10        5
  3         3         3
  4         1         1
  5         150       1
  6         2         1

Then I would merge all of those to a single shapefile of all [name] polygons. But I can't imagine how to modify my process so it can calculate the number of unique owners for each polygon area and [name].

Does anyone have an idea?


texts: Using Geoserver to serve Tiles in zoom/x/y.png format?

texts:
I have an application that receives tiles from an OSM tile server using the url format URL/zoom/x/y.png . I'm looking into using GeoServer instead of the OSM tile server. I understand that GeoServer uses the WMS standard.

Is it possible to serve tiles using GeoServer with the endpoint url format similar to the OSM tile server?


texts: Subset a SpatialLinesDataFrame in R based on the name field containing a certain string

texts:
I have a SpatialLinesDataFrame in R that is a road network dataset. The names of certain roads begin with the string "NH" to indicate that these are national highways. E.g., NH44, NH75 etc. Roads that are not national highways have names beginning with other strings.

I want to subset the roads dataset to include just roads that are national highways. Is there any way to do this in R?


texts:
I'm trying to create a vector grid for a municipality. Do you know why does it crash? The scale I'm using is 1:1000000.

Any tips?


texts: Automatically chopping lumps off polygons

texts:
I have a load of polygons that were originally digitised from CAD. These unfortunately include traingles that were used to sybolise the direction the polygon is facing. As can be seen here:

The trangles are of uniform shape and size and are too numerous to remove manually, so I am trying to automate a sytem that finds them (based of angle and size?) and clips them off the polygons leaving the unsymbolised polygons. I was planing on either doing this using model builder and FME, but I'm no sure which tools/transformations to use. Has anyone got any experience of doing this, and if so how did they go about it and what tools did they use?

Thanks


texts: Create Shapefile with auto increment primary key in QGIS

texts: Is there a way to create an auto increment primary key column when creating a new Shapefile in QGIS?


texts: Error 000732 - Input Layer "does not exist or is not supported" when using 'ApplySymbologyFromLayer_management'

texts:
Using Arc 10.2

I am getting the following error message...

  ExecuteError: Failed to execute. Parameters are not valid.
      ERROR 000732: Input Layer: Dataset Q:\Data\Infrastructure\Infrastructure.gdb\LiftStations_Active does not
  exist or is not supported
      Failed to execute (ApplySymbologyFromLayer).

...when I run a script containing the following code snippet:

dest = arcpy.GetParameterAsText(1)
destWorkspace, destName = os.path.split(dest)
arcpy.FeatureClassToFeatureClass_conversion(source, destWorkspace, destName, "", mapping)
in_symbology_layer = r"Q:\Data\lyrfiles\LiftStations.lyr"
arcpy.ApplySymbologyFromLayer_management (dest, in_symbology_layer)

There are no problems with the first four lines of my snippet (all the parameters are legit and I'd rather spare you the entire or extraneous sections of my code).  The error message says nothing about those lines and I've confirmed the "in_symbology_layer" exists (Q:\Data\lyrfiles\LiftStations.lyr).

I have double checked - the FC mentioned in the above error message does exist.  In fact, since it comes directly from the tool's 2nd GetParameterAsText, the tool warns me that it "already exists" with the exclamation symbol beside the parameter field.

Why am I getting this error?

Further Thoughts

Am I not allowed to use "ApplySymbologyFromLayer_management" on a Feature Class?  On a newly-created FC?

Does it have anything to do with the fact that since "dest" is an Output Parameter it gets deleted as soon as the script begins?  (reference: "Why is “arcpy.GetParameterAsText” deleting Target Feature Class?")

Does it have anything to do with the fact that "dest" then gets recreated only two lines prior to 'ApplySymbologyFromLayer_management'?


texts: Getting uploaded shapefile into geometry field on PostGIS table?

texts:
I am working with a PHP web application where we have a user submitting a ShapeFile.

I am able to convert that ShapeFile then to a format of geojson.

What I am wondering now is how I can best get that geojson into a PostGIS table as a geometry column?


texts: Selecting features and saving selection as shapefile using PyQGIS?

texts:
For some reason my code is not doing what I expect it to do.

from qgis.core import *
from qgis.utils import iface
import processing

# Load the layer
layer = iface.mapCanvas().currentLayer()

bluespaces = QgsExpression ( " \"ITEM2012\" = 'Wetlands' OR \"ITEM2012\" = 'Water'")
selection = layer.getFeatures(QgsFeatureRequest(bluespaces))
ids = ([k.id() for k in selection])
layer.setSelectedFeatures( ids )
newlayer = QgsVectorFileWriter.writeAsVectorFormat(layer,
               'U:\PythonTest\bluespaces.shp', 'System',
               QgsCoordinateReferenceSystem(25382), 'ESRI Shapefile', bool(True))

This should give me a shapefile containing the blue spaces from the input layer. It works fine with other selection criteria and I checked if something was wrong with the selection criteria used above numerous times, but to no avail.

Can someone tell me what I am doing wrong?


texts: Deleting features with identical table attributes in QGIS

texts:
How do I delete features with identical table attributes in QGIS? Otherwise, is there a tool like the 'Delete Identical Tool' in ArcGIS?
For example, if I had 4 rows, and 3 columns, where I wanted to delete identical features based on 2 columns:

id
value 1
value 2

1
A
m

2
A
m

3
A
p

4
B
m

I want it to be reduced to only 1 feature for each identical set of attributes:

id
value 1
value 2

1
A
m

3
A
p

4
B
m

In ArcGIS Desktop, this is a simple tool shown below called 'Delete Identical', where I can select one or more columns at the same time. I use this tool all the time, am switching to QGIS.


texts: determining what polygons intersect another layer

texts:
I am attempting to establish what spatial relationship exists (if any) between two layers.

I have two layers - one is a road network with a buffer of 500m, and the other is a collection of polygons relating to environmental features.

The road network buffer layer comprises of multiple roads and I would like to know what envrironmental features intersect what roads. Or more specifically I would like a column in the environmental layer attribute table which describes what road it interacts with?

Apart from using the 'select by location' spatial tool and then manually populating an empty column in the environmental layer, is there another spatial tool or method that I haven't considered?

The data itself is being stored in PostGIS and used in QGIS, but it would be possible to use ArcGIS if necessary.

Thanks for your time.


texts: Features lost when creating "Local copy for editing"

texts:
I want to edit a hosted feature service from AGOL (no server). I followed the steps to create a local copy for editing, however when I do this - not all of the features exist in the local copy. I cannot check to see if there are any errors or messages since it was not a geoprocess.  Does anyone know why this could be happening and how to resolve this issue? It happened multiple times.

STR:

Sign into AGOL from ArcMap
Add hosted feature service to map (drag and drop from Catalog)
Navigate to service in TOC (Listed by drawing order)
Right click, editing, create local copy for editing
Open attribute table for one layer... Holds 20 features, when it should have 68


texts: How to Make a Thematic Map with ORACLE using MapXtreme?

texts:
I want to make a thematic map with Oracle source longitude and latitude values.
I  have made thematic maps with tab file source.  Part of the code looks like this:

MapInfo.Engine.Session.Current.MapFactory[0].Layers.Insert(0, new MapInfo.Mapping.FeatureLayer(MapInfo.Engine.Session.Current.Catalog.GetTable("SPA_POINTS"), "TempThematicLayer", "TempThematicLayer"));
 MapInfo.Mapping.FeatureLayer lyr = harita.Layers["TempThematicLayer"] as MapInfo.Mapping.FeatureLayer;

    MapInfo.Mapping.Thematics.RangedTheme ranthm = new MapInfo.Mapping.Thematics.RangedTheme(lyr, "MAGNITUDE", "SPA_POINTSThematic", 10, MapInfo.Mapping.Thematics.DistributionMethod.EqualCountPerRange);

What can I use  instead of:

GetTable("SPA_POINTS")

or how can I make a thematic map with an Oracle table?


texts: Database compatibility between QGIS and ArcMap

texts:
I'm picking up ArcMap having started off in QGIS. I've been struggling for a while now to make databases that both QGIS 3.10 and ArcMap 10.7 can read and write. It seems ArcMap imposes much stricter conventions on data. I've run into the same problem with geopackages and now PostgreSQL. This isn't an exhaustive list of the headaches but as far as I can tell ArcGIS:

Only accepts lowercase letters, numbers and underscores in table names;
Will throw a SQL syntax error on loading if a table name starts with a number;
Only accepts lowercase letters, numbers and underscores in variable names - it'll show the geometry of layers that break this rule but throws an error when you try to view the attributes; and
Won't edit or display attributes for layers that include invalid geometries.

While I haven't tried naming a table or variable #, @ or " in QGIS, it seems like anything goes. How can I open these databases in ArcMap?


texts:
We're running our own CartoDB stack and all is working well, but in writing an API testing script, that creates and subsequently deletes a dataset, we've discovered we're not sure how to programmatically delete a dataset completely.

Existing questions I've found on the net about this all lead to using the SQL API to simply DROP the dataset's table, but that leaves references to that table in place in the primary Carto database, and that means the dataset is still listed in the UI but errors out when trying to view it -- basically the Carto UI application still thinks the dataset exists, even without its table.

That primary database is also outside the range of the SQL API so a DELETE query isn't plausible.

Has anyone found a more "correct" way of doing this? It seems like an incredibly huge gaping hole in API functionality not to be able to delete something you've created.

Addendum:
I see that the SQL API has the "cdb_cartodbfytable" function that registers a table (created with the SQL API) so that it's visible to the editor, but doesn't seem to have an equivalent function to deregister upon dropping.

Edit: After a few days of no response I've also added it as an issue in the software's repo. https://github.com/CartoDB/CartoDB-SQL-API/issues/401


texts: Efficient ways to read table data from a large feature layer?

texts:
I have a feature layer that contains about 460,000 records and it currently takes about 20 minutes to read that table using the arcpy.da.TableToNumPyArray() tool. Is there a more efficient way to read the rows to then be able to manipulate these data? It seems like there should be a more "C-like" way to access the rows. Here is the function in it's entirety, though I'm focused on the line near the bottom where I'm calling arcpy.da.TableToNumpyArray() to read out the data:

def import_surface(surface, grid, nrow, ncol, row_col_fields, field_to_convert):
    """
    References raster surface to model grid.
    Returns an array of size nrow, ncol.
    """
    out_table = r'in_memory\{}'.format(surface)
    grid_oid_fieldname = arcpy.Describe(grid).OIDFieldName

    # Compute the mean surface value for each model cell and output a table (~20 seconds)
    arcpy.sa.ZonalStatisticsAsTable(grid, grid_oid_fieldname, surface, out_table, 'DATA', 'MEAN')

    # Build some layers and views
    grid_lyr = r'in_memory\grid_lyr'
    table_vwe = r'in_memory\table_vwe'
    arcpy.MakeFeatureLayer_management(grid, grid_lyr)
    arcpy.MakeTableView_management(out_table, table_vwe)

    grid_lyr_oid_fieldname = arcpy.Describe(grid_lyr).OIDFieldName
    # table_vwe_oid_fieldname = arcpy.Describe(table_vwe).OIDFieldName

    # Join the output zonal stats table with the grid to assign row/col to each value.
    arcpy.AddJoin_management(grid_lyr, grid_lyr_oid_fieldname, table_vwe, 'OID_', 'KEEP_ALL')

    # Take the newly joined grid/zonal stats and read out tuples of (row, col, val) (takes ~20 minutes)
    a = arcpy.da.TableToNumPyArray(grid_lyr, row_col_fields + [field_to_convert], skip_nulls=False)

    # Reshape the 1D array output by TableToNumpy into a 2D structured array, sorting by row/col (~0.1 seconds)
    a = np.rec.fromrecords(a.tolist(), names=['row', 'col', 'val'])
    a.sort(order=['row', 'col'])
    b = np.reshape(a.val, (nrow, ncol))

    return b


texts:
When i try to delete a feature class , after or during the use of an Arcgis engine app which uses the same feature class i got an error concerning the lock  , so  i can't get ride of this lock until i restart the computer.

i would like to know a way to remove locks using arcpy or arcobjects 10.1 , it would be perfect if i can get ride of the lock inside the arcengine app which also locks the geodatabase , i have seen different posts in SE about locks but no one of them gives the solution to free locks using ArcObjects.


texts: Creating a new point shapefile when vector grid features have a common attribute?

texts:
I have a very large shapefile that is composed of many cells. If there are 4 cells together that have the same attribute value, I must locate a point geometry, like below:

Each side of the cell is 625 meters long. Any ideas are welcomed (QGIS, ArcGIS, or code in PyQGIS...).


texts: Google Earth Engine: Problem with the updateMask function

texts:
I'm trying to use the updateMask function in GEE but I'm not quite sure how to use it in this occasion. So i did a trend analysis with formaTrend() to analyze the loss of forest cover in an area:

  var aoi = ee.Geometry.MultiPolygon([
  [6.753654920661802,49.98597828191218],
  [6.7461018200758645,49.902020897367365],
  [6.941795789802427,49.897598033936816],
  [6.945915662849302,49.99083449205881],
  [6.753654920661802,49.994807390209445],
  [6.753654920661802,49.98597828191218]]);

var addVI = function(image) {
  image = image.addBands(image.normalizedDifference(['B4', 'B3']).rename('NDVI'));
  return image;
}

var coll = ee.ImageCollection('LANDSAT/LT05/C01/T1_TOA')
  .filterDate('1995-01-01', '1995-12-30') // Zeitraum vor den Sturmereignissen
  .filterBounds(aoi)
  .filter(ee.Filter.lte('CLOUD_COVER', 4))
  .sort('DATE_ACQUIRED')// Collection sortiert nach Cloudcover

var coll2 = ee.ImageCollection('LANDSAT/LT05/C01/T1_TOA')
  .filterDate('1996-01-01', '1996-12-30') // Zeitraum vor den Sturmereignissen
  .filterBounds(aoi)
  .filter(ee.Filter.lte('CLOUD_COVER', 8))
  .sort('DATE_ACQUIRED')

print(coll, coll2);

var collection = coll.merge(coll2);
collection = collection.map(addVI);

// Forma Trend
var ndvistack = collection.select('NDVI');

var trend = ndvistack.formaTrend();
var ltTrend = trend.select(['long-trend']).clip(aoi);

var loss = ltTrend.mask(ltTrend.lt(0));

Map.addLayer(loss, {color: 'red'}, 'Loss');

var img = collection.median();

// img = img.updateMask(loss); doesnt work ):

After I figured out the areas which are affected by loss (stored in there variable "loss"). Now i want to calculate another vegeatation index exclusively on the areas where the algorithm detected a loss.

I tried using the updateMask() function like this:

img = img.updateMask(loss);

but i think i'm missing something, which i can't find out.
Is there a way to use the loss variable as a mask for other images?


texts: Extracting coordinates of points from grid cells

texts:
Within QGIS, I have two layers: (1) a layer of points and (2) a grid layer. Each grid cell has its own ID.

For each grid cell with greater than 10 points, I'd like to extract the coordinates for the points. At the moment, there is no common variable between the grid layer and the points to do so.

Is there a simple way to do this with the processing toolbox?


texts: Browser-based user login (ESRI dojo Javascript 3.14)

texts:
I have a web application (several custom Javascript modules plus a single HTML page) hosted on our local server. This works fine when the services that the application consumes are shared to 'Everyone' on ArcGIS.com.

To meet a business requirement I've now shared the tile and feature services that the application consumes only with a Private group on ArcGIS.com. I'm about to register the Javascript application on ArcGIS.com as a Ready to Use Web Mapping application. However, I don't understand the Redirect URL that is required as part of the registration process. In fact, I'm to sure I understand how to authenticate a user via ArcGIS.com, open our application, and consume the private group's services.

So, if our web site was ficticioussite.org.uk/sites.htm (this URL is set as the application's URL when the application is added to ArcGIS.com)  what would be the redirect URL? Is it ficticioussite.org.uk/sites.htm again?

I've studied the documentation and samples but I don't understand how to authenticate the user (each user is a member of the Private group) and consume the private services. Would the user open the application via
ficticioussite.org.uk/sites.htm and login to the private group using IdentityManager?


texts: Points within polygons in different projections

texts:
I have a polygon and a point  feature class in some projection. The points are contained within the polygons

The points are always inside the polygons no matter the projection destination I'm projecting points and polygons to (including 'unprojecting' to geographic).

Putting visually JGH example I have a big polygon in WKID =4326 with a point outside of it...

...when I project them (say in 3347 to follow JGH example) I get this...

where the point now is inside the polygon. I used the ArcGIS Pro 'project' function.

When I use project setting the optional parameter 'Preserve Shape' as in the following picture.

...I get something definitely better like...

where (hopefully) the point stays outside the polygon.
'Preserve shape' parameter did what outlined by J.R. Extra vertices that provide a better discretization of the polygon in the destination projection....otherwise everything is reduced to projecting 4 vertices and drawing straight lines in the destination space to connect them (producing the indesiderable output of the second picture above)

Am I getting this right ?


texts: libqgis_core.so.3.4.1 is missing (ubuntu 18.04)

texts:
Wile I was trying to amateurish solve the problem of that post: QGIS 3.4.1 crash without error after opening "options" (18.04. Bionic Beaver)
I got a new, much more dramatic problem:

  /usr/lib/qgis/crssync: error while loading shared libraries:
  libqgis_core.so.3.4.1: cannot open shared object file: No such file or
  directory

I purged, autoremoved, installed and uninstalled, changed the repositorys, reconfigured and so on... But I have no glue, what to do next. I just need to get QGIS3 running and have no time to re-setup my system!

Where do I get the libqgis_core.so.3.4.1? Google isn't telling me...

UPDATE: result after trying to install via 'sudo apt-get install qgis python-qgis qgis-plugin-grass':

>     Setting up qgis-providers (1:3.4.1+28bionic) ...
>     /usr/lib/qgis/crssync: error while loading shared libraries: libqgis_core.so.3.4.1: cannot open shared object file: No such file or
> directory
>     dpkg: error processing package qgis-providers (--configure):
>      installed qgis-providers package post-installation script subprocess returned error exit status 127
>     dpkg: dependency problems prevent configuration of qgis:
>      qgis depends on qgis-providers (= 1:3.4.1+28bionic); however:
>       Package qgis-providers is not configured yet.
>
>     dpkg: error processing package qgis (--configure):
>      dependency problems - leaving unconfigured
>     dpkg: dependency problems prevent configuration of qgis-plugin-grass:
>      qgis-plugin-grass depends on qgis (= 1:3.4.1+28bionic); however:
>       Package qgis is not configured yet.
>
>     dpkg: error processing package qgis-plugin-grass (--configure):
>      dependency problems - leaving unconfigured
>     No apport report written because the error message indicates its a followup error from a previous failure.
>                                                                                No apport report written because the error message indicates its a
> followup error from a previous failure.
>                                                                                Errors were encountered while processing:
>      qgis-providers
>      qgis
>      qgis-plugin-grass
>     E: Sub-process /usr/bin/dpkg returned an error code (1)


texts: ApplySymbologyFromLayer works in ArcMap but not in python script

texts:
So the ApplySymbologyFromLayer_management() tool doesn't seem to like me, at least when I try to use it in a script. I have a Field_Boundary2.lyr file that has a "hollow" fill color and a black outline for a polygon boundary. In ArcMap, when I import this layer's symbology onto another layer, it looks like this:

But if I try to use it with arcpy.ApplySymbologyFromLayer() in a script, it gives me this:

Here's the code:

        df = arcpy.mapping.ListDataFrames(mxd,"dfElevation")[0]
        fc1 = arcpy.mapping.Layer(fcGrid)
        arcpy.ApplySymbologyFromLayer_management(fc1, arcpy.mapping.Layer(symbologyTemplatePath + "Elevation.lyr"))
        arcpy.mapping.AddLayer(df,fc1,"BOTTOM")
        fc2 = arcpy.mapping.Layer(fcBound)
        arcpy.ApplySymbologyFromLayer_management(fc1, arcpy.mapping.Layer(symbologyTemplatePath + "Field_Boundary2.lyr"))
        arcpy.mapping.AddLayer(df,fc2,"TOP")

From the picture above you can see that the "Elevation.lyr" symbology is being imported correctly, but not the "Field_Boundary2.lyr". Any idea why this might be happening?


texts: Exporting rasters in GCS in Google Earth Engine

texts:
I have a raster that covers an entire country with the base data having a 0.05° resolution. The country spans two UTM zones, so I am reluctant to reproject the file to a projected coordinate system. The issue with this is that exporting the raster in Google Earth Engine requires the user to set a scale, which defaults to 1000m per pixel. This doesn't really work for a raster with a resolution that depends on the latitude. Is there a work-around for this issue? Below is an excerpt script of my data (though the AOI in the larger script is much bigger):

var aoi = ee.Geometry.Polygon(
        [[[101.17563162758336, 18.896757345850304],
          [101.15365897133336, 18.70955737155526],
          [101.39535819008336, 18.69394798887354],
          [101.39535819008336, 18.912347927872148]]]);

var chirps = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY")
            .filterDate('1984-06-01','1984-06-30')
            .filterBounds(aoi)

var chirpsMean = chirps.mean().clip(aoi)

Map.addLayer(chirpsMean)

Export.image.toAsset({
  image: chirpsMean,
  description: 'LaosExample',
  //scale: 5000,
  region: aoi
})


texts: How to translate unknown old coordinates from a 1934 map of Tanzania?

texts:
I have a old map depicting Tanzania, dated 1934, with X and Y coordinate values. I wish to translate these locations precisely into latitude and longitude (degrees, minutes, and seconds), or any international coordinate system. Any suggestion will be most welcome.

My XY coordinate values look like this:

X = 12106.7

Y = 72776.2

Those particular coordinates should be located near 3° 21' S / 35° 29' E in lat/long.


texts: How to calculate "water accumulated in each cell" using QGis?

texts:
I create channel network and catchment area from a DEM with SAGA GIS. I know that flow is starting from a cell to an other cell. Firstly, a cell is filled with water and then it is starting to flow into an other cell.

"Flow accumulation" create accumulated flow into each cell. It is about  number of cells.

But I want to know, How can I determine how much water accumulate in cells? Do I need a volume layer to solve this?


texts: WFS: DescribeFeatureType query with typenames?

texts:
I am struggling with the DescribeFeatureType query with WFS KVP encoding.

I have this WFS: https://inspire-wfs.maanmittauslaitos.fi/inspire-wfs/au?request=GetCapabilities&service=wfs&version=2.0.0 and I can see in the GetCapabilities that there is a FeatureType called au:AdministrativeUnit.

I could of course use DescribeFeatureType to give me the description of all Feature Types of the service, but I would like to use &typeName=XXXX in the kvp to choose only this FeatureType au:AdministrativeUnit

So I basically do the following request: https://inspire-wfs.maanmittauslaitos.fi/inspire-wfs/au?&request=DescribeFeatureType&service=WFS&version=2.0.0&typeName=au:AdministrativeUnit

But the result is not filtering at all to give me only the Feature Type with the name au:AdministrativeUnit.

I saw that in the OGC specifications, there is sometimes written typename or typenames. I tried both of them with changing the version and I don´t get it!

Am I missing something important here?

With an other WFS, I can do the query with filtering typenames without issue!


texts: Web-based thumbnail image in 'Item Description' for a Python toolbox [ArcGIS 10.1]

texts: I'm using locally stored images as thumbnails in the Item Description for a python toolbox [ArcGIS 10.1].  I'd like to load the images from an online source so other users of the toolbox can still see the images.  Is there a way to do this?


texts: Permanent add Geoserver Server on Geoexplorer different URL

texts:
I have my Geoexplorer run on Open shift (asume http://geoexplorer-aaa-rhcloud.com) using this Git (https://github.com/jason-callaway/geoexplorer-on-openshift). My geoexplorer run perfectly,

Problem is my geoserver are on another gear on openshift (asume http://geoserver-aaa-rhcloud.com/web). I have Github File for making Geoexplorer.war (https://github.com/boundlessgeo/suite/tree/r4.1/geoexplorer)  and I know where should I put my new "Edited" Geoexplorer.war on 1st Git to push on open shift.

My question is how to make my geoexplorer permanently add my geoserver url, at which part of the Git source I should add / edit so that I dont have to repeat step "Add New Server" each time I open / reload my geoexplorer site?


texts: How to have more than one output on iterate tables (modelbuilder)?

texts:
I'm doing a model in modelbuilder where I have as input 5 excel tables and what I want is to have 5 outputs in the end of the model. How can I do that?

I want to have 5 outputs after Project tool.
The model currently is like this:


texts: How to dock the "Identify" panel to the top or bottom of the QGIS window?

texts:
I'd like to dock the "identify" window to the top or bottom of the QGIS window. Docking to the left and right seems to work fine (a blue area indicates the new position of the panel after releasing the mouse button). However, when dragging the panel to the top or bottom of the window, no blue area is shown.

However, docking other panels (e.g. "View Results") to the top or bottom works?
I'm using QGIS 3.2.3.


texts:
I have clipped a raster file using a shapefile as a mask.

And I tried to export the new layer and modified the projection from WGS 84 EPSG:4326 to EPSG: 102022 Africa Albers Equal Area Conic.
When export was processing QGIS suddenly crashed.
What went wrong?


texts: Creating a custom widget like Form Aware Value Relation Widget

texts:
I'm trying to get to know how to make my own widgets to be used in the Edit Widget.
I have looked at this plugin here and it is very usefull:
https://github.com/elpaso/qgis-formawarevaluerelationwidget

I'm trying to make a custom widget which will have a combobox and a few small buttons (like flash feature, zoom to, pan to)

In the example i have created from the form aware relation widget i can succesfully create a combobox.

class FormToolsRelationWidgetWrapper(QgsEditorWidgetWrapper):

    def __init__(self, vl, fieldIdx, editor, parent):
        """
        QgsVectorLayer* vl, int fieldIdx, QWidget* editor, QWidget* parent
        """
        some init stuff here...

    def createWidget(self, parent):
        #self.mComboBox = TestWidget(parent)
        #self.mComboBox = TestWidget()
        self.mComboBox = QComboBox(parent)
        QObject.connect( self.mComboBox, SIGNAL( "currentIndexChanged( int )" ), self, SLOT( "valueChanged()" ) )
        return self.mComboBox

The result is this:

Looks very good. And i know i have not yet populated the combobox with data.

If i instead activate a widget i have created with this code (with TestWidget(parent) ):

class TestWidget(QWidget):
    def __init__(self, parent=None):
        super(TestWidget, self).__init__(parent)

        # create a combobox
        combobox_1 = QComboBox(self)
        combobox_1.addItem("Item 1")

        # create a combobox with a parent widget
        combobox_2 = QComboBox(self)
        combobox_2.addItem("Item A")
        combobox_2.addItem("Item B")
        combobox_2.addItem("Item C")

        # layout
        layout = QHBoxLayout(self)
        layout.addWidget(combobox_1)
        layout.addWidget(combobox_2)
        layout.addStretch()

        # sets the layout
        self.setLayout(layout)

I get his result:

I'm not sure what is wrong. I have test by not assigning QComboBox and QHBoxLayout to self. Without any luck.
Here the widget does not replace the ID field. i thought this would be replaced with my widget. Like the example with qcombobox.

Any suggestions to make this work?

I can not find any examples/tutorials explaning this. Maybe if it get this sorted i will make a small introduction to this.


texts: Is nearest search done on Euclidean plane accurate

texts:
I am writing a program which will find all the points within 1km. I am planning to to do this on a Euclidean plane as my database only supports Euclidean plane or sphere indexing (http://docs.mongodb.org/manual/applications/geospatial-indexes/). I chose the plane due to performance. I was wondering how accurate using this plane would be all over the earth. I am not too considered about EXACT matches, but I would like no more than 50m of error when searching for all points within 1km of a certain point. Would this work all over the earth or would I lose accuracy if I move throughout the earth (poles and stuff)?

PS: I am not expert at GIS stuff, wouldn't even consider myself a beginner.


texts: Clip/Intersect but not split features that are on the border

texts:
When I use the clip or intersect tools in QGIS it splits the polygon clipper features that are on the border of the clipper, which produces polygons that are only part of the real polygons.

Is there a tool that allows me to specify that features that are on the border are either in or out of clipping result set?


texts: Aggregating multiple time series to match MODIS timestep in R

texts:
I have multiple time series data for precipitation and evapotranspiration that are irregularly spaced.  All but one of the datasets are daily products, but sadly one is a MODIS product that has irregular spacing (8 days most of the time, 5-6 days at the end of the year depending on whether it's a leap year).  But I need all products to match the 8 day temporal resolution of the MODIS product.  For instance:
dates        prod1     prod2
             (daily)   (variable)
10-01-2001   1          NA
10-02-2001   1          NA
10-03-2001   1          NA
10-04-2001   1          NA
10-05-2001   1          NA
10-06-2001   1          NA
10-07-2001   1          NA
10-08-2001   1          8

What I am looking for is a way to consistently get the aggregate of product 1 to match prod2, but prod2 has a variable number of NAs depending on the period.
Obviously, rollapply and aggregate are out because I don't have consistent spacing. Is there another way to do this without a very complicated ifelse statement/mutate combination?


texts: Running GDAL function in R to read File Geodatabase

texts:
I was able to install GDAL through Conda and run readOGR with the following steps.
Programs needed:
R (download available here with RStudio: https://rstudio.com/products/rstudio/download/)
Anaconda3 (download available here: https://www.anaconda.com/products/individual)
You will need to install GDAL in order to run the R functions to open the .gdb files. One way to install GDAL is through Conda, which I will describe here:
Installing GDAL in Conda:
Install Anaconda3 at the linked provided up.
-Search for Anaconda Navigator in windows search box.
-Click on Environments on the left hand side of the Anaconda Navigator.
-Create a new Environment with the + buttom on the bottom left.
-On the top bar of the environment, click on the drop bown and select "not installed".
-Search for gdal in the "not installed" packages.
-Install r-gdalutils, r-rgdal, gdal, and libgdal. (I'm not sure if all of these are actually necessary, I just installed anything with the name gdal). The installation can take a few minutes.
The readGDBlayer function worked for the 2014 .gdb files on Marine Cadastre (where the sample data originated), but was unable to open .gdb files on 2009 - 2013. I am not sure why this is. The readOGR function was able to open the .gdb files from 2009 - 2013
Here is a link to the 2009 AIS data that you can download https://coast.noaa.gov/htdata/CMSP/AISDataHandler/2009/index.html.
Does anyone know why readGDBlayer function will not work with this data?
read_GDB_Layer <- function(dsn, layerName, overwrite = T) {
   conversionDir <- tempdir()

   gdalUtils::ogr2ogr(src_datasource_name = dsn,
                      dst_datasource_name = conversionDir,
                      f = "ESRI Shapefile", layer = layerName,
                      verbose = T, overwrite = overwrite)

   df <- foreign::read.dbf(file.path(conversionDir, paste0(layerName, ".dbf")))

   return(df)
}


texts:
I am trying to talk to the PyQGIS API for use in a 3rd party application using the offical .dmg image provided on the QGIS homepage.

After setting PYTHONPATH to /Applications/QGIS3.4.app/Contents/Resources/python and initializing everything as stated in https://docs.qgis.org/3.4/en/docs/pyqgis_developer_cookbook/intro.html#using-pyqgis-in-custom-applications, I am facing the following error when running from PyQt5.QtGui import or from qgis import *:

  ImportError: dlopen(/Applications/QGIS3.4.app/Contents/Resources/python/PyQt5/QtGui.so, 2): Library not loaded: @executable_path/../Frameworks/QtGui.framework/QtGui
  Referenced from: /Applications/QGIS3.4.app/Contents/Resources/python/PyQt5/QtGui.so
  Reason: image not found

Looks like I need to set some linkers or similar?

Side facts:
- I have no problems when using the homebrew version of QGIS).
- I am on macOS 10.14.6 (Mojave)


texts: Appending interpolated or nearest raster cell value to vector point on a cell with no data?

texts:
I have been searching for a solution to this problem, but was not satisfied by the answer given in Finding nearest raster cell value based on vector point?

I have some points that land on cells with NoData in a raster of continuous data, and would like to assign values to the points based on either the closest cell with data or calculate a value with bilinear interpolation from neighboring cells with data.  However, using tools in ArcGIS 10.2 Spatial Analyst such as Extract Multi Values to Points, I get a NoData result whether the bilinear interpolation of values is checked or not, despite neighboring cells having data, apparently because the tool is not overlooking the NoData cell.  Is there a way to ignore the NoData cell usign this method or another way of getting values for these points?

Here is a screenshot example of what I'm looking at, where white cells have no data:


texts: Hide/Show Leaflet Routing Machine when click on button (Easy Button plugin)

texts:
I tried to put this code :
var rlayer = null;

L.easyButton('fa-level-up',
    function () {
        if (rlayer) {
            map.removeLayer(rlayer);
            rlayer = null;
        } else {
            var routing = L.Routing.control(...);
            rlayer = L.layerGroup([routing]);
            map.addLayer(rlayer);
        }
    }, 'Display Route').addTo(map);

on my map but when I click on the button, the console tells me 2 errors :

Uncaught TypeError: Cannot read property 'call' of undefined
Uncaught TypeError: Cannot read property 'appendChild' of undefined


texts: Use PyQGIS With IronPython : No module named qgis.core

texts:
I need to create a Qgis python script called by IronPython for application purpose. I've already read Using PyQGIS in standalone scripts & Running application of Qgis Cookbook (https://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/intro.html).

After set environment path & run my script, it gave me an error :

  IronPython.Runtime.Exceptions.ImportException : 'No module named
  qgis.core'

  IronPython.Runtime.Exceptions.ImportException   HResult=0x80131500
  Message=No module named qgis.core   Source=Microsoft.Dynamic
  Arborescence des appels de procédure :    à
  Microsoft.Scripting.Runtime.LightExceptions.CheckAndThrow(Object
  value)    à
  IronPython.Runtime.Operations.PythonOps.ImportStar(CodeContext
  context, String fullName, Int32 level)    à
  Microsoft.Scripting.Interpreter.ActionCallInstruction`3.Run(InterpretedFrame
  frame)    à
  Microsoft.Scripting.Interpreter.Interpreter.Run(InterpretedFrame
  frame)    à
  Microsoft.Scripting.Interpreter.LightLambda.Run2[T0,T1,TRet](T0 arg0,
  T1 arg1)    à
  IronPython.Compiler.PythonScriptCode.RunWorker(CodeContext ctx)    à
  IronPython.Compiler.PythonScriptCode.Run(Scope scope)    à
  IronPython.Compiler.RuntimeScriptCode.InvokeTarget(Scope scope)    à
  IronPython.Compiler.RuntimeScriptCode.Run(Scope scope)    à
  Microsoft.Scripting.SourceUnit.Execute(Scope scope, ErrorSink
  errorSink)    à
  Microsoft.Scripting.Hosting.ScriptSource.Execute(ScriptScope scope)
  à Microsoft.Scripting.Hosting.ScriptEngine.ExecuteFile(String path,
  ScriptScope scope)    à
  LivrableChecker.Checker.FileDepencyChecker.QGISChecker() dans
  C:\Users\BVCS5568\Desktop\Livrable Checker\3 - Source
  code\LivrableChecker\LivrableChecker\Core\Checker\FileDepencyChecker.cs
  :ligne 141    à ConsoleApp1.Program.Main(String[] args) dans
  C:\Users\BVCS5568\Desktop\Livrable Checker\3 - Source
  code\LivrableChecker\ConsoleApp1\Program.cs :ligne 71

Utilities: Ironpython 2.7, Qgis 3.0, Visual Studio

Here is my code For C# call and Python script:

https://pastebin.com/FpnncwSj


texts: How to calculate raster statistics for each values in a band?

texts:
I have one raster scene that contains one band with multiple pixel values (0-14).  Each values represent deforestation in certain year (0=no tree cover loss, 1=loss in 2001, 2=loss in 2002, and so on) and I need to calculate the deforestation statistics for each years (e.g. statistics for each pixel values in that scene).

I've tried using Zonal Statistics, but it only gives me a lump of all values (total count in zonal statistic = sum of all values). What I need is to count how many pixels have value '1', '2', and so on. Any suggestion on what tool(s) to use?


texts: Ways to visualize multiple large LiDAR tiles?

texts:
I want to load approximately 10 large LiDAR tiles for analysis. I have tried using ArcScene but it appears there is a limit on how much data ArcScene can handle at one time. Each tile is 65-90MB in size. The only returns that I am using (and that are available) are ground returns.

My computer has 8GB ram, an amd 6 core processor, and a really good video card. I am fairly certain that it is more than sufficient to handle large amounts of data.

I am not against using something other than an ESRI product to complete this task. What are other options to visualize multiple large LiDAR tiles?

For visualizing LiDAR/point cloid data in web browsers, see Viewing LiDAR data from web browser?


texts:
I'm aware of several books on the history and developement of GIS, but not one specifically about SYMAP.

SYMAP was developed in 1963 and is considered the first automated computer mapping program.

Has there been a book written specfically about the development of SYMAP, it's impact and influence?


texts: Exporting a raster object with a categorical legend

texts:
I created a raster object (package raster in R) using the function rasterize starting from a shapefile of polygons. I passed as field a shapefile attribute (CODE) containing 37 landuse categories.
This is what I did in code:
myrast <- rasterize(x = shapefile, y = another_raster, field = shapefile$CODE)

The command converts the categorical variable into numbers from 1 to 37 and correctly create the raster.
Now I want to export the raster to see it in a GIS but I would like that a legend would be shown with the CODE and not with the number from 1 to 37.
How can I do this?


texts: Editing featureclass in Oracle SDO_GEOMETRY

texts:
I'm building an empty featureclass in Oracle and want to enable editing in arcmap/py with having the featureclass user SDO_GEOMETRY for the shape column. However, I'm running into difficulties on how the information is being stored. I'm receiving this error when trying to edit a few vertices and save.

The idea here is that Oracle could populate this table's SDO_GEOMETRY column from other sources/tables/applications and then this table could be viewed in ArcMap and also edited if needed.
The featureclass is Registered as Versioned.


texts:
I am trying to convert an ArcGIS based GIS setup to Open Source - QGIS.

In my old setup I had a lot of domains on an ArcSDE Database. I am struggling a bit to find a similar functionality in QGIS.

I have tried Value Maps, but they seem to be project file specific so whenever I am opening my project file on a new PC I have to define them again. This is not working for me as I am setting up a solution where more than one person needs to access the same data and work on it.

Another possible solution is Value Relations. I have succesfully created a table with my "domain values" and joined these to the layers in QGIS, but one big problem is that when I try to label my features they only show the key coloumn eg. a number. Instead I want to see what has been written in the value coloumn.

Does anyone know how to fix this? Or perhaps know a better equivalent to domains in QGIS?


texts: Saving an output file during merging on QGIS

texts: I am using QGIS for Landsat 7 data and under raster> miscellaneous > merge I am trying to save the output file, however I don't know where to save it as everything I press asks me to 'open' it rather than 'save' it.


texts: Opening PDF w/ action based on variable crashes QGIS

texts:
I found this script which I modified a bit to open a pdf based on a page number given in the attribute list. Upon trying to execute the action it works just find but QGIS (3.4.15) is crashing!
Any ideas why QGIS is crashing or if there is something not working properly within the code?
My modified code:
import subprocess
import os
path_to_pdf = os.path.abspath(r'[% "PDFpath" %]')
path_to_acrobat = os.path.abspath(r'C:\Program Files (x86)\Adobe\Acrobat Reader DC\Reader\AcroRd32.exe')
process = subprocess.Popen([path_to_acrobat, '/A', 'page=[%PDFpage%]', path_to_pdf], shell=False, stdout=subprocess.PIPE)
process.wait()


texts: Downloading all images in time interval with WMS getmap request

texts:
I'm trying to download sentinel-2 satellite images from sentinel-hub (https://www.sentinel-hub.com/) via WMS requests. The request allows me to specify a time interval with the time parameter but the I don't seem to be getting all the images in the specified time interval.
The problem I'm trying to solve is the following:
Given a bounding box (latmin,longmin,latmax,longmax) I would like to get all the satellite images taken in the box from a date yyyy-mm-dd to another YYYY-MM-DD. The returned object should be some sort of an array of images, ideally with a timestamp on every image.
I'm building a wms getmap request with python and this is what I've already done.
dist = 5
#i = 21800
for i in range(0,len(fireBC)):
    if fireBC["FIRE_YEAR"][i] >= 2018:

center = np.array(fireBC['geometry'][i].representative_point())
date = fireBC['FIRE_DATE'][i]
firedateB = dt.datetime.strptime(date, '%Y%m%d')
firedateE = firedateB + dt.timedelta(days=10)
firedate_str = firedateB.strftime("%Y-%m-%d") + "/" + firedateE.strftime("%Y-%m-%d")

fire_imgnew = wms.getmap(
                    layers=['SWIR'],
                     styles=['default'],
                     srs='EPSG:4326',
                     bbox=boundingBox(center[0],center[1],dist),
                     size=(512, 512),
                     #geometry = ade,
                     format='image/png',
                     #Transparent = False,
                     time=firedate_str)

print(type(fire_imgnew))
decoded = cv2.imdecode(np.frombuffer(fire_imgnew.read(), np.uint8), -1)
cv2.imwrite("Images/fire"+str(i)+".png", decoded)

FireBC is a geopandas files of all historic forest fires in British Columbia Canada.
The number of requests that I can make is limited.
Is there no way to get an array of images with one request?


texts: Decimal disappears after calculation with CalculateField_management

texts:
I try to change the field name with the following steps:

Add new fields
Calculate new fields by making them equals to old ones
Delete old ones

As the new fields are not fixed, for example, some features may have three fields and some may have five, I cannot refer to the exact field names. Therefore, I use a for loop to iterate fields. It can run but the result deletes the decimal comma, making 100,101 to 100101 for example.
Here is the code I have. To make it simple, I only attached the part related to the field calculation.
for items in field_names:
    if len(items) <= 7:
        new_field_name = items+str(n)+'b'
    else:
        new_field_name = items[:7]+str(n)+'b'
    arcpy.AddField_management(track_point, new_field_name, "DOUBLE")
    arcpy.CalculateField_management(track_point, new_field_name, '\'!'+item+'!\'', "PYTHON_9.3")
    arcpy.DeleteField_management (track_point, item)

If I use the field name directly, everything is Ok. Can anyone tell me how to fix this problem or is there another better way to change the field name?


texts:
When using the following PyQGIS line, I write a layer as a CSV into an output path.

QgsVectorFileWriter.writeAsVectorFormat(lyr, outputPath, 'utf-8', lyr.crs(), "CSV",
                                        layerOptions=['GEOMETRY=AS_XYZ'])

I would like to write it as type "TXT". What are alternative ogr drivers than "CSV"? Perhabs there is a list of all available ogr drivers?

Unfortunately, the PyQGIS Documentation doesn't give me any options.


texts: Create line feature from point feature attributes: azimuth and distance

texts:
I've used Shape Tools plugin's Azimuth Distance Digitizer tool, to point and click, and then enter azimuth and distance in order to get a line.

This case is the same, but I have many features, so it would be nice an automated way to create these line features.

I've read Draw line with lon, lat, distance and azimuth

Can someone point me in the right direction to make it with QGIS?


texts: PostGIS (or rpostgis) issue on Ubuntu 18.04

texts:
I've installed PostgreSQL 10 and PostGIS on Ubuntu following instructions here. Now I'm migrating an R interface from Windows to said Linux box and have an inconsistency I can't sort out.

When I pgInsert a normal table I have no trouble: the first insert creates my table, second adds to the table, etc.

When I try to pgInsert a geometry object into a new table, PostGIS fails, complaining:

  Creating new table... Error in postgresqlExecStatement(conn,
  statement, ...) :    RS-DBI driver: (could not Retrieve the result :
  ERROR:  relation "spatial_ref_sys" does not exist LINE 1: SELECT srid
  FROM spatial_ref_sys WHERE auth_name = 'EPSG' AN...

Seems there is a subsystem containing CRSs which I am missing?


texts: How to define the type of a field in the field calculator so that the rasterize (vector to raster) tool recognize it?

texts: I have a shapefile (polygon) that I want to rasterize by an attribute but the rasterize tool does not recognize the right field. When I create a new field using the field calculator it does not matter what type I choose (integer, real, ...) the result is qlonglong (integer64), a type that the rasterize tool apparently does not recognize.


texts: How do I fix "name 'iface' is not defined"?

texts:
i try to export a map but at the end of my process I always get the following message: "name 'iface' is not defined".

What is wrong with my script?

I am working with QGIS 2.18.12

from PyQt4.QtGui import QColor, QImage, QPainter
from qgis.core import QgsExpression, QgsStatisticalSummary, QgsSymbolV2,   QgsDataDefined, QgsRuleBasedRendererV2, QgsMarkerSymbolV2
from PyQt4.QtCore import Qsize

...

map_settings = iface.mapCanvas().mapSettings()
c = QgsComposition(map_settings)
c.setPaperSize(400,160)
c.setPrintResolution(dpi)

...


texts: On QGIS Composer how to filter the atlas coverage layer using a project variable (rather than just a value)?

texts:
A QGIS Composer atlas can be generated from a subset of the coverage layer features using the coverage layer filter control. I can filter for coverage layer features by comparing an attribute of the features to a static value, but not by comparing to a project variable (with the same value). Can this be done? What is the syntax of the expression? Is it just not possible for some reason? Is there a work-around?

Specific situaion

I have QGIS print composer set up to generate an atlas. This uses a coverage layer 'Locations'. In the 'Locations' data is a field 'Project'. I want to print atlas pages for any feature in Locations where "Project" = '1'. Entering "Project" = '1' into the atlas coverage layer filter field works fine - I get atlas pages.

So far so good, but I'd like to be able to change the value '1', using a project variable @myvariable in its place.

Entering "Project" = @myvariable into the filter field results in no atlas pages (an error appears) even if @myvariable has the value 1. I've tried numerous changes in syntax without success.

Images to clarify question:

A Composer atlas can use select features (rather than all features) from the coverage layer - using a filter - as shown (the filter is set to "ProjectID"='1'):

Using a variable in the filter seems not to work - the error is shown here when filtering using ("ProjectID"="@currentprojectUPK):

It can be seen that the variable @currentprojectUPK has the value '1' in the list of variables as shown here:

I wondered if somehow the Composer couldn't access the variable values, but this can be used elsewhere (for example in a text box) successfully. I also wondered if there was an issue with the variable being a string not an integer, but changing it to an integer in the statement didn't help - and anyway I'm using the value successfully elsewhere.