texts: Avoiding overlapping charts on a map with Leaflet?

texts:
I'm currently developing an online GIS with Leaflet and need to plot bar charts on a map. I'm using leaflet.minichart (https://rte-antares-rpackage.github.io/leaflet.minichart/) for the charts.
My problem is that some of the bar charts are overlapping - usually due the proximity of the centroids of my polygons. You can see a good example in the attached printscreen on the darkest polygon.
I would like to avoid such overlaps and at this stage I'm fairly agnostic to the type of solution - I just want to ensure a nice rendering along with user-friendliness. Here below my code:

function add_poly_charts(features){
            var feature; // contains the polygon centroid and name, while polygon_charts is an array containing the data for each polygon

            // size simply changing with the number of charts, this si what we woudl like to automate to avoid collisions
            var size=80;
            if(n_poly_charts>5){
                size=30;
            }

            for (var i=0; i< features.features.length; i++){
                feature=features.features[i];
                var chart=L.minichart([feature.properties.centroid_y,feature.properties.centroid_x],
                                      {data: polygon_charts[feature.properties.adm],type:"bar",height:size,width:size, labels: names_list,
                                      colors:chart_colors});

                poly_charts.push(chart);
                mymap.addLayer(chart);
            }

        }


texts: Unable to add PostGIS to database due to undefined symbol: proj_get_type

texts:
I am setting up a PostgreSQL10 server with pgRouting on Fedora 30 VM. by running these command as root:

dnf install https://download.postgresql.org/pub/repos/yum/reporpms/F-30-x86_64/pgdg-fedora-repo-latest.noarch.rpm

dnf install postgresql10

dnf install postgresql10-server

/usr/pgsql-10/bin/postgresql-10-setup initdb
systemctl enable postgresql-10
systemctl start postgresql-10

dnf install https://kojipkgs.fedoraproject.org//packages/pgRouting/2.6.2/2.fc30/x86_64/pgRouting-2.6.2-2.fc30.x86_64.rpm

When I trying to add the first extension to the database, it failed:

create extension postgis;

>ERROR:  could not load library "/usr/pgsql-10/lib/postgis-3.so": /usr/pgsql-10/lib/postgis-3.so: undefined symbol: proj_get_type
SQL state: XX000

It seems that it is postgis libraries linking problems. However, I did't have any legacy libraries installed prior this, so I starts getting lost here. Attempting to update the packages mentioned (proj, geos) did nothing, removing them via dnf would also remove postgis as dependency.

I am not sure that I had to manually build and install everything one by one or possible to solve it by editing some configure file sitting somewhere unknown. I never tried to Build a RPM Package from Source.


texts: Change coordinate system to add buffer in meters

texts:
I have am trying to add functionality to allow users to find points that are within a buffer distance inside or outside the parameter of a polygon. The map is in IProjectedCoordinateSystem so to my understanding the the ITopologicalOperator.Buffer is in degree units. I want to be able to use the users provided buffer in meters.

I figured I could project the polygon to another coordinate system, apply the buffer in meters and then convert back. This does not seem to be working correctly. Is this the correct or only way to get a buffer in meters when working in degrees?

private static IGeometry AddBufferToSearchFromGeometry([NotNull] IGeometry originalShape, double buffererSizeInMeters)
{
    IGeometry projectedShape = (originalShape as IClone).Clone() as IGeometry;
    ISpatialReferenceFactory3 spatialReferenceFactory = new SpatialReferenceEnvironmentClass();
    IProjectedCoordinateSystem projectedCoordinateSystem = spatialReferenceFactory.CreateProjectedCoordinateSystem((int)esriSRProjCSType.esriSRProjCS_NAD1983UTM_11N);
    projectedCoordinateSystem.SetFalseOriginAndUnits(-180, -90, 1000000);
    projectedShape.Project(projectedCoordinateSystem);
    IGeometry5 bufferedGeometry = (projectedShape as ITopologicalOperator).Buffer(buffererSizeInMeters) as IGeometry5;
    bufferedGeometry.Project5(originalShape.SpatialReference, 0);

    return bufferedGeometry;
    }


texts: Opening location folder for data files using QGIS?

texts:
I am working in a corporation that publishes maps with ArcGIS but I use QGIS for working with data.  I need to find the data files that show the information I want to publish and quickly send them via email to the office that has the ArcGIS facility.  Using the layer properties, I can see where the data is located, but there is no open folder button, which is what I am looking for.  I use Win 10, so Explorer will not open the path that I cut and paste from QGIS unless I replace the back slashes with forward slashes.

Is there an easy way to open the source folder from within QGIS 2.18.2?


texts: Image Mosaic JDBC plugin for multiple tiles

texts:
I want to publish a large raster (10000 km², 0.5m resolution). My original data is several tiles of 5km².

I've merge all tiles and published an unique raster file with Image Mosaic JDBC Plugin. But is it possible to indicate in my mosaic table in postgres several tiles ? And would it be more efficient than publish an unique raster file, even if it pyramidize ?

I tried with 2 tiles, but just the first was published, even if the coordinates was computed on the mosaic table.


texts: PostGIS -> LineString + PointZ = LineStringZ - Misjoined segments

texts:
I have used modified version of a querry from this post (I work on linestrings) and have came across few issues:

if I run the querry on my table it doesn't join the segments back in the right order
if I modify the querry modify ORDER BY from multiline_seq to multiline_seq, lineseq it joins the segments back in the right order but creates ca. 25000 features rather than expected ca. 1500
if I limit the noumber of rows with WHERE statement it works fine

MY QUERY:

DROP TABLE IF EXISTS util_div_lines_3D;

CREATE TABLE util_div_lines_3D AS

WITH
    util_div_lines_vertices AS (
            SELECT
                    gid,
                    dxf3d,
                    layer,
                    color,
                    ROW_NUMBER() OVER() AS multiline_seq,
                    (ST_DumpPoints(geom)).path[1] AS line_seq,
                    (ST_DumpPoints(geom)).geom AS geom
            FROM util_div_lines
    )

SELECT DISTINCT
        gid,
        dxf3d,
        layer,
        color,
        multiline_seq,
        ST_MakeLine(
            ST_SetSRID(
            ST_MakePoint(
            ST_X(geom),
            ST_Y(geom),
            nearest_survey_point.elevation
            ),
            27700)) OVER (PARTITION BY multiline_seq ORDER BY multiline_seq) AS geom
FROM util_div_lines_vertices
JOIN LATERAL (
    SELECT
            (
            CASE ST_DWithin(geom, util_div_lines_vertices.geom, 0.01)
            WHEN TRUE THEN elevation
            ELSE -999
            END
            ) AS elevation
    FROM util_div_survey_points
    ORDER BY geom <-> util_div_lines_vertices.geom
    LIMIT 1
) AS nearest_survey_point
ON true
WHERE dxf3d = 'YES'

EXPECTED OUTPUT

OUTPUT FROM MY QUERRY


texts: ArcGIS tool with the same result like r.sun in GRASS

texts: I'm interested if there is any ArcGIS tool which give us the same result like GRASS command r.sun. I need to get drop shadows. You can suggest some other software?


texts: Stopping Ribbon from Auto Collapsing in ArcGIS Pro

texts:
In ArcGIS Pro 2.7 how do you stop the ribbon from collapsing every time you click outside of the ribbon?
For example, to edit a polygon I must click the Edit menu option, this opens the ribbon and I can select the selection tool.  Then I can select my polygon but the ribbon collapses. To select an editing tool I have to click Edit again to open the ribbon and select some editing tool.  This results in a lot of extra clicking. I do not see a setting in the Project Options to fix the ribbon.


texts: Aligning line with closest line segment in PostGIS

texts:
I've got a road network dataset and a dataset that consist of multiple disconnected lines of varying lengths in PostGIS.

Most of these lines are (largely) parallel to the road network, but the distance between these lines and the road network can differ. Since these lines tell something about the condition of the road segments that are parallel to these lines, I would like to create a new line segment that is aligned with the closest road segment. The blue lines in the example below is what I eventually would like to achieve.

Does anyone have an idea how to solve this in PostGIS?


texts: How to create a heatmap from overlapping polygons using ArcGIS Server JavaScript API?

texts:
Using the ArcGIS Server JavaScript API version 4.18, how can I visually highlight areas with high concentrations of overlapping polygons, such as these historical cyclone tracks?

In addition I need to give higher weighting to more severe cylones based on a numerical attribute (HazardValue, as shown in the legend).
The ArcGIS Server JavaScript API contains a HeatMapRenderer, which looks promising but is only suitable for points.
If possible I would like to find a front-end solution rather than creating a raster on the server (I plan to add the ability to filter the polygons and would like the heatmap to regenerate dynamically based on the selection, which seems more suitable to front-end processing).


texts: Token doesn't work when applied to map service - Silverlight

texts:
I am trying to implement a Silverlight Web Mapping application with ArcGIS server and I'm trying to use the Token authentication with active directory as the user/roles store.

I have the service running so that I can successfully obtain tokens but when I try to put them in a map service the map service fails with the 'Unauthorized Access' exception.

I have made sure that the account is in the permissions for that service.

I've been on the Esri forums and I've found about 4 other posts with people with the same problem, but no solutions.


texts: Zonal Statistics calculator differences between QGIS Lyon and Wien

texts:
I am trying to implement a process that uses the Zonal Statistics menu function. In Lyon, the menu that pops up had a field with a list of calculations, each with a check box. Here you can select what to calculate, specifically the mean of values in a polygon. However, in the Wien version, the zonal statistics function does not have the same field.

How would you go about doing this in PyQGIS, so that it is unaffected by QGIS version changes? Specifically, how could you use zonal calculator in PyQGIS to calculate mean of areas in a raster using polygons from a shapefile?


texts: Convert lat/lng to x/y pixel 2D PNG image

texts:
I have a 2D PNG map which I know the corners of my image and I also have the image size all as below:
image size: 600x600 pixel
top-right coordinate: (35.74706129,52.08012633) => (lat,long)
bottom-left coordinate: (35.66569172,52.00198088) => (lat,long)

Also I have a lot of points in lat/long that I need to find them on my map.
How can I convert lat/long to pixel to find the exact right place of my points on the map?

i did it with calculate distance between top and bottom but i got below:


texts: Obtaining yealy maximum NDVI values from Landsat, Earth Engine code editor vs rgee package

texts:
I am trying to obtain the Landsat yearly maximum NDVI values for a region using Earth Engine through the R package rgee (an implementation of Earth Engine Python API for R). So far, I had success in achieving this in the Earth Engine code editor following a similar case in https://mygeoblog.com/2017/05/18/annual-rainfall/ . However, I have not been able to achieve this in the rgee package using a similar code.
This is my code in the code editor:
function addNDVI(image) {
  var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
  return image.addBands(ndvi);
}

var region = ee.Geometry.Polygon([[-105, 39.05], [-104.5, 39.05], [-104.5, 38.75], [-105, 38.75], [-105, 39.05]]);

var landsat = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                  .filterDate('2014-01-01', '2019-12-31')
                  .filterBounds(region)
                  .map(addNDVI)
                  .select('NDVI');

var years = ee.List.sequence(2014, 2019);

var yearlyLandsat = ee.ImageCollection.fromImages(
  years.map(function (year) {
    var annual = landsat
      .filter(ee.Filter.calendarRange(year, year, 'year'))
      .max();
    return annual
      .set('year', year)
      .set('system:time_start', ee.Date.fromYMD(year, 1, 1));
  }));
print(years)
print(yearlyLandsat)

This is my code on rgee so far:
library(rgee)

ee_Initialize()

addNDVI <- function(image) {
  ndvi <- image$normalizedDifference(c('B5', 'B4'))$rename('NDVI');
  return(image$addBands(ndvi))
}

region <- ee$Geometry$Polygon(list(c(-105, 39.05), c(-104.5, 39.05), c(-104.5, 38.75), c(-105, 38.75), c(-105, 39.05)))

landsat <- ee$ImageCollection('LANDSAT/LC08/C01/T1_SR')$
                  filterDate('2014-01-01', '2019-12-31')$
                  filterBounds(region)$
                  map(addNDVI)$
                  select('NDVI')

years <- ee$List$sequence(2014, 2019)

landsatYearly <- ee$ImageCollection$fromImages(
    years$map(function(year) {
        yearly <- landsat$
            filter(ee$Filter$calendarRange(year, year, "Year"))$
            select("NDVI")$max()
        return(yearly$
               set("year", year)$
               set("system:time_start", ee$Date$fromYMD(year, 1, 1)))
}))

This code for rgee keeps giving me this error: Evaluation error: argument "year" is missing, with no default.
I believe the answer to this issue may be related to Map function over ImageCollection in Python 2.7 in Earth Engine . However, I have not been able to make it work. How this part for calculating the yearly maximum NDVI should be written for rgee? The logic should be similar to applying this in the Earth Engine Python API.


texts: Getting field index from input vector layer in PyQGIS?

texts:
I am very new to QGIS3 and I am currently trying to convert a python script from QGIS 2.18 to QGIS 3.
Therefore, I am trying to integrate my python script into the template from underdark -> see here

I am currently stepping from one error into another so I seek help as I was not able to find any other question on Stack Exchange or any documentation.

Here is some part of my script so far:

def processAlgorithm(self, parameters, context, feedback):
    inEdges = self.parameterAsVectorLayer(parameters, self.INPUT_VECTOR_LAYER_EDGES, context)
    inNodes = self.parameterAsVectorLayer(parameters, self.INPUT_VECTOR_LAYER_NODES, context)
    outMerged = self.parameterAsOutputLayer(parameters, self.OUTPUT_VECTOR_LAYER_MERGED, context)

    #DO SOME CALCULATION

    # get field indices for inputlayer
    idxEdgeId = inEdges.fieldNameIndex(ID)
    idxPriority = inEdges.fieldNameIndex(PRIORITY)
    idxLength  = inEdges.fieldNameIndex(LENGTH)
    idxTime    = inEdges.fieldNameIndex(TIME)
    idxFrom  = inEdges.fieldNameIndex(FROM)
    idxTo  = inEdges.fieldNameIndex(TO)

    idx_id = inNodes.fieldNameIndex(NODEID)
    idx_x = inNodes.fieldNameIndex(COORD_X)
    idx_y = inNodes.fieldNameIndex(COORD_Y)

    # Set the provider to accept the data source
    outdp = outMerged.dataProvider()

    # define attributes and add them to the vector layer
    outdp.addAttributes([QgsField(EDGEID, QVariant.String),
                        QgsField(NODE1, QVariant.String),
                        QgsField(NODE2, QVariant.String),
                        QgsField(PRIORITY, QVariant.Int),
                        QgsField(VTO, QVariant.Double),
                        QgsField(VFROM, QVariant.Double)])
    outMerged.updateFields()

    # get field indices for outputlayer
    outEdgeId = outdp.fieldNameIndex(EDGEID)
    outNode1 = outdp.fieldNameIndex(NODE1)
    outNode2 = outdp.fieldNameIndex(NODE2)
    outPriority = outdp.fieldNameIndex(PRIORITY)
    outVTo = outdp.fieldNameIndex(VTO)
    outVFrom = outdp.fieldNameIndex(VFROM)

    listOfEdgeIds = []
    count = 0
    flag = 0

    (...)

For your information: I have 2 input vector layers (one point, and one line layer) and out of that, I want to create one outputlayer (which should be a line-vector layer as well).

But currently I have trouble with the "fieldNameIndex(QString)" - method in the first section of the code. I have already tried to call the method with the dataProvider() of the input layer, but that hasn't worked out. I get error:

  Traceback (most recent call last): File
  "C:\Users\tlichten\AppData\Roaming\QGIS\QGIS3\profiles\default\processing\scripts\merging.py",
  line 89, in processAlgorithm idxEdgeId = inEdges.fieldNameIndex(ID)
  AttributeError: 'NoneType' object has no attribute 'fieldNameIndex'

Although I have to input vector layers, I cannot select any of them if I want to start the script. The execute script window looks like this:

In QGIS2.18, there would be two drop-down menus in which to choose between the input layeres.

I have adapted the initAlgorithm() to the following:

    def initAlgorithm(self, config=None):
    self.addParameter(QgsProcessingParameterVectorLayer(
        self.INPUT_VECTOR_LAYER_EDGES,
        self.tr('Edges'),
        [QgsProcessing.TypeVectorAnyGeometry]))
    self.addParameter(QgsProcessingParameterRasterLayer(
        self.INPUT_VECTOR_LAYER_NODES,
        self.tr('Nodes')
        [QgsProcessing.TypeVectorAnyGeometry]))
    self.addParameter(QgsProcessingParameterRasterDestination(
        self.OUTPUT_VECTOR_LAYER_MERGED,
        self.tr("Output Vector Layer Merged"),
        [QgsProcessing.TypeVectorAnyGeometry]))

But I can still not select any input layers if I execute the script.


texts: Displaying cross sections/3D volumes in QGIS 3.X

texts:
Is there is a way to display vertical images (cross sections) in a 3D view window in QGIS 3? I want to display some cross sections underneath a geological map overlain on the topography.

I can do this in ArcMap & ArcScene, but would like to know if its possible in QGIS. My input data would be a series of images (the vertical profiles) hanging off 3D polylines which are their surface profiles.

Something like this:

Image courtesy of the USGS


texts:
My WMS layer(EPSG:3847) looks fine in GeoServer

However, it is distorted in browser:

I found someone hade the same problem(Openlayers and Geoserver: OSM/Google Maps and WMS Overlay)but I still couldn't solve the problem. I also tried to add code like

var geographic = new OpenLayers.Projection("EPSG:3847");
var mercator = new OpenLayers.Projection("EPSG:4326");
var lund_center = new OpenLayers.LonLat(13.2, 55.7).transform(
                geographic, mercator

and

map.setCente(lund_center, 10);

but the WMS layer is still distorted.
I also tried this code but nothing's changed (EPSG3847 is not included in proj4.js, is that why the code didn't work?):

EPSG4326 = new OpenLayers.Projection("EPSG:4326");
EPSG3847  = new OpenLayers.Projection("EPSG:3847 ");

var option = {
            projection: new OpenLayers.Projection("EPSG:3847 "),
            displayProjection: new OpenLayers.Projection("EPSG:4326")
            };

How can I solve this problem? Any help please?
*the code I'm trying now are as follows, still,nothing changed：

    var lon = 12.1;
    var lat = 48.4;
    var zoom = 10;
    var map, vectors, controls;
 var options = {
        controls: [],
        projection: new OpenLayers.Projection('EPSG:3857'),
        displayProjection: new OpenLayers.Projection('EPSG:4326'),
        maxExtent: new OpenLayers.Bounds(
        -20037508, -20037508,20037508, 20037508.34),
        //modify parameter
        maxResolution: 5000,
        //modify parameter
        //setting zoom levels
        numZoomLevels:15,
        units: 'm',

    };

 function init(){

    var map = new OpenLayers.Map({
    div: "map",
    allOverlays: true
    });

    var roads_wms = new OpenLayers.Layer.WMS( "Roads",
            "http://localhost:8090/geoserver/wms/wsNYCRoad", {
            projection: map.displayProjection,
            layers: 'vagar_alla',
            styles: '',
            srs:'EPSG:3006',
            transparent:true,
            tiled: 'false',
            tilesOrigin : map.maxExtent.left + ',' + map.maxExtent.bottom
        },
        {
            buffer: 0,
            displayOutsideMaxExtent: true,
            isBaseLayer: false,
            reproject: true
        });
        Proj4js.defs["EPSG:3006"] = "+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs";
        Proj4js.defs["EPSG:4326"] = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs";
        var src = new OpenLayers.Projection('EPSG:3006');
        var dst =  new OpenLayers.Projection('EPSG:4326');
        Proj4js.transform(src, dst, roads_wms);
        roads_wms.projection=dst;


texts: Running 64-bit background geoprocessing inside ArcMap?

texts:
I'm running ArcMap 10.2 and just set up 'ArcGIS for Desktop Background Geoprocessing (64 bit)' and now I have both the ArcGIS10.2 folder and ArcGISx6410.2 folder inside my python folder. I have background processing enabled as well. When I open the python window inside of ArcMap and do import sys, print sys.version it shows that it is running the 32 bit version.

How do I get the 64 bit version to run inside of ArcMap?


texts: Checkbox becomes unchecked when map pans

texts:
I have a table of contents that is used to toggle visibility of layers. I also have a dialog box that populates with the visible layers so that the user can pick from the visible layers what layer they would like to identify on. They use a checkbox to do this. However, when the map is panned or a feature is zoomed to the checkbox becomes unchecked.  here is a screen shot.

Here is the code i am using for creating the visible layer checkbox:

function buildLayerList() {
  require(["dojo/on", "dojo/dom", "dojo/_base/array"], function (on, dom, arrayUtils) {
    var mapLayer = map.layerIds;
    var myItems = [];
    arrayUtils.map(mapLayer, function (layerName) {
      var myLayer = map.getLayer(layerName);
      if (myLayer.id !== "NAIP2012" && myLayer.id !== "NAIP2009" && myLayer.id !== "NAIP2006" && myLayer.id !== "trLayer" && myLayer.id !== "Label SEO wells by actual yield" && myLayer.id !== "Label SEO wells by MWBZ top") {

        if (myLayer.visibleLayers) {
          //console.log(myLayer.id + " : " + "has vis Layers");
          var items = arrayUtils.map(myLayer.layerInfos, function (info, index) {
            if (info.name === "Label  SEO wells by depth" || info.name === "Label  SEO wells by actual yield" || info.name === "Label  SEO wells by MWBZ top") {
              // console.log(info.name);
            } else {
              if (myLayer.visibleLayers.indexOf(info.id) > -1) {
                return "<input type='checkbox' class='list_item'" + (info.defaultVisibility ? "checked=checked" : "") + "' id='" + info.name + "CB'' />" + info.name + "
";
              }
            }
          });
          myItems = myItems.concat(items);
        }
      }
    });
    var ll = dom.byId("legendDiv");
    ll.innerHTML = myItems.join(' ');
    document.getElementById("legendDiv").style.display = "block";
  });
}


texts: Creating evenly distributed points within an irregular boundary

texts:
I need to create an evenly distributed series of points within a series of oddly shapped polygons (formerly squares, but now squares with donut holes).

The way I have solved this problem so far is to create a fishnet of the polygon and then use the centroid of each unit that the fishnet creates.

However, the problem has become more complex and I now have more complex polygons. The centroids of the fishnet units are no longer good enough.

I was trying to convert the polygons to a raster and then use the Split tool for rasters, and create an output with a specified number of equal area units, but that won't work, as my input vector data doesn't have the necessary vaues for that raster process to run properly.

I am working with Arc 9.3 (but also have access to several other software packages)


texts: Making shapefile Read-Only to ArcGIS for Desktop?

texts:
How do you make the shapefile a Read-only so the client can not make any edits using 10.2.2 ArcGIS.  It is not within a geodatabase, just a stand alone shapefile.

I am using Windows 8


texts: PostGIS - raster - get multiple values at once

texts:
Is there a way, how to obtain more points from raster in one query faster than this?

SELECT x, y, ST_Value(raw_data, 1, x, y) AS val_0,
      ST_Value(raw_data, 2, x, y) AS val_1
FROM data
CROSS JOIN generate_series(1, 256) As x
CROSS JOIN generate_series(1, 256) As y
WHERE ((x = 184 AND y = 192) OR (x = 185 AND y = 192) OR ....)

Speed of this is the same as if I call ST_Value for every pair manually (but that is not, what I want)

Also, I have two bands and need value from both. But ifthere is something faster for only one band data, I can separate them.


texts: What are the necessary correction/calibration on Landsat 8 imagery for land cover classification?

texts:
I downloaded calibrated/corrected Landsat 8 data (Level 2) from ESPA (https://espa.cr.usgs.gov/) to do land cover classification. But there are many corrections and calibrations can be done as listed below and maybe more, such as:

 - Surface Reflectance
 - Top of Atmospheric
 - sensor zenith and azimuth
 - solar zenith and azimuth
 - geometric
 - radiometric

I am not sure all of them were processed on the scene. What are the must do corrections/calibrations for Land cover classification analysis?


texts:
After running a SAGA line dissolve process in QGIS, an MSHP file is created (amongst other things). Fine by me, though I'm not sure what it's used for (metadata in XML format?). Problem is that it can't be deleted without closing QGIS. This is a problem, as I'm running the process repeatedly in a Python script.

Here is the relevant code:

lines_layer = some_vector_polygon_layer
dissolved_lines_path = "SOME SHAPEFILE PATH HERE"
unused_field = lines_layer.dataProvider().fields()[0].name() # Name of a field, not used
processing.runalg("saga:linedissolve", lines_layer, unused_field, unused_field, unused_field, 1, dissolved_lines_path)

I then load the resulting layer and add it to the map registry & table of contents. Shortly after, I delete the layer from the map registry (& table of contents).

I can delete the resulting SHP quite fine, but I am unable to delete the MSHP. There are no errors in the resulting shapefile or in the processing log. Any ideas?


texts: How to add a bottom bar for citations in QGIS Print Composer?

texts:
I want to add something like the bottom bar in this picture in order to include my citations.

Is there a way to do so?


texts: Changing python code, to only import points from dxfs (or to only import lines)?

texts:
edit -  This code no longer works in the new version of QGIS (QGIS 3.0)
If I paste it into the python console it generates a series of requests for the CRS but then imports nothing.
QGIS 3.0 has moved to a newer version of Python and the code I have is incompatible.

How do I edit it to make it work again?

Original question:

I just merged several hundred layers into QGIS which had been imported using this code:-

import glob, os
path = "H:/aOsmap/Miltkeyn/2012 ntfs/all/"
for layer in glob.glob(path + "*.dxf"):
    vlayer = QgsVectorLayer(layer, 'name', 'ogr')
    subLayers = vlayer.dataProvider().subLayers()
    for subLayer in subLayers:
        geom_type = subLayer.split(':')[-1]
        uri = "%s|layername=entities|geometrytype=%s" % (layer, geom_type,)
        dfx_file_name = os.path.splitext(os.path.basename(layer))[0]
        layer_name = "%s - %s" % (dfx_file_name,geom_type,)
        sub_vlayer = QgsVectorLayer(uri, layer_name, 'ogr')
        QgsMapLayerRegistry.instance().addMapLayer(sub_vlayer)

Merging them was difficult because the merge window (vector-merge layers) mixes up the linestring layers and the point layers so selecting just the point layers or just the linestring layers, from about over 600 layers, is laborious and prone to clicking on the wrong layer.

It would be easier and quicker to import all the points from my dxfs in one go and merge those (saving the merged file) and then import all the linestrings from the dxfs and merge those.

However, I don't understand the code so I have no idea what to alter to just select points or just select lines.

Can anyone help me?


texts: Calculate area of polygons using OGR in python script

texts:
I have a shapefile of many thousands of polygons. I'm trying to append an area field to the attribute table and calculate the area of each polygon (in sq. meters)

driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(path_to_shp_data, 1)
layer = dataSource.GetLayer()
new_field = ogr.FieldDefn("Area", ogr.OFTReal)
new_field.SetWidth(32)
layer.CreateField(new_field)

for feature in layer:
    geom = feature.GetGeometryRef()
    area = geom.GetArea()
    feature.SetField("Area", area)
    layer.SetFeature(feature)

dataSource = None

The code runs successfully but the resulting "Area" field contains all 0's.

Projection info is as follows:

PROJCS["WGS_1984_UTM_Zone_48N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",105],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]


texts:
I am using following code to take subset of all files in a folder. The code is not giving any error but it is not creating output files. How to solve it?

import numpy as np
import re
import gdal
import ogr
from osgeo import gdal_array
from gdalconst import *
import os
import osr
import subprocess
e=1
indir = '/home/vigna/Documents/polynyas/SIC/2017/Sept'
directory=os.fsencode(indir)
for file in os.listdir(directory):
    file1=os.fsdecode(file)
    if file1.endswith(".tif"):
        ds = gdal.Open(os.path.join(indir,file1),GA_ReadOnly)
        out=str(e)+'_subset.tif'
        inputf = file1# input raster
        #print(len(inDS))
        inDS=str(inputf[0:33])
        #print(inDS)
        outDS = out# output raster
        #print(outDS)
        ulx = -161255
        uly = 2971737
        lrx = 437613
        lry = 2475328
        translate = 'gdal_translate -projwin %s %s %s %s %s %s' %(ulx, uly, lrx, lry, inDS, outDS)
        #print(translate)
        os.system(translate)
        e+=1


texts: Adding attribute features from the intersection of a grid with a point shapefile?

texts:
I want to calculate the frequency of observation of an animal in a forest habitat within the different zones of my study area. To do this I will use the simple equation (number of individuals)/(survey effort). I want to do this within a grid of 150mx150m, where a greyscale gradient will indicate the varying population densities.

So far, I have a polygon grid (fishnet) clipped to an area of study of 150mx150m. The survey effort per grid has been calculated using the 'sum line' function. (vector->analysis tools->sum line length)

Picture: a visual representation of the field site (green polygon), with survey trails (pink lines) and observation locations (dark red points).

Next, I want to calculate the number of individuals observed per grid square. The point shapefile representing the observation locations contains a feature in the attributes table called 'cluster size' which is the number of animals that were observed for any observation point. I want to have the total addition of 'cluster size' for all points within each grid cell. Preferably as a feature (column) in the 'sum line' shapefile so that I can then make a simple calculation using the equation (number of individuals)/(survey effort) for population density within each cell of the grid.

I have tried using the point count function found in vector->analysis but this provides only the number of points per grid cell, while I would want the sum of cluster sizes per grid cell. Can anyone suggest how to do this?

I am using QGIS 2.18.7


texts: Removing small "salt & pepper" polygons from layer in QGIS?

texts:
I have a map made up of multiple polygons, with a few hundred thousand points very small polygons on within it.

I need the larger polygons, but not the points small polygons.  How do I select and remove all the points small polygons, while preserving the larger polygons?  (I obviously do not want to do this manually - there are a few too many points for that!)

I'm totally new to QGIS, and I've Googled extensivey with no luck.  Probably because I'm unsure what the correct terminology to use in a Google search is!

The files I am using are a .dbf, .prj, .qpj, .shp and .shx.  I'm simply dragging the .shp file on to QGIS, and the map appears as one layer.  Basically a map of the UK divided into regions (which I want) and covered in points (which do not want).

This is what it looks like up close.  I did a dissolve to see if they'd go in to the surrounding layer, with no luck:


texts: Inserting spatial data into existing PostGIS geometry column from shapefile based on ID

texts:
I have an existing PostgreSQL database table "city_data" with Non-spatial data in the following format:

City_ID as varchar, City_name as varchar, pop as int,....

I have enabled the PostGIS extension to this existing database and added a geometry column "the_geom" to "city_data" table with SRID and other essential details same as my shapefile.

Now I have a shapefile with city boundary as polygon and it also has city_id which is exactly same as PostgreSQL table.

I would like to load the polygon data from shapefile and insert/update into corresponding rows in PostGIS table based on city_id.

Is it possible to do directly or is there some workaround for doing this task?


texts: Saving QML styles directly into QGS project file

texts:
I have a QGIS Model (saved in directly in QGIS 3 project file ".qgs") which creates some layers and assign predefined styles to them. I would like to save those QML styles into my .qgs and reference the styles there directly in my QGIS model.

Is it possible to store QML Style files directly in .qgs project? Or is there other solution to this? The new version of QGIS Project now allows to store a lot of kind of data diretcly in it.


texts: Select shapes in one layer that overlap with shapes from another layer?

texts:
This question is a duplicate of this: How to select only polygons that overlap with another layer?

But in my question I want to know how to do this in QGIS not in ArcGIS.
I want to be able to select shapes in one layer if they intersect /overlap shapes in another layer.

Intersection does not work because it gives me these odd silver polygons:

Because both the shapefiles are not perfectly aligned, it gives those small odd shapes, I would like to instead select entire shapes and not intersect the layers.


texts: i want to Search for points within a raster layer using ST_Contains. Keep getting no function matches given error

texts:
I have made a Simple points shapefile with seven points and I would like to be able to query which of these points are contained within the boundary of the raster layer.

SELECT "points inside raster_test".name
FROM "points inside raster_test", myraster3
where ST_Within ("points inside raster_test".geom, myraster3.rast);

I keep getting the following error-

ERROR:  function st_contains(geometry, raster) does not exist
LINE 3: where ST_contains ("points inside raster_test".geom, myraste...
              ^ HINT:  No function matches the given name and argument types.
You might need to add explicit type casts.

Do I have to write this query differently because there is a raster instead of another vector layer? I have tried using a vector polygon with the same points and wrote it exactly the same and it has given me the results I have desired.


texts: Is it possible to have a Leaflet control with large images / icons

texts:
Im trying to create a map where the layer control have large icons or images, and the "buttons" are located in the bottom of the map and not on the sides. I have used the simple controls until now.

I want something similar to this, where the control is located in along the bottom and not in the side.

https://leaflet-extras.github.io/leaflet-providers/preview/

Does anyone have a good example?


texts: inaccurate distance measurements in Python

texts:
I've tried several ways of measuring the distance between two points in my Django application and comparing the results to a reliable measurement. My numbers are way off. For example, I'm assuming the real distance between Dublin and Liverpool is ~217 km as reported by Google Maps:

Using geopy (should be extremely accurate):

>>> from geopy.distance import distance
>>> dublin = (-6.270447, 53.339791)
>>> liverpool = (-2.991028, 53.402061)
>>> distance(dublin, liverpool).km
362.70989455939394

Using Django's GEOS API (less accurate linear calculation):

>>> from django.contrib.gis.geos import Point
>>> dublin = Point(-6.270447, 53.339791, srid=3857)
>>> liverpool = Point(-2.991028, 53.402061, srid=3857)
>>> dublin.distance(liverpool)*100
328.00101418228877

EDIT:
Using a better projection for this area (UTM 30N) yields almost the same result:

>>> dublin.transform(32630)
>>> liverpool.transform(32630)
>>> dublin.distance(liverpool)*100
328.32200116442266

In both cases I'm off by over 100km! Measuring small distances (<1km) is just as inaccurate. What am I doing wrong here!?


texts: QGIS print composer, export as image, rest ESRI raster as basemap

texts:
I am trying to export an image from the print composer using this layer [ESRI rest layer][1] as a basemap. http://services.arcgisonline.com/arcgis/rest/services/ESRI_Imagery_World_2D/MapServer

I connected to the layer using the python plugin and this tutorial. http://hub.qgis.org/wiki/quantum-gis/Arcgis_rest

When I go to export my map made in print composer as an image the basemap doesn't show up and I get a grey hatching with 'no map data' on it. It worked a few times and then quit for some reason. The print composer version looks great but the only way I can export seems to either screen grab or save the maps as images and then add through the print composer.

Any thoughts on how to solve this issue? It worked so I imagine i've twiddled something and need to untwiddle it to resolve the issue.


texts: Raster statistics - incorrect and apparently arbitrary minimum value

texts:
I have projected a raster with gdal as part of an R workflow:

gdalwarp -s_srs EPSG:4283 -t_srs EPSG:3577 -r bilinear -multi -tr 1000 1000 -dstnodata -9999 in.tif out.tif

When I add out.tif to ArcMap, the minimum value indicated beside the colour ramp is incorrect (268.048, while the true minimum is 91.2931), and the colour ramp applied to the data is stretched between this incorrect minimum and the (correct) maximum.

I can calculate statistics (data management > raster > raster properties toolbox) to remedy this, but I'm concerned that if I fail to do this for one of the many rasters I'm working with, I (or my colleagues) might misinterpret a map.

The trick of switching to Classified symbology and back to Stretched returns the error: This dataset does not have a valid histogram required for classification method.... Attempting the same with Unique Values returns the error: The number of unique values exceeds the default limitation (>65536)....

Why has ArcMap chosen this seemingly arbitrary minimum of 268.048? Is there anything I can do to prevent it from happening (e.g. to have ArcMap either calculate the correct minimum when the raster is added, or to use default statistics that span the possible range of values given the pixel type and depth)?

NB: I am using a fresh installation of ArcGIS 10.2, and "Calculate statistics" is checked in Geoprocessing > Environments > Raster Storage.


texts:
I'm trying to generate evenly-spaced points (1m spacing) along a line shapefile in QGIS 3.4.3. This works with certain shapefiles. However, with other shapefiles, the units automatically change from "meter" to "degrees" and yellow triangles appear (see image below).

Does anyone know why this happens, and how I can make it work?


texts: How to geospatial cluster FeatureCollection of polygons in GEE?

texts:
How do you geospatial cluster a FeatureCollection of polygons in Google Earth Engine? For example, the FeatureCollection fc consists of polygons.
var fc = ee.FeatureCollection('GLIMS/current')
  .filter(ee.Filter.eq('geog_area', 'French Alps'))

Map.addLayer(fc, {},'unclustered', false)
Map.centerObject(fc, 7)

In other words, I want each polygons to have an additional property 'cluster' which can be used to group geospatially close polygons.


texts: How school catchment area analysis is done

texts:
In fact i am asking for the method and the data needed for determining school catchment area.Example is http://www.qgso.qld.gov.au/maps/edmap/.

Probable Prerequisites:

Catchment boundary must not cross buildings and should follow transport network.
There not be any un-allocated area (void) i.e. area is not under any  catchment.
I need some sort of automation at most extent.

Overview-- see below

When zoomed in a bit more-- see below

N.B. I think Example website related one can give elaborate answer.


texts: Associating polygons with overlapping polygons by area of overlap (known as allocation factor)

texts:
I am attempting to associate School Attendance Boundary (SAB) polygons with Zip Code Tabulation Area (ZCTA) polygons. They are in separate .shp files.
One approach is to weight the allocation of a SAB to a ZCTA by the proportion of the SAB area that lays in the ZCTA area, i.e. if SAB_1 is 100ft^2 and 10ft^2 lay in ZCTA_1, then the allocation factor to be used in associating measures from SAB_1 to measures from ZCTA_1 is 0.1.
I note that these are common values for associating geographic areas (see MCDC Geocorr for examples). A previous Q&A (Calculating percentage of overlap between two different layers) provided a helpful walkthrough, but was not sufficiently detailed for me.
For example, one instruction was to "execute a union" and set a target layer. To my knowledge, unions don't involve targets.
How can I obtain the allocation factors I have described?


texts: Data from REST service shows symbology for items not on my map

texts:
I am using QGIS 3.4 to display eelgrass data from this website:
https://geo.wa.gov/datasets/8ef27802b6cf4b61b814c80dd91d6e9e.
To get the specific layer I need I use this link to connect to an ArcGIS mapping layer and add layer #42:
https://gis.dnr.wa.gov/site3/rest/services/Public_Aquatics/WADNR_PUBLIC_AQR_NS_HBTT/MapServer
This adds the data to my QGIS project but the data are classified into two classes called continuous and patchy.  I would like to remove continuous from the symbology and legend as this class does not appear in my area of interest.  I cannot figure out how to alter the text or remove continuous from the legend in QGIS using the layers symbology.  I have tried to download the data as a TIFF file for my area of interest and this creates a TIFF file of junk with only 2.5KiB of space.
How do only display the patchy class of data in QGIS?


texts:
I'd like to save the polyline coordinates. I'm creating it with Leaflet.Draw. The problem is that the coordinates are stored twice.

49.907018,16.324997,49.907018,16.324997,49.818936,16.340103,49.818936,16.340103....

this is my code

map.on('draw:created', function(p) {
   map.addLayer(p.layer);
    var points = p.layer.getLatLngs();
  puncte1=points.join(',');
  puncte1=puncte1.toString();
  puncte1=points.join(',').match(/([\d\.]+)/g).join(',');
           document.getElementById('coordinates').value=puncte1 ;

});


texts: GeoPandas to_file() saves GeoDataFrame without coordinate system

texts:
I have

GeoPandas 0.2.1 py27_0
fiona 1.7.0 np110py27_1
Python 2.7.10

installed in Anaconda 2-4.1.1-Windows-x86 environment. I am able to construct GeoDataFrame by reading input datasets and operate the data yet saving an output dataset does not preserve the coordinate system.

import geopandas as gpd
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.to_file(driver='ESRI Shapefile',filename=r'C:\GIS\Temp\world_out.shp')

The world_out.shp does contain the features and the schema of the source shapefile, but coordinate system is undefined (<Unknown>). The world_out.prj file is 0 KB and doesn't contain anything.

out = gpd.read_file(r'C:\GIS\Temp\world_out.shp')
print out.crs
>>> {}

I have tested to perform this operation using pure Fiona and it works fine, preserving the crs information:

import fiona

datasrc_path = r'C:\Program Files (x86)\ArcGIS\Desktop10.4\TemplateData\TemplateData.gdb'

with fiona.drivers():
    with fiona.open(datasrc_path,layer='cities') as source:
        meta = source.meta
        meta['driver'] = 'ESRI Shapefile'
        meta['layer'] = source.name

        with fiona.open(r'C:\GIS\CitiesFiona.shp', 'w', **meta) as sink:
            for f in source:
                sink.write(f)

The CitiesFiona.prj file does contain the information about the coordinate system.

So, it seems as GeoPandas is not able to read/write crs of the data frame into the output dataset when calling the to_file() method. I have tried calling to_file() on GeoDataFrames constructed from various sources and datasets and the coordinate system information was never present in the output dataset.

Has anyone run into this issue before? Can anyone test this?


texts: Using Split By Attributes in ArcGIS ModelBuilder?

texts:
I wish to use "split by attribute" tool in ArcGIS ModelBuilder to split a shapefile (points with metadata) into different groups based on a parameter. The tool works and creates for example 4 .dpf files in the given folder but it does not seem to want to add these new files (tables) into ArcGIS automatically even though I specify "add to display", but more importantly the output of the tool in ModelBuilder seems to be a dead-end that I cannot run any analysis on. This is a bottleneck in my model as I want my model to carry on working on these data sets.

How can I do that?


texts: How to add top and bottom margins to leaflet popup with scroll

texts:
I've been customizing Leaflet in order to let popup display more content by adding scrolling to them. Everything works fine except my content touches top and bottom of the popup window.

How can I add a margin at top and bottom ?

Here's my custom CSS:

.custom-popup, .leaflet-popup-content-wrapper {
    background:#e7e7e7;
    color:#504e4e;
    font-family: 'Molengo', sans-serif;
    font-size:12px;
    line-height:10px;
    border-radius: 2px;
    height : 10 px ;
    max-height:300px;
    overflow:auto;
    }

.leaflet-popup {
    position: absolute;
    text-align: center;
    }

.custom-popup, .leaflet-popup-tip {
    background: #e7e7e7;
    border: none;
    box-shadow: none;
    }

.leaflet-popup-content {
     min-width: 100 px !important;
    }

And here's a capture of what I have and what I would like to get:

Scrollbar is on the right side of the popup :

EDIT : Blur effect I try to obtain (red arrows)


texts: Joining CSV-Generated Datapoints with Shapefiles in ArcGIS Desktop?

texts:
I have a question about joining data points (with a large amount of information) with shapefiles in ArcGIS Desktop:

I have imported a CSV file into GIS. This CSV File is quite large: each row contains roughly 75 columns. As such, each row in the CSV file contains quite a bit of information. I was able to create an "Events" layer with the CSV document and display the information visually as a series of data points in GIS (each row of the CSV document has latitude and longitude coordinates to make this possible).
It is my goal to spatially-join the data points with the polygon shapefile that that they fall inside of (please see the picture below). For a bit of reference the polygon is an island and the data point describes information regarding that island. Unfortunately GIS does not allow for this to happen because there is no OBJECT_ID attached to the data points. Additionally, exporting the data to create a new layer does not work because the amount of information is too large!

Any thoughts on how to overcome this issue?

I am quite new to GIS.


texts: How to avoid multiple label in zoomin in tilemill

texts:
Now i am working on tilemill ,and i am facing a problem in a tilemill when i am styling zoomin a road layer that time display road layer label multiple times.
How can avoid disply multiple road layer labels in tilemill.


texts: How to hide label when"count" is equal to 1 in pointstacker?

texts:
Is there a way to hide labels when the "count" value is 1 (one), using pointstacker?

"count" = <ogc:PropertyName>count</ogc:PropertyName>

SLD Code:

<?xml version="1.0" encoding="utf-8"?>
<StyledLayerDescriptor version="1.0.0"
 xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
 xmlns="http://www.opengis.net/sld"
 xmlns:ogc="http://www.opengis.net/ogc"
 xmlns:xlink="http://www.w3.org/1999/xlink"
 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <NamedLayer>
    <Name>Pontos de Ônibus</Name>
    <UserStyle>
      <Title>Pontos Empilhados</Title>
      <FeatureTypeStyle>
        <Transformation>
          <ogc:Function name="vec:PointStacker">
            <ogc:Function name="parameter">
              <ogc:Literal>data</ogc:Literal>
            </ogc:Function>
            <ogc:Function name="parameter">
              <ogc:Literal>cellSize</ogc:Literal>
              <ogc:Literal>100</ogc:Literal>
            </ogc:Function>
            <ogc:Function name="parameter">
              <ogc:Literal>outputBBOX</ogc:Literal>
              <ogc:Function name="env">
                <ogc:Literal>wms_bbox</ogc:Literal>
              </ogc:Function>
            </ogc:Function>
            <ogc:Function name="parameter">
              <ogc:Literal>outputWidth</ogc:Literal>
              <ogc:Function name="env">
                <ogc:Literal>wms_width</ogc:Literal>
              </ogc:Function>
            </ogc:Function>
            <ogc:Function name="parameter">
              <ogc:Literal>outputHeight</ogc:Literal>
              <ogc:Function name="env">
                <ogc:Literal>wms_height</ogc:Literal>
              </ogc:Function>
            </ogc:Function>
          </ogc:Function>
        </Transformation>
        <Rule>
          <Name>Pontos de Ônibus</Name>
          <Title>Ícone SVG</Title>
          <PointSymbolizer>
            <Graphic>
             <ExternalGraphic>
              <OnlineResource xlink:type="simple" xlink:href="v1.1_pontos.svg" />
              <Format>image/svg</Format>
             </ExternalGraphic>
             <Size>
              <ogc:Literal>30</ogc:Literal>
             </Size>
           </Graphic>
          </PointSymbolizer>
          <TextSymbolizer>
            <Label>
             <ogc:PropertyName>count</ogc:PropertyName>
            </Label>
            <Font>
              <CssParameter name="font-family">Arial</CssParameter>
              <CssParameter name="font-size">12</CssParameter>
              <CssParameter name="font-weight">bold</CssParameter>
            </Font>
            <LabelPlacement>
              <PointPlacement>
                <AnchorPoint>
                  <AnchorPointX>0.5</AnchorPointX>
                  <AnchorPointY>0.8</AnchorPointY>
                </AnchorPoint>
              </PointPlacement>
            </LabelPlacement>
            <Stroke>
              <CssParameter name="stroke">#ffffff</CssParameter>
            </Stroke>
            <Halo>
               <Radius>1.5</Radius>
               <Fill>
                 <CssParameter name="fill">#000000</CssParameter>
               </Fill>
            </Halo>
             <Fill>
               <CssParameter name="fill">#ffffff</CssParameter>
             </Fill>
          </TextSymbolizer>
         </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>


texts: I am using the latest installation of QGIS (3.2) in a Macos Sierra 10.12.6 and I get those graphic issues on the top left of the software (see image). Also when I go to full screen the tool tabs get duplicated (second image). Does anyone else has the same issue? Is it related with my graphics card? I have tried to toggle off/on the toolbars but this doesn't fix the issue.


texts: Finding difference in 2 shapefiles and deleting the difference

texts:
I have 2 shapefiles: shapefile A (yellow) and shapefile B (purple). I already did some basic data cleanup on shapefile A, but now want to use shapefile B because it has less overlaps between polygons.

How can I detect the difference between the 2 shapefiles? Basically I want to remove the differences. So in the following image, the purple polygons on the right (shapefile B) should not be there as they were already removed in shapefile A.

I tried in the Difference tool in QGIS but it didn't help.


texts: Plugin "processing" does not work - not possible to enable - QGIS 3.4.3 Madeira - Windows 10 Pro

texts:
The "Spatial data processing framework for QGIS" does not work anymore.

Also known as "processing".... this tool can not be uninstalled...

and it is disabled.

QGIS 3.4.3 Madeira - Windows 10 Pro

When i try to enable it, QGIS crashes...

I tried to uninstall QGIS and also cleaned the registry and install QGIS again, but the same happens again...it seems that QGIS saves all the adjustments...


texts: Is there a style sheet available for OS Open Roads layer?

texts:
I've added this vector layer into QGIS but unlike other layers available from the OS I haven't been able to find a standard style to work with it. Is there anything available out there?

https://www.ordnancesurvey.co.uk/business-and-government/products/os-open-roads.html


texts: Python script with forward slash in list causing errors?

texts:
Here is my script:

import arcpy

# set environment, workspace, and output location

arcpy.env.overwriteOutput = 1

arcpy.env.workspace = r"Database Connections\Transportation.sde"

output_loc = r"C:\Users\johndiaz\Documents\Outerspatial_Shapefiles"

agency_list = ['CASP', 'CIC', 'MRCA/SMMC', 'LACODPR', 'NPS', 'PCTA']

# assign variables to layers

th_fc = r"Database Connections\Transportation.sde\eGIS_Transportation.EGIS.DPR_TRAILHEADS"
poi_fc = r"Database Connections\Transportation.sde\eGIS_Transportation.EGIS.DPR_TRAILS_POI"
trails_fc = r"Database Connections\Transportation.sde\eGIS_Transportation.EGIS.DPR_TRAILS"

# shapefiles for trailheads and poi
for list in agency_list:
    print "generating trailhead and poi shapefiles"
    where_clause = "agency = '" + list + "'"
    arcpy.MakeFeatureLayer_management (th_fc, "th_lyr")
    arcpy.SelectLayerByAttribute_management ("th_lyr", "NEW_SELECTION", where_clause)
    output_th = "\\" + list + "_TRAILHEADS.shp"
    arcpy.FeatureClassToFeatureClass_conversion ("th_lyr", output_loc, output_th)
    arcpy.MakeFeatureLayer_management (poi_fc, "poi_lyr")
    arcpy.SelectLayerByAttribute_management ("poi_lyr", "NEW_SELECTION", where_clause)
    output_poi = "\\" + list + "_POI.shp"
    arcpy.FeatureClassToFeatureClass_conversion ("poi_lyr", output_loc, output_poi)
    print "trailhead and poi shapefiles complete"

for list2 in agency_list:
    print "generating trails shapefiles"
    where_clause = "trail_type = 'Existing' AND agency = '" + list2 + "'"
    arcpy.MakeFeatureLayer_management (trails_fc, "trails_lyr")
    arcpy.SelectLayerByAttribute_management ("trails_lyr", "NEW_SELECTION", where_clause)
    output_trails = "\\" + list + "_TRAILS.shp"
    arcpy.FeatureClassToFeatureClass_conversion ("trails_lyr", output_loc, output_trails)
    print "trails shapefiles complete"

print "SCRIPT COMPLETE"

If i remove the list item 'MRCA/SMMC' the script runs fine.  Otherwise, I get this error:

Message File Name   Line    Position
Traceback
    <module>    C:\Users\johndiaz\Documents\trails_app_create_shapefiles.py 36
    FeatureClassToFeatureClass  c:\program files (x86)\arcgis\desktop10.5\ArcPy\arcpy\conversion.py 1891
ExecuteError: ERROR 000210: Cannot create output C:\Users\johndiaz\Documents\Outerspatial_Shapefiles\\MRCA/SMMC_TRAILHEADS.shp
Failed to execute (FeatureClassToFeatureClass).

Any advice?


texts: How to save and load elements to/from a layer in correct z-order?

texts:
I copy elements from one layer to another like this:

IGraphicsContainer container = _mapControl.Map as IGraphicsContainer;
IGraphicsContainerSelect containerSelect = container as IGraphicsContainerSelect;

containerSelect.SelectAllElements();
IEnumElement enumElement = containerSelect.SelectedElements;
enumElement.Reset();

IElement element = enumElement.Next();
while (element != null)
{
  int zorder = 0; // How to get this?
  layerContainer.AddElement(element, zorder);
  element = enumElement.Next();
}

This works fine so far but I obviously like to preserve the element's zorder. I tried to get the zorder with graphicsContainer.GetElementOrder(enumElements) which returns an object (COM object). Well, but what to do with that one, eh?

How can I get the zorder of an element?

Update
The enumElement is already in the correct order (Thanks, Kirk). I now took a deeper look into it.
The elements get saved in the correct z-order, but when I load the layer from file (with IMapControl.AddLayerFromFile) they appear in reversed z-order every time.
I tried to save them with reversed z-order but on the loaded layer they still appear in the wrong z-order.

How can I load the layer with the elements in correct z-order?


texts:
I'm rendering  a MVT (Mapbox Vector Tiles) layer on my Leaflet map. I'm using custom icons to represent my points, so my vector tile options look like this:

var vectorTileOptions = {
maxNativeZoom: 14,
vectorTileLayerStyles: {
    'vw_meteo_fc': function(properties, zoom) {
        console.log(properties.label);
        return {
            icon: new L.Icon({
                            iconUrl: 'icons/icon.png',
                            iconSize: [40, 40], // size of the icon
                            iconAnchor: [20, 20], // point of the icon which will correspond to marker's location
                        })
        }
    }
},
interactive: true,  // Make sure that this VectorGrid fires mouse/pointer events
getFeatureId: function(f) {
    return f.properties.osm_id;
}
};

Is there a way to add labels to these points based on feature subproperty? It can be accessed while rendering the point with: properties.label (console logged in code above).


texts: I have many buildings in CityGML format. Currently cesium does not parse citygml. I tried Geojson and KML formats as in Cesium official tutorials using Sandcastle. In geojson buildings can be seen only 2d (as footprints) and in kml experience Cesium gives an error  "Unsupported geometry type: model". I want to render both 3d geometries and semantic attributes in Cesium. I have approximately 2000 buildings. So, What is the best and easiest way to visualize buildings and their attributes in Cesium?


texts: List all MXDs that a feature class is in using arcpy

texts:
I know it is possible to list the paths of feature classes that are in an MXD using arcpy. But is it possible to do the reverse -- to list all the MXDs that a feature class is in?

I am trying to organize very full file geodatabases that have >1000 feature classes in them. My client wants to keep those feature classes that are used in the most frequently used MXDs. To do this, I was thinking that maybe it was possible to run through every feature class and get a list of every MXD it is in and go from there. Can I do this?


texts: Prevent OGR saving OSM points as features when exporting to GeoPackage

texts:
I'm attempting to export data from the points layer imported from an Open Street Maps pbf file to a GeoPackage file.

In both the script I write to automate the export as well as when I export using QGIS, the data_type column in the gpkg_contents table is set to 'features'.

This is unexpected because when I usually do any export to GeoPackage using any of the other layers, this data_type column is set to one of the geometry types defined by the GeoPackage standard.

My script to export the data to a new DB:

labelDB = QgsVectorLayer("Point", "labels", "memory")

labelDB.startEditing()
labelDB.addAttribute(QgsField("label", QVariant.String))
labelDB.commitChanges()

pointsLayer = QgsVectorLayer("%s/%s|layername=%s" % (input, file, "points"), "points", "ogr")

if pointsLayer.isValid():
    pointsLayer.setSubsetString("place = 'country'")

    #move the resulting features to the labels db
    for feature in pointsLayer.getFeatures():
        labelFeature = QgsFeature()

        labelFeature.setFields(labelDB.fields())
        labelFeature.setGeometry(feature.geometry())

        #get the english name for the country
        name = ""
        otherTags = feature["other_tags"].split(",")

        for tag in otherTags:
            kvp = tag.split('=>')

            if kvp[0].replace("\"", "") == "name:en":
                name = kvp[1].replace("\"", "")

            labelFeature.setAttribute("label", name)

            labelDB.startEditing()
            labelDB.addFeature(labelFeature)
            labelDB.commitChanges()

    QgsVectorFileWriter.writeAsVectorFormat(labelDB, labelsOutDir + "/labels", "utf-8", labelDB.crs(), "GPKG")

How would I go about ensuring that the data is stored as Points instead of features?


texts:
I have a file, Council_Districts.topojson, with the following structure:

{ type: 'Topology', arcs: [...], transform: {...}, objects: {
  Council_Districts_2016: {
    type: 'GeometryCollection',
    geometries: [
      {
        arcs: [...],
        type: '...',
        properties: {
          'DISTRICT': 0,
          'ID': 0,
          'OBJECTID': 0,
        }
      },
      ...
    ]
  }
}

I also have a `Sample_Data.csv' file with the following structure:

DISTRICT, varA, varB, ....,
0, .1, .2, ...

I'm looking to combine these two files such that my output file has the following structure:

{ type: 'Topology', arcs: [...], transform: {...}, objects: {
  Council_Districts_2016: {
    type: 'GeometryCollection',
    geometries: [
      {
        arcs: [...],
        type: '...',
        properties: {
          'DISTRICT': 0,
          'ID': 0,
          'OBJECTID': 0,
          'varA': .1,
          'varB': .2,
        }
      },
      ...
    ]
  }
}

Googling around, I'm finding the previous answer has been to use the external properties flag on the topojson cli, e.g. topojson -e blah, but this flag does not seem to be available within the current topojson 2.0


texts: Alter Field & Add Portion of Filename to Field Name

texts:
I have 350 tables in a geodatabase.
Using a python script I want to iterate through all 350 tables and alter the field names of 3 fields in each table.

I'd like the new field names to contain part of the table's file name. For example, for a table named, chirps_1982_02_stat, I could separate out parts of the file name so that the new field names only had the year and month (1982_02) added to them.
The desired field names would change from:

SUM_grid_code -> Sum_precip_1982_02,
MEAN_grid_code -> mean_precip_1982_02
STD_grid_code -> std_precip_1982_02.

I'm very new to python so I've found instructions on altering a field name but I'm not sure how to parse out the elements of the file's name and add it to the new field name. I am also not sure how to format the code so that I can do this process 3 times, 1 time for each field name I want to change

What I've found for altering field name

>>> arcpy.env.workspace = r'D:\myfolder\mygdb.gdb'
>>> fcList = arcpy.ListFeatureClasses() #get a list of feature classes
... for fc in fcList:  #loop through feature classes
...   fieldList = arcpy.ListFields(fc)  #get a list of fields for each feature class
     for field in fieldList: #loop through each field
...         if field.name.lower() == 'MEAN_grid_code':  #look for the name MEAN_grid_code
...            **arcpy.AlterField_management(fc, field.name, 'Mean Precip %fileyear%', 'Mean Precip XXXX_XX')**


texts: I've created a perl script to create multiple line files from points. It iterates through, creates everything as expected. When you hover over a line, I get a balloon with all of the attribute information that I'm pushing.  Only issue, I don't want 300 separate kml files. So, I modified the script to write to a kml "folder". All of the lines display as expected. Looks just like the individual files, but I can open the one file and see everything. EXCEPT, now the attribute information does not show up when I hover over a line/click on a line/even when I go to properties. It's the same exact  setup that is in the individual file...Is it impossible to show several balloons in one file, even if they are separate "documents"?


texts: QGIS algorithm to copy vector layer with only selected fields

texts:
I am trying to find out if there is a way - or even an algorithm in QGIS - that enables you to effectively create a new feature layer that is a copy of an original layer but with only a selection of key fields that one might want to limit the layer to, for further analysis.
In other words, my original layer has, for example, 30 columns or fields but I want to make a copy of that layer that only has 6 of those fields retained within it. I am working with QGIS 3.14


texts:
I have satellite swath data in the form of an HDF5 file containing 3 arrays: 'data_array', 'latitude' and 'longitude'.
All three are 32-bit floats, are arrays of the same dimensions (X-by-Y) and the names are self-explanatory. All three datasets use a fill value of -32767

What I want to do it reproject data_array from native satellite swath onto a latitude/longitude grid with a pixel size of 0.009 degrees. I want to save the resulting data array as a geoTIFF containing project information (for use in QGis/ENVI/etc).

I have tried using gdal based upon this but gdalwarp segfaults (maybe because there are some fill value pixels in the lat/lon grid? Unsure how to deal with this).

I also tried using python and the pytroll libraries, but find the documentation somewhat impenetrable. All my attempts thus far have failed.

So, can anyone recommend a simple method to perform this reprojection?

I assumed it would be something quite common, but am having difficulty finding a solution.

I can plot project the data onto a grid using NASA Panoply, but this appears to be for visualisation only. I can't see how to save a GEOTiff output.


texts: Differentiate between 10-based Locator and Other Locators

texts:
In ESRI documentation, it states that the IAdvancedIntersectionGeocoding interface is only available for use in older locators (9.3.1 and below). I have verified this statement. One major problem is that 10-based locators still implements IAdvancedIntersectionGeocoding even though you cannot use it.

I want to be able to differentiate between locator versions to either be able to implement IAdvancedIntersectionGeocoding or not. Documentation states that ESRIFDOAddressLocator represents a pre-10 locator, while ESRIGen2AddressLocator represents a 10 locator.

I have been comparing the locator returned from the currently active locator to these two objects to determine versions, but they both return true (even if I select the 9.3.1 locators).

Is there any way to differentiate between locator versions? If not, is there any other interface or object I missed to be able to return standardized intersections?

Diagram:

EDIT:

Checking against ESRIGen2AddressLocator and ESRIFDOAddressEditor seems to be a problem only after Service Packs are installed. I have not drilled down to which one, but if ArcMap does not have a service pack, these checks work fine. I have logged a bug with ESRI to see if we can get this resolved.


texts: How to specify custom CSS to an OpenLayers Zoom Control in Angular 6?

texts:
I am attempting to modify the Zoom buttons on OpenLayers 5.3 in an Angular 6 component. And I'm using this OL3 question as a starting point: Creating new buttons and removing default ones - Openlayers-3.

However, my custom CSS class for my component is being ignored (or overwritten with the original ol.css for the zoom buttons).

this.map = new OlMap({
  target: 'map',
  layers: [
    new OlTileLayer({
      source: new OlOSM()
    })
  ],
  view: this.view,
  controls: [
    new OlZoom({ className: 'custom-zoom' })
  ]

And I've placed .custom-zoom and .custom-zoom button in my scss for my component, but they do not render for the Zoom buttons.

Any ideas on an Angular 6 solution to custom CSS classes for the Zoom buttons?


texts:
I just started in the R field. I installed R project and then R studio.

I want to execute a process with library (water) and I need before calling RGDAL and I get the following error:

library("rgdal", lib.loc="~/R/win-library/3.4")
rgdal: version: 1.2-8, (SVN revision 663)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
 Path to GDAL shared files: C:/Users/MaGregor/Documents/R/win-library/3.4/rgdal/gdal
 Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
 Path to PROJ.4 shared files: C:/Users/MaGregor/Documents/R/win-library/3.4/rgdal/proj
 Linking to sp version: 1.2-5


texts: Take pixel locations from rasterio raster and write out shapefile

texts:
I'm pulling a satellite image into python with rasterio, running some analysis to find a specific bounding box within that image, and then want to write those points out as a shapefile. I'm struggling with the Fiona documentation, much online seems out of date and I can't convert the official docs well to my example.

So basically I've got

import rasterio

dataset = rasterio.open(image_filepath)
print(dataset.crs)

# this outputs '+init=epsg:4326'

"""
run analysis to find points on raster
"""

# example of a single rectangle that we get from that analysis
geolocated_xs = (-113.21645410141903, -113.21645410141903, -113.2167506315605, -113.2167506315605)
geolocated_ys = (26.827766731681987, 26.828098563030775, 26.827766731681987, 26.828098563030775)

"""
write out to shapefile
I'm not totally sure what this should look like...
"""

# I assume I need shapely somewhere here?
from shapely.geometry import mapping, Polygon

import fiona
from fiona.crs import from_epsg
crs = from_epsg(4326)

schema = ?

# so from the docs

with fiona.open(
     'foo.shp',
     'w',
     driver='ESRI Shapefile,
     crs=crs,
     schema=) as c:
        c.write( ??? )

But I'm not sure what to write out and how I need to convert those tuples to have them go into this Shapefile.

I seek suggestions for pointers to good examples that I can use to learn.


texts: Adding map current location(Latitude and Longitude) in ViewParams of Geoserver WMS using Leaflet

texts:
I am drawing Geoserver layers on google map using Leaflet. Base database for Geoserver is PostgreSQL.In PostgreSQL there is a table that stores Millions of polygon.These polygons are in the form of square of 300 sq. feet each. I want to search polygons within 50 Miles in circular radius from map's center location. For this I need center location of map whenever user zoom in/out or drag the map.So that I can do a query in PostgreSQL on the basis of location. Here is the code.

   //layersearcharea = user will enter this value.
  //layersearcharea is area that will be used in circular radius search
  //algorithm or say KNN algorithm

    var viewparams = ['searcharea:' + layersearcharea,'centerlat:'+map.getCenter().lat,'centerlong:' + map.getCenter().lng];

                    var baseURL="http://localhost:8080/geoserver/wms?";
                     var baseURL=encodeURI(baseURL+"viewparams="+viewparams.join(';'));

                var wmsLayer = L.tileLayer.wms(baseURL, {
                    layers: 'MyWorkSpace:GetCellLayer',
                    transparent: true,
                    format:'image/png',
                    tiled:true
                }).addTo(map);

Question is How and where can I update centerlat and centerlong in viewparams whenever user change zoom level of map ?


texts: Conversion problem with 'Feature to Raster': resulting Raster not visible in ArcMap

texts:
I want to convert a Shapefile consisting of numerous point-features that each represent a specific height in my DEM to GRID-format so that i can then analyze it with Spatial Analyst.

Mysteriously, the conversion with "Feature to Raster" only takes 3 minutes (even tough the merging of the shape took 50) and the resulting GRID-layer is not visible in ArcMap or ArcCatalog.

There seen to be only two points in the resulting GRID, even though I picked height as value to assign information to the output raster.

What am I doing wrong here?

(When viewed completely, my Shape appears as a black area, but I assume that's because of the high density of point features. If I zoom in they are visible as separate points.

Shapefile-PRJ:

  PROJCS["DHDN_3_Degree_Gauss_Zone_3",GEOGCS["GCS_Deutsches_Hauptdreiecksnetz",DATUM["D_Deutsches_Hauptdreiecksnetz",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Gauss_Kruger"],PARAMETER["False_Easting",3500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",9.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

The desired resolution of my output-raster is 1x1 m.


texts: Progress bar does not show progress during processing

texts:
I am trying to show the progress in the progress bar, but it is not running, I have seen several solutions but none has solved me.

This is part of the code that I am using :

def progress_changed(self, progress):
    self.progressalg.setValue(progress)

def AdicionarDatosparcelas(self):
    self.usuario_origen = self.cfgdlg.v_usuario_o.text()
    self.contrasena_origen = self.cfgdlg.v_contrasena_o.text()
    self.usuario_destino = self.cfgdlg.v_usuario_d.text()
    self.contrasena_destino = self.cfgdlg.v_contrasena_d.text()

    progressMessageBar = self.iface.messageBar().createMessage("Insertando parcelas en la base de datos...")
    self.progressalg = QProgressBar()
    self.progressalg.setAlignment(Qt.AlignLeft | Qt.AlignVCenter)
    progressMessageBar.layout().addWidget(self.progressalg)
    self.iface.messageBar().pushWidget(progressMessageBar, Qgis.Info)

    f = QgsProcessingFeedback()
    f.progressChanged.connect(self.progress_changed)

    geom_column = None
    uri_origen = QgsDataSourceUri()
    uri_origen.setConnection(self.servidor_origen, str(self.puerto_origen), self.basedato_origen, self.usuario_origen,
                             self.contrasena_origen,
                             QgsDataSourceUri.SslDisable)
    uri_origen.setParam('key', 'mi_prinx')
    uri_origen.setParam('checkPrimaryKeyUnicity','0')

    uri_origen.setDataSource("catastro", "parcelas", geom_column)
    parcela_origen = QgsVectorLayer(uri_origen.uri(), 'parcela_origen', "postgres")
    QgsProject.instance().addMapLayer(parcela_origen, True)

    uri_destino = QgsDataSourceUri()
    uri_destino.setConnection(self.servidor_destino, str(self.puerto_destino), self.basedato_destino,
                              self.usuario_destino, self.contrasena_destino,
                              QgsDataSourceUri.SslDisable)

    uri_destino.setParam('key', 'mi_prinx')
    uri_origen.setParam('checkPrimaryKeyUnicity', '0')
    uri_destino.setDataSource("catastro", "parcelas", geom_column)
    parcela_destino = QgsVectorLayer(uri_destino.uri(), 'parcela_destino', "postgres")
    QgsProject.instance().addMapLayer(parcela_destino, True)

    processing.run("etl_load:appendfeaturestolayer", {
        'SOURCE_LAYER': parcela_origen,
        'SOURCE_FIELD': 'mi_prinx',
        'TARGET_LAYER': parcela_destino,
        'TARGET_FIELD': 'mi_prinx',
        'ACTION_ON_DUPLICATE': 1}, feedback=f)
    self.iface.messageBar().pushSuccess('parcelas exportadas', '')


texts: How to center a user-defined CRS on an area of interest?

texts:
My study extent in central Iran straddles two UTM zones (39 and 40):

I would like to create a user-defined CRS in Transverse Mercator projection with central meridian inside this area. Through Custom CRS>> Add new CRS, I try to add the following proj4 string:

+proj=tmerc +lat_0=Latitude of natural origin
          +lon_0=Longitude of natural origin
          +k=Scale factor at natural origin
          +x_0=False Easting
          +y_0=False Northing

How can I calculate these parameters from my study area? I use QGIS 2.6.1.


texts: Create history map with inkscape with each country as a distinct object

texts:
I want to create history maps (maps showing historical borders) with inkscape and to have an object for each country.

I want to draw borders between countries only once and have this part of a path duplicated and assigned to the bordering countries.

E.g. I want to draw the border between France and Spain only once and have it assigned to both countries.

I know that I can simply draw each border as a separate path, cut out a part or use the whole path, copy it and merge each copy with the rest of the respective country's outline.

Is there a better approach to the problem?

Is there a tool, e.g. an Inkscape addon, that can identify all path part combinations that surround and therefore define areas?


texts: How to install QGIS 2.99.x master on Raspbian GNU/Linux 8 (jessie)

texts:
I want to install QGIS 2.99.x master on a Raspberry Pi 3, but there is something going wrong. First I've added the following lines to /etc/apt/sources.list:

deb     http://qgis.org/debian-nightly jessie main
deb-src http://qgis.org/debian-nightly jessie main

Unfortunately, this installs QGIS 2.4 Chugiak after running sudo apt-get update and sudo apt-get install qgis python-qgis qgis-plugin-grass.

I noticed the following error in the terminal when running sudo apt-get update:

W: GPG error: http://qgis.org jessie InRelease: The following signatures couldn't be verified because the public key is not available: NO_PUBKEY 073D307A618E5811
W: Failed to fetch http://qgis.org/debian-nightly/dists/jessie/main/binary-armhf/Packages  404  Not Found [IP: 104.27.139.220 80]

E: Some index files failed to download. They have been ignored, or old ones used instead.

I tried to fix this by adding the key manually (see http://hub.qgis.org/issues/8404):

gpg --keyserver keyserver.ubuntu.com --recv 073D307A618E5811
gpg --export --armor 073D307A618E5811 | sudo apt-key add -

Now sudo apt-get update throws a new error:

W: Failed to fetch http://qgis.org/debian-nightly/dists/jessie/InRelease  Unable to find expected entry 'main/binary-armhf/Packages' in Release file (Wrong sources.list entry or malformed file)

E: Some index files failed to download. They have been ignored, or old ones used instead.


texts:
New to programming here. My legend (bottom left corner) keeps showing up blank even though there is technically a legend element there. Am I missing something in my code? I've included the relevant code and an image of my map.

lyrFile = arcpy.mp.LayerFile(meaninc)
aprx = arcpy.mp.ArcGISProject("E:/2016Fall/ProgrammingFundamentals/Final/Maps/Maps.aprx")
# Ask user to input their desired file name and location
location = input("Input file path:")
map = aprx.listMaps("Map")[0]
# Remove previous layer files
for l in map.listLayers():
  if l.name in names:
    map.removeLayer(l)
# Add desired layer file
map.addLayer(lyrFile)
# Turn on the legend???
lyt = aprx.listLayouts("Layout")[0]
legend = lyt.listElements("LEGEND_ELEMENT", "Legend")[0]
legend.visible = True
# Export PDF layout
lyt.exportToPDF(location)

aprx.save()
del aprx


texts: I have two large polygon data-sets, one highlights the building footprints and one indicating parcel outlines. I want to find all parcels without a building on them. Essentially all parcel A's which do not contain a parcel B. Is the Spatial query, select by location the best way to do this? The files are very large and often crash.


texts: Dissolve a table with no geometry QGIS return Null if any are Null

texts:
I have a table without geometry that contains three columns: Item, Sub-item, and amount. I would like to create a new table with the sum of the amount of the sub-items that have the same Item attribute. (These categories are invented)

Normally I would use the GDAL Dissolve with the Item as the dissolve field. But in QGIS it doesn't accept layers without geometry as an input.

I have created dummy point geometry to get around this but feel like I must be missing something. Could I just be using an sql query?
Also I would like it if any of the sub-items are null it would give me a null sum. Currently it only gives a null if all of the items are null.


texts: Sentinel 2A Classification in Google Earth Engine

texts:
I tried to get classify image for my region of interest using Sentinel 2 data. I got the result, but I am getting an error.

var roi = ee.Geometry.Rectangle([82.5,18,83,18.5]);

//this comment for filtering the sentinel2 images in order to cloud cover , date,roi
var data = sentinel2.filterBounds(roi)
      .filterDate ('2020-05-01','2020-06-01')
      .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',10))
      .select('B2','B3','B4','B8')
      .median();

print(data);

Map.addLayer(roi);
Map.centerObject(roi);
// Map.addLayer(data);

var cliproi = data.clip(roi);

Map.addLayer (cliproi);

var training = waterbody.merge(agriculture).merge(forest).merge(buildup);

// taking training samples
var training = cliproi.sampleRegions({
  collection: training,
  properties: ['landcover'],
  geometries:true,
  scale: 10
  });

print (training);

// Extract training data from select bands of the image, print to console
var bands = ['B2', 'B3', 'B4', 'B8'];
var training = cliproi.select(bands).sampleRegions({
  collection: training,
  properties: ['landcover'],
  scale: 10
});

// Train classifier
var classifier = ee.Classifier.randomForest().train({
  features: training,
  classProperty: 'landcover',
  inputProperties: bands
});

//Run the classification
var classified = cliproi.select(bands).classify(classifier);

//Add the classification to the map view, specify colours for classes
Map.addLayer(classified,{min: 0, max: 3, palette: ['#0f1ad6', '#98ff00', '#037b3b','#ff4d4d','#ffb570']},'Classified');


texts: Performing logarithmic data transformation for Sentinel-1 using Google Earth Engine?

texts:
I am trying to calculate the logarithmic ratio of VV and VH bands of Sentinel-1 images according following equation: log(Ratio) = log(VV) - log(VH). The code is running fine (no error messages), but I don't get any values for the ratio. Unfortunately, I can not find the error. I exported the test image 'LogRatio' and checked in GIS software, also as .csv in Excel, and there are no values at all.

// Filter to get S1-images with VV and VH dual polarization.
var S1_asc = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filterDate(start, finish)
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'))
    .filter(ee.Filter.eq('resolution', 'H'))
    .filter(ee.Filter.eq('resolution_meters', 10))
    .filterBounds(ROI)
    .select('VV','VH');
print('Collection: S1', S1_asc);

// log VV/VH ratio

// select one image from time series to create function.
var imageS1 = ee.Image(
  S1_asc.filterDate('2017-05-06', '2017-05-07')
    .first()
);

// create function on band math and expression()
var addLogRatio = function(imageS1) {
  var logVV = imageS1.select('VV').log();
  var logVH = imageS1.select('VH').log();
  var logratio = logVV.subtract(logVH).rename('logVV/VH');
  var logratio_float = logratio.toFloat();
  return imageS1.addBands(logratio_float);
};

// Test the addLogRatio function on a single image.
var logratio_float = addLogRatio(imageS1);
print('Test Image: LogRatio', logratio_float);

// Apply addLogRatio function to all image of image.collection.
var S1_asc_logRATIO = S1_asc.map(addLogRatio);
print('Collection: LogRatio', S1_asc_logRATIO);


texts: Creating a selection that selects points that are within a certain distance of other points in QGIS

texts:
So I am in a beginner Geospatial Analysis class which is focused on learning Arc and ESRI-owned GIS software but I wanted to be able to apply what I am learning in something that would be more applicable to me. The instructions are meant primarily for ArcMap but I am trying to follow them in QGIS.

I need to be able to select and count how many points are within one mile distance of a point in another set of points. The instructions for ArcMap say to use the select by location tool, but the one in QGIS doesn't seem to have the same functionality as in ArcMap. I have seen some talk about the distance matrix but I don't know anything about that.


texts: Getting strange exception since switching to GeoServer 2.16 from 2.8 with WFS queries in OpenLayers 2

texts:
Here is the error in the GeoServer log (as it is returned to the caller):

 2019-10-28 17:03:59,021 INFO [geoserver.wfs] -
 Request: getServiceInfo
 2019-10-28 17:03:59,053 ERROR [geoserver.ows] -
 org.geoserver.wfs.WFSException: Could not locate {http://my.namespaceurl.com}MyLayerName in catalog.

Here is the code I used to register a new OpenLayers.Protocol.WFS control:

var theProtocol = new OpenLayers.Protocol.WFS({
    version: "1.1.0",
    url: "http://my.geoserver.sample:8080/geoserver/wfs",
    featurePrefix: "MyWorkspaceName",
    featureType: "MyLayerName",
    srsName: "EPSG:900913",
    featureNS: "http://my.namespaceurl.com",
    geometryName: "the_geom"
});
var theFeatureControl = new OpenLayers.Control.GetFeature({
    protocol: theProtocol,
    maxFeatures: 2,
    clickTolerance: 10,
    box: false
});
theFeatureControl.events.register("featureselected", this, function (e) {
    ...
});
theFeatureControl.events.register("featureunselected", this, function (e) {
...
});
theMap.addControl(theFeatureControl);
theFeatureControl.activate();

And in GeoServer, I have a workspace named "MyWorkspaceName" with a namespace of "http://my.namespaceurl.com" and a layer named "MyLayerName" ...

This all worked with a GeoServer 2.8, now I've installed a new server with the lastest version (2.16 at this time) and created the same workspaces, stores, and layers.  But this one function seems to be getting tripped up somehow, instead of finding "MyWorkspaceName:MyLayerName" it seems to be using the namespace name??

Anybody ever had a similar thing happen?  Any hints/ideas?


texts: How can I get an average value for area using 30 different vector shapefiles as input?

texts: I have 30 different vector boundary polygon shapefiles for a study area that don't exactly match up.  Each polygon has a risk score.  I want to overlay all 30 maps into a single map that represents the average risk score of the areas.  I think I need to use vector to raster in ArcMap on the 30 shapefiles, but then I'm not sure how to overlay them.  I've used zonal statistics before, but I'm not sure if that is relevant to this kind of project. I think the final output will have to be a raster since the boundaries of the 30 input files do not match up exactly.


texts: Getting Canada shapefile to overlay model on QGIS?

texts:
How can I get Canada to overlap with the rest of my map?

I think this is a CRS problem but I really do not know how to solve it.


texts: Quickmapservices plugin not showing Google or Bing layers

texts:
Ive installed QGIS 3.10.14 on a new windows machine and then the Quickmapservices plugin. However, when I go to the menu, it doesn't show google or bing options. Were they discontinued along the way? Is this a known issue? Is there a solution?
TYVM!


texts: Replacing urls in QGIS attribute table with images using python

texts:
With the help of python, I managed to build a small script that replaces all URLs in a text column with images using a delegate. The images where loaded from Microsoft SharePoint using an adapted version of request_ntlm.

import PyQt4
from PyQt4 import QtGui,QtCore
import requests
from requests_ntlm import HttpNtlmAuth, HttpNtlmSspiAuth
from PyQt4.QtGui import QApplication, QTableView

class ImageDelegate(QtGui.QStyledItemDelegate):
    def __init__(self, parent):
        QtGui.QStyledItemDelegate.__init__(self, parent)

    def paint(self, painter, option, index):
        # Get Item Data
        url = index.model().data(index, QtCore.Qt.EditRole)

        if url != '':
            style = QtGui.QApplication.style()
            opt = QtGui.QStyleOptionGraphicsItem()
            opt.rect = option.rect
            r = requests.get(url, auth=HttpNtlmSspiAuth())
            data = r.content

            pixmap = QtGui.QPixmap()
            pixmap.loadFromData(data)
            if option.rect.x() > 250:
                pixmap=pixmap.scaledToWidth(250)

            painter.drawPixmap(option.rect.x(), option.rect.y(),  pixmap)

dialog = iface.showAttributeTable(iface.activeLayer())
targetColumnIndex = 8 # how to find Index of column?
delegateImage = ImageDelegate(iface)

dialog.findChildren(QTableView)[0].setItemDelegateForColumn(targetColumnIndex, delegateImage)
dialog.findChildren(QTableView)[0].setColumnWidth(targetColumnIndex, 250 )
dialog.findChildren(QTableView)[0].verticalHeader().setDefaultSectionSize(200)

What I couldn't manage by now: how to retrieve the index of a given column (see hardcoded targetColumnIndex in the code above) in the attribute table (QTableView) even when "colxy" is moved or some columns are hidden?


texts: What web server is embedded in ArcGIS 10.1 for Server?

texts:
Esri seems to have a new embedded webserver in ArcGIS 10.1 for Server.

Does anyone know what technology they have used to build it?


texts:
Assume we have a well populated PostGIS database with lots of metadata about raster and feature based data.

What open source technologies would one use to create a end-user friendly catalog search and download system (web based).

Something similar to what is located at: https://catalog.data.gov/dataset

(I realize they just provide pointers to the original source, we would be caching the rasters and the feature in a datalake, and using PostGIS for the catalog/index)
yea, the data.gov uses ESRI, but let's take that off the table.

If I had to do it by hand, I would use Leaflet or OpenLayers, but since this is such a common need, I wonder if someone has not already put that together).


texts: Connecting to lines, which are not spatially connected to each other, using FME?

texts:
I have a shapefile, which contains lines. These lines are not spatially connected to each other.

More specifically the end of one linestring has a slightly different coordinate from the start of the next linestring (5th or 6th decimal digit).

Is there a way to connect these lines together using FME 2014?


texts: How to compute orthodromic distance in QGIS?

texts:
I have opened a geotiff file with QGIS. The CRS of this file is EPSG:3857 (web mercator).
When I use the distance tool in the menu, the computed distance is towards web mercator coordinates and not the "real" distance on the ground.
I have found the GeoServer extension to compute the orthodromic distance. It works but it's not really easy to use.
Is there a way to change QGIS parameters or anything else that allows to get orthodromic distance with the native tool?


texts: Making QGIS attribute table field automatic?

texts:
I'm working in a way to make hydrologic projects using QGIS and an Excel worksheet that I have. To do this, I want to extract some information of the lines, included in a vector layer, which represents a pipe section.

The information that I need to extract are:

ID Number
Length
X, Y start and end coordinates

I've found a way to capture this field using the "$length" and another algorithm for the X and Y coordinates, but for that I need to open the Attributes table, put the expressions in each attribute column and click to update fields.

Is there a way to when I draw a line these fields are filled in automatically? That is, I draw / edit a line (start editing or end node ) and when I open the attribute table the length fields and coordinates X , Y are filled and updated.


texts:
When using the heatmap plugin, I am getting some insanely large numbers for the cell size field when setting the radius to 5000.  I know that rows and columns are related to cell size x and cell size y, but I don't know how exactly.

Anyway, when I click OK and run the plugin, I get "nan" as a layer on the raster.  Check the pictures below to see what I mean.


texts: Sharing classes and methods code among processing scripts

texts:
I am writing a set of processing scripts for seismic data acquisition (2D and 3D) planning and geometry processing.   They are very crude scripts (I'm not a developer) that performs the tasks I meant to.   At first they ware developed in the processing environment of QGis 2.x, and I'm porting now to QGis 3.X (working with 3.10 now).

There are a number of functions/classes/methods that are used by some of them.  I would like to be able to put them inside a directory and just import were needed.  Is there any variable in the processing environment pointing to the scripts directory?
e.g.: ~/.local/share/QGIS/QGIS3/.../scripts

My point is that if I could put those utility code in a directory in the scripts folder, it would be possible to insert that directory in PYTHONPATH to have it available during the script run.

Let me know if there is any other approach that can make it possible.  The goal is to have a dynamic way to change it according to the platform or personal configurations.


texts: Reproject raster catalog - ArcGIS 9.3/FME 2012

texts:
I have a couple of raster catalogs in different file geodatabases and I need to re-project them to a new coordinate system. Now there is no option to project a whole raster catalog, nor can I use the 'Project Raster' tool, as it won't accept a raster catalog as file type (for selecting single rasters within the catalog).

There are two options that I thought about so far:

use tool 'Raster Catalog to Raster Dataset', reproject the raster dataset and bring it back into a catalog. Now the problem is on the one hand that this process takes ages (incl. building of pyramides + calculating statistics) and I loose the single items within the raster catalog (after loading the dataset back into a raster catalog), as it just mosaics the datasets together.
use Python and loop through raster datasets in raster catalog and reproject them individually. I have not investigated this version, but I doubt whether that will work. In case it would, that could be a solution. Anyone did something like that before?

My question therefore is, whether there is a proper/efficient way of reprojecting a raster catalog or the raster datasets within a raster catalog respectively from one coordinate system to another? I basically want to end up with raster catalogs that look exactly as the originals (incl. file names of raster datasets within them), but they do need to have a different projection.

Any tips appreciated as my problem is quite urgent!

UPDATE:

I have access to FME 2012 now, never used it before though. Does anyone have an idea whether this exercise could be accomplished using FME? (it looks like FME can only access SDE raster catalogs, so that may not be the full solution either..)


texts: Assigning spatial info for NetCDF files using ArcGIS for Desktop?

texts:
I have quite a few NetCDF files that arrive from a air quality model without any spatial reference.  However, I do have the spatial reference information that is required.

Within Arc I use "make a netcdf raster layer", it has "undefined" spatial info in properties.
Within the undefined netcdf properties I can view several of the chemical variables like NH3, C02, etc in the drop down menu.

Defining the projection does not change the spatial info and returns the same "undefined"

I have been able to assign a spatial reference by converting that layer to ASCII and then creating a prj and manually adding the spatial info needed to the files.  This creates a chemical layer that is spatially correct...
However this process does not keep the several chemicals in the netcdf drop down mentioned earlier but only one chemical...and this is a slow process.
Is there an easier way to go about this instead of  displaying in Netcdf in arc as raster, then converting to ASCII and then manually adding in spatial info?  And still be able to keep all the variables in the drop down menu in the Netcdf tab?


texts: How to export point geomtery geojson from QGIS?

texts: I want export geojson (from shapefile), for leaflet project. But QGIS 2.6 exported geojson is multipoint geometry(leaflet doesn't read that), so its way to export point geometry geojson from QGIS?