texts: Attribute table shows histogram/pixel count for a thematic raster image, but on enquiring, every pixel's value is zero

texts:
I have a thematic raster image containing 18 land use classes (like built-up, water bodies, plantations, forest, etc.). When I display this image in a GIS application (Arcmap/Erdas IMAGINE), the map is displayed with all the classes in the attribute table. The attribute table has the histogram/count column for each class and the area of each class that I have computed. But when I enquire a pixel, its value is zero. When I try to run an algorithm on this image, the output is a blank image (all zeros), because the pixel values in my original image are zeros.

The original values were present to begin with, because I had worked with the image in the past. Only now, the image with the colors I had assigned to different classes is visible and the attribute table has the original data, but the pixel values are all zeros.

What could be the problem? And is there a way to restore this raster image?


texts: Loading directory of JSON files with PyQGIS

texts:
I'm trying to load a directory of JSON files into QGIS 3.x

They all have the same prefix of 'land'. My code to target all the files starts with:

# 1. import glob library
import glob

# 2. Targets files with prefix
layers = glob.glob('V:/GIS - Files/3. Data/Agriculture/land*.json')

Each json file contains a line and polygon sub-set of data. I'm trying to load both geometry types separately. Is there a way of doing this?

Geometry types for each JSON file:


texts: Interpolate line and profile graph error, not sure how to fix

texts:
I am using the "Interpolate Line" tool in 3D Analyst to create a profile using "Profile Graph". However, I get two error windows:

 1.

  Unable to create temporary shapefile

 2.

  To create a profile graph, the selected graphic or feature has to be a
  3D line. To create a 3D line, use the Interpolate Line tool or Line of
  Sight tool.

My DEM has been selected.

Workflow:


texts: Getting junction id for roads in network dataset using ArcPy?

texts:
I have a network dataset with roads and junctions.

I am wondering what class I can use in ArcPy to extract FromJunctionId and the ToJunctionId for each road segment in the network dataset.


texts:
I'm using gdal_grid to interpolate IDW mode with points in GPKG format but I don't know how to write the sequence. I'm using the following code:

gdal_grid -zfield "z" -a invdist:power=2.0:smoothin=1.0 -txe 465000 750000  -tye 4852900 4623700  -outsize 6000 5000 -of GTiff -ot Int32 -l "deliverynote_specie_name_E" C:\coberturas\temp\deliverynote_specie_name_E.nitens.gpkg "C:\coberturas\temp\Enitens.tif"

But gdal gives me the following error:

  ERROR 1: bad application_id on
  'C:\coberturas\temp\deliverynote_specie_name_E.nitens.gpkg' Unable to
  open input datasource
  "C:\coberturas\temp\deliverynote_specie_name_E.nitens.gpkg". bad
  application_id on
  'C:\coberturas\temp\deliverynote_specie_name_E.nitens.gpkg'

I think the problem is how to write the name of point layer of my GPKG.

on the other hand, how could I define the spatial resolution of 100 meters?


texts:
I ST_Dwithin the function returns TRUE if the geometry A is completely inside B. I happen to have many polygons (representing regions) and need to know which streets are within each of these polygons. But it has some streets that hath part in a region and part in another.
Which function is most appropriate to select postgis geometries that have any part inside another?


texts: r.viewshed - Can't find with Grass toolbox

texts:
I'm trying to use r.viewshed and have got Grass installed but can't find it anywhere!  under toolbox, I click 'Grass commands (152 geoalgorithms) and have a choice of imagery, raster, vector, etc, etc. under raster I have a long list of tools, but no r.viewshed... anyone know where I'm going wrong?

David.


texts: Seeking method for automatic aggregation of polys with sparse data?

texts:
I am part of a project on presenting health related data on maps. In order to maximize the analytical usefulness of the data, we wish to divide them into as small areas as possible. However, to protect individuals we also have to make sure that each area contains enough cases to secure anonymity. We are therefore looking for methods to automatically aggregate an area with too sparse data to a neighboring area – or perhaps a bit more advanced – to a non-neighboring  area in the map, which according to sociodemographic parameters resembles the area in question.

I seek any input on methods for automatically:

Aggregate an area with too few cases to a neighboring area – for example the neighboring area with the lowest number of cases.
Aggregate an area A to the area B – neighboring or non-neighboring – which on a number of parameters has the closest resemblance to area A.

The work will be done in QGIS and/or ArcGIS Desktop, so any readily implemented methods in any of these systems are welcome, but general theoretical descriptions of statistical methods – probably most relevant for the second case – are also of interest.


texts: Looking for a functionality in QGIS to find all neigbours in a polygon layer

texts:
I am looking for a functionality in QGIS to find all neigbours for each polygon in a polygon layer. The result should be a table with 2 ID-Columns: one for the polygon and the second for the neighbours.

I tried "spatial query" in QGIS with the layer as source and its copy as reference, but to no avail.

In Mapbasic of Mapinfo it would be a simple expression:

select t1.id1, t2.id2 from t1, t2 where t1.obj intersects t2.obj

(where t1 and t2 are two copies of the polygon layer).


texts:
I have a table which contains an field zoneShape of the type Polygon.

When inserting some geometry in this field, in phpmyadmin it says: affected 1 row. The row is added, but the geometry field says it's size is 0B. Which can't be accurate, cause it has multiple points.

This is the Insertion code:

SET @g = 'POLYGON(50.866753 5.686455, 50.859819 5.708942, 50.851475 5.722675, 50.841611 5.720615, 50.834023 5.708427, 50.840744 5.689373, 50.858735 5.673923, 50.866753 5.686455)';
INSERT INTO zone SET zoneShape = PolygonFromText(@g)

As you can see; the first and last point are the same, so the polygon is closed.

Can anyone help me with adding this polygon to the database?
Thanx


texts: Removing file from os results in WinError 32

texts:
Using PyQGIS for a standalone plugin, I try to remove a file I created before as an output file for QGIS processing:

    layers = [layer for layer in QgsProject.instance().mapLayers().values()]
    inputIndex = self.dlg.layerlist.currentIndex()
    input = layers[inputIndex]

    output = self.dlg.directory.text()

    # run algorithms
    params_snap = {
        'input' : input,
        'type' : 1,
        'tool' : 1,
        'threshold' : 0.000095,
        'output': output,
        'error' : 'TEMPORARY_OUTPUT',
        'GRASS_OUTPUT_TYPE_PARAMETER': 2
    }
    processing.run('grass7:v.clean', params_snap)
    context = QgsProcessingContext()
    input = QgsProcessingUtils.mapLayerFromString(output, context, True)
    QgsProject.instance().removeMapLayer(output)
    if os.path.exists(output):
        os.remove(output)

But I get a [WinError 32] The process cannot access the file because it is being used by another process. As you see in the code I tried adding QgsProject.instance().removeMapLayer(output) to remove the file from QGIS first, but it doesn't work. Any advice?

What I try to do by deleting the output file from the algorithm is to run a second (and third...) algorithm after the first one, using the same file as the output.


texts: Drawing polygon and summarizing values in map on Portal for ArcGIS?

texts:
I have a map application created in Portal for ArcGIS. I have a points layer in the map. Each point has a numeric value. I want to be able to draw a polygon on the map, containing some of the points. It is fine if it is just a markup. I want to sum the numeric values from the points in the polygon and show the result in the app.

Are there any widgets or other ready made functionality I can make use of?


texts:
I have to reproject a shapefile's coordinate system. To do that I am using Arcpy and I want to get an output to use it as input in other tool (in Model Builder).

Now I get an output, but I cannot use it as input in model builder's tool. I am not programing well when I process the output shapefile.

The code is:

def getParameterInfo(self):
        """Define parameter definitions"""

 # First parameter

    param0 = arcpy.Parameter(
        displayName="Shape entrada a reproyectar",
        name="in_features",
        datatype="DEShapefile",
            parameterType="Required",
        direction="Input")

# Second parameter

    param1 = arcpy.Parameter(
        displayName="Imagen satelite de la cual se utilizara el CRS",
            name="raster",
        datatype="GPRasterLayer",
            parameterType="Required",
        direction="Input")

# Third parameter
     param2 = arcpy.Parameter(
        displayName="Shape de salida",
        name="out_features",
        datatype="DEShapefile",
            parameterType="Derived",
        direction="Output")

    parameters = [param0, param1, param2]
    return parameters`

def execute(self, parameters, messages):
        """The source code of the tool."""

shape= parameters[0].valueAsText
raster= parameters[1].valueAsText
output=parameters[2].valueAsText
#-----------------------------

desc = arcpy.Describe(shape)

out_CRS = arcpy.Describe (raster).spatialReference

output=arcpy.Project_management(shape, desc.path + "\\" + desc.baseName + "_Rprj.shp", out_CRS)

return


texts: Getting max() from list of dates using ArcPy

texts:
Using Python 3.6.8

I had the below portion of my script working previously on some test data, but now that I've tried to spool it up to the production data, it's not working. This is a snippet of a much larger script, but I believe my error is somewhere in these lines of code below. The field that the SearchCursor is searching is a date type field and every record has a date in it (no nulls/none).

It is looking at a list of dates and should be assessing the 'max' (most recent) but on line 25 it's returning:

  TypeError: '>' not supported between instances of 'NoneType' and
  'datetime.datetime'

Here is my pared-down code:

import arcpy
import datetime
from datetime import time, timedelta

# Environment
arcpy.env.overwriteOutput = True
arcpy.env.workspace =r"C:\GIS\CAD_reporting\CAD_reporting.gdb"
arcpy.env.scratch = r"C:\GIS\CAD_reporting\CAD_reporting.gdb"

workspace = r"C:\GIS\CAD_reporting\CAD_reporting.gdb"
ptFC = r"C:\GIS\CAD_reporting\CAD_reporting.gdb\RMI_pts"
dateList=[]
nowDT = (datetime.datetime.today()).strftime("%m/%d/%Y %H:%M:%S")
lowDate = datetime.datetime.today()

# START
print ("Start: " + nowDT)

# Get date for the query start(low date)
dtCursor = arcpy.da.SearchCursor(ptFC, "ReportQueryDate")
for row in dtCursor:
    dateList.append(row)
del dtCursor

maxDateTuple = max(dateList) # TypeError: '>' not supported between instances of 'NoneType' and 'datetime.datetime'

lowDate = maxDateTuple[0] - timedelta(days=200) # subtracts 200 days
lowDateString = lowDate.strftime("%m/%d/%y %H:%M:%S")
print (lowDateString)

I can tell this error has to do with the data types, but I'm not sure how to handle my variables, to get them into the same/correct format for comparison.


texts: Coordinates of bounding box in Map composer

texts:
I would like to find a way to export the coordinates of the bounding box inside map composer. I know i can create labels with coordinates but I'm not able to define the lat/long of the corner of the map (fig 1)

I also followed this guide (https://gisunchained.wordpress.com/2014/09/09/coordenadas-dos-cantos-do-mapa-em-qgis-map-corner-coordinates-in-qgis/) but it gives the entire shapefile corners, aka the whole world in my case, instead of the area selected. I would like to have a similar output to this (fig 2) that i have in  tilemill.


texts:
I am writing a arc engine application that can take a shapefile, and merge a set of polygons together.

Here is the code that I am using to merge the polygons:

private static void mergeFeatures(IFeature merger, IFeature mergee)
{
    ITopologicalOperator union = (ITopologicalOperator)mergee.Shape;
    mergee.Shape = union.Union(merger.Shape);
    merger.Delete();
    mergee.Store();
}

Is this the best solution to merge polygons? In my application this part of my code takes 60% of the execution time.

Is it slow because I am using a shapefile? Does a file geodatabase help with this? Or is this even just the totally wrong approach for merging polygons.

Just to clarify, I am taking a large set of polygons and am trying to merge ones under a certain number of acreage. From these polygons that are small, I need to merge polygons with similar attributes (for example on a vegetation layer, ones that are similar vegetation.) So this can incur a large number of small polygons getting merged together, and eliminating some of the merges can be difficult. If there is an easier way to bulk merge, then that would be best, but thus far, I have not found one.


texts:
I am trying to add a legend to my leaflet map and no matter what I try I cannot get it to display properly.  I have imported a point level geojson into my map and I have styled the points in a similar way to how it's done within the Leaflet Choropleth example:

function getColor(d) {
        return d === 'Road Surface'  ? "#de2d26" :
               d === 'Signage'  ? "#377eb8" :
               d === 'Line Markings' ? "#4daf4a" :
               d === 'Roadside Hazards' ? "#984ea3" :
                            "#ff7f00";
    }

    function style(feature) {
        return {
            weight: 1.5,
            opacity: 1,
            fillOpacity: 1,
            radius: 6,
            fillColor: getColor(feature.properties.TypeOfIssue),
            color: "grey"

        };
    }

I can't add the code for the data because it contains sensitive information but I have added the geojson and have added some popups to it.

When it comes to the legend here is my code:

var legend = L.control({position: 'bottomleft'});
    legend.onAdd = function (map) {

    var div = L.DomUtil.create('div', 'info legend');
    labels = ['Categories'],
    categories = ['Road Surface','Signage','Line Markings','Roadside Hazards','Other'];

    for (var i = 0; i < categories.length; i++) {
            div.innerHTML +=
            labels.push(
                ' ' +
                (categories[i] ? categories[i] : '+'));
        }

        div.innerHTML = labels.join('
');
    return div;
};

legend.addTo(map);

This code brings up a box that contains a title and the categories that I have defined, but no colours, which is the thing that I can't work out.


texts: How to merge layers into a single QgsVectorLayer?

texts:
I have a set of different layers (and their respective shapefiles) with information of electric lines. All the layers have the same Fields in their Attribute Tables, but contain information about different features. I want to merge only the layers that I choose and save them into a new QgsVectorLayer (Or something similar).

To get all of the layers, I'm using this:

    lyrs = self.iface.legendInterface().layers()

And I want to get something like:

    new_layer=lyrs[1]+lyrs[3]
    new_layer.getFeatures() #Do something else from here, but is irrelevant to the question

(Obviously that last part won't work because "TypeError: unsupported operand type(s) for +: 'QgsVectorLayer' and 'QgsVectorLayer'", but that is the idea of what I need to get)

Some details:

1: I would prefer a method that does not create a new shapefile like "qgis:mergevectorlayers" (The method will be called multiple times, and it would create a lot of unnecesary shapefiles)


texts: How to save in a given XML schema in QGIS?

texts:
Related to this closed question: https://gis.stackexchange.com/q/184033/64092, I'd like to know if there is a way to save data in a particular GML/XML schema (say, any of these) using QGIS.

A step-by-step procedure would be best, as I'm not even sure I understand correctly whether those are "application schemas" (and whether it really matters).


texts:
I just started to play with leflet/geojson a little. But my coordinates are not rendered properly and I have no clue what is going on.

My coordinates are: 52.23943, 4.97599. They work correct with the setView function.

var map = L.map('leaflet_map').setView([52.23943, 4.97599], 15);

But using a geojasonFeature they are, hmmm, 'projected', somewhere east of Somalia.

var geojsonFeature = {
            "type": "Feature",
            "properties": {
            "name": "Coors Field",
            "amenity": "Baseball Stadium",
            "popupContent": "This is where the Rockies play!"
            },
            "geometry": {
                "type": "Point",
                "coordinates": [52.23943, 4.97599]
            }
        };
        var myLayer = L.geoJson().addTo(map);
        myLayer.addData(geojsonFeature).bindPopup("I am a gjson point.");

Anyone who knows what is happening here?

EDIT

Out of pure curiosity I changed the coordinates around [4.976143930893815,52.23925499011473] and the point appears at its correct location. A known bug!?


texts: (How could one use GIS software to determine...) What latitude has equal land areas north and south?

texts:
There's obviously more land in the northern hemisphere than the southern, but at some point, there must be a line of latitude where the area of land above and below is exactly the same.  Is it the same latitude as the centroid of Earth's landmass?

See: What is Centroid of all lands of Earth?


texts: Splitting a polyline at vertices and calculating start and end measurements for new features based on a known measurement of the original polyline

texts:
I have a route network which has a length-attribute defining length of each feature. This length-attribute comes from a data register outside of gis.

Now I need to split each original feature at vertices and ad a start and end measure for each new feature. For example:

A polyline feature is 1000 meters long as a true coordinate based length but in the data register it's length is 1100 meters. Now this polyline consists of 5 lines which are 200, 250, 150, 300, 200 meters long. Now I need to split this polyline at vertices (forming 5 new features) and ad start-end measurement for each based new feature on the length of 1100 meters.

What tool should I use in ArcMap 10.3.1?


texts: Select only the lines that have the a group letters

texts:
I have the following code:

from qgis.core import *

import processing
layer = QgsVectorLayer('C:\Test\Shapes\box.shp',
'box', 'ogr')
QgsMapLayerRegistry.instance().addMapLayer(layer)

query ='"Name" = "AAA"' and '"Func" != 1'

selection = layer.getFeatures(QgsFeatureRequest().setFilterExpression(query))
layer.setSelectedFeatures([k.id() for k in selection])

And Im trying to select only the lines that have the letters "AAA" in Name. The problem is that, this select everything that have that have similar code, like "AA" and "AAB, but i only want the "AAA" ones.

Can any of you help me?

I have tried different things but nothing work.


texts: Can not get Leaflet map to show via link in Rails 4

texts:
I'm stumped here, and wondering if anybody can help point me in the correct direction.  For that matter, I'm wondering if anybody has ever bumped into this?

When I enter my url in manually the map loads just fine(http://localhost:3000/page/map), but when I create a link in rails the map does load <%= link_to "Map", page_map_path %> or eve <%= link_to "Map", "http://localhost:3000/page/map" %> the page loads, but doesn't show the map until I manually do a refresh.

Here's a link to the repo on bitbucket.

From my testing this happens if I use a gem (either ver 1.0, or 0.7), or load leaflet directly.


texts: Using Intersect and dissolve to divide district into certain regions with data of interest without losing spatial data in ArcGIS Desktop?

texts:
I am using an intersect/dissolve feature workaround in ArcGIS Desktop 10.5.

I have a layer of districts of a city and another layer of parks in the whole city /layer of noise amount around big streets. I want to compute the coverage of that two features in comparison to the whole district area.

For doing so, I first intersect my layer of parks/noise-layers with the districts layer to get the parks/... in that certain district. After that, I dissolve the outcome for the districts to only have the districts as my main features, as these are the main interest.

My problem arises, when there are no parks or noise layers in a certain district because then the intersect logically cannot intersect anything, so no feature for a certain district is created. When I later want to show the quantities of the coverage to compare the values I have, that certain districts shows up as empty field.

For the comparison I join the dissolved parks/noise-layer on the first district layer over the single districts. As described above, the intersect and with that also the dissolved parks/noise-layer, has for certain districts no data. So in the end these data show up in the districts layer as null values, because they don't exist.

I do not have any idea how to preserve a certain district even if he has no parks/noise-layers inside. I should have value zero but I cant find any option for keeping this, and because that's the scope of the intersect function I would not make sense.

My Main Focus is every single district and if it in the end does not show up I lose information.

Are there other functions that match my need?

Or did I forget something in my workaround?

the image shows the districts where no noise-layer is intersecting with the district layer so in the end it does not have any data to be shown
the table shows the joined district layer with the coverage layer from the noise-layer


texts:
Has anyone added the Ocean Base map to QGIS?

I was able to add the OpenLayers plugin and get the "land" basemap (among others), but not the ocean.


texts: Reversing the order of features in a layer using OGR

texts:
I'm working on rasterizing some polygons using gdal_rasterize. This works great, but I want the rasterize command to place features with lower FIDs on top of those with higher FIDs. (Gdal_rasterize currently takes the latest polygon and burns that value onto the raster).

I can think of a way to do it by manually swapping the rings for the polygons, but is there a faster way reverse the order of features in a layer?

Hoping for a simple reverse sort function? I'm working with OGR/GDAL on Python 2.x


texts:
Is an SLD filter like this valid?

<ogc:Filter>
        <ogc:PropertyIsEqualTo>
        <ogc:PropertyName>aaaaaa</ogc:PropertyName>
        <ogc:PropertyName>bbbbb</ogc:PropertyName>
        </ogc:PropertyIsEqualTo>
</ogc:Filter>

or is a Literal always required?


texts:
I'm testing OpenLayers 6 beta.
ol.inherits looks missing in OpenLayers 6.

ol.inherits(ol.source.NowMap, ol.source.OSM);

This causes error

Uncaught TypeError: ol.inherits is not a function

What should I use instead?


texts:
I'm trying to parse some funky address data into separate fields either using ArcMap's Field Calculator or a Python Script. Currently, my address data including street number, street name, and unit/apt/suite number are all in one field. Some examples of the diversity in my data include:

1111 2ND ST
1753 21ST ST 2ND FLOOR
2125 ARIZONA AVE #303
2424 12TH ST SUITE 100

Thankfully, there are only a handful of suffixes which include:

AV | AVE | BL | CT | DR | LANE | PKWY | PL | RD | ST | TER | WAY

What python methods could I employ to search for these unique identifiers and then split whatever follows them (e.g., "2ND FLOOR", "#303", "SUITE 100")?


texts: How to compare two string field values within a given string value in ArcGIS10

texts:
I have a polygon shapefile which displays villages [NAMING]. On its attribute table I have, in one Field, their names and for each feature and I have its municipality too [POLYGON_NM].

Now what I want to do is to find if I have features repeated into the same municipality.

I wanted to work with ModelBuilder, and I've tried with some scripting code into Field Calculator, but it's a disaster more than success.

Here is a screenshot:


texts: Filling raster gaps in ArcMap based zone layer?

texts:
I have the following raster layer that has data gaps. I can use focal statistics to fill them:

Con(IsNull("raster"), FocalStatistics("raster", NbrRectangle(5,5, "CELL"), "MEAN"), "raster")

However, I want to use a zonal layer and limit the neighbors based on that layer so that the function only uses the mean of all the neighboring values that are in the same zone.


texts: How to group individual groups of pixels into a single class

texts: I have a raster image with certain pixel groups of equal area. Each group represents a single house. I want all the pixels near a group of houses to have the same pixel value, i.e., I want to name that area "residential area". How can I do this? I want something like this done: this image should change to:


texts: Defining the correct AEA proj4string for FEWSnet Rainfall Data (Southern Africa) in R

texts:
I have some general questions about how to go about solving this problem:

I am trying to incorporate FEWSnet's monthly Rainfall Estimate data from 2014-2015 for Southern Africa into an analysis where I am using a combination of survey data that I collected myself (projected in long/lat coords, WGS84 ellipsoid), the FEWSnet data, and some other rasters/SPDF's that already have projections defined (these are not my concern - I have been able to use projectRaster() and spTransform() to coerce these to new coordinate systems as needed).

The FEWSnet data are challenging because according to their documentation (http://earlywarning.usgs.gov/fews/product/120) they are published in an "Albers Equal Area Conic projection that uses the Clarke 1866 spheroid". But, when I download and read the files into R and call crs() on them, I get "CRS Arguments: NA" in return.

Still being a newbie to these things (and being largely self-taught) I'm confused about exactly what information to use to inform how to correctly set the projection of these data. I've tried looking for an EPSG code for AEA projections in Africa, as well as simply searching the internet for a solution, but am coming up empty. Here is the code that I used to import the data:

# Assume that data are downloaded from the above link, and stored in
# a directory on the desktop entitled "Map_Layers". Set this as wd.
# "Map_Layers", in my case, contains data from various sources including FEWS
# and within FEWS subdirectory are further directories for the raster files
# representing year, month, and dekadal periods of rainfall data.

DecRFE1 <- raster("./FEWS/2014_Dec_dd1_RFE/s14121rf.bil")
DecRFE1 # printed output below...
#    class       : RasterLayer
#    dimensions  : 530, 640, 339200  (nrow, ncol, ncell)
#    resolution  : 8000, 8000  (x, y)
#    extent      : -1456000, 3664000, -4224000, 16000  (xmin, xmax, ymin, ymax)
#    coord. ref. : NA
#    data source : /Users/HenryHWells92/Desktop/Map_Layers/FEWS/2014_Dec_dd1_RFE/s14121rf.bil
#    names       : s14121rf
#    values      : 0, 255  (min, max)

list.files("./Map_Layers/FEWS/2014_Dec_dd1_RFE")
# [1] "s14121rf.bil" "s14121rf.blw" "s14121rf.clr" "s14121rf.hdr"

So, consider this question a general plea for help on several things:

How to choose the correct CRS given the information in the files themselves and in the documentation on FEWSnet's website.
What are the uses of each of the four file types that the FEWSnet data downloads give you? (The .bil files are not in a format that I'm familiar with...)
If the best way forward is to set a custom projection, what information do I need to do so? Is there an EPSG crs I could take advantage of?


texts: So I have lines(road network) and buffers around them. I need to find "empty bufers" whose lines is missing (pic below). Lines have attributes of buffer width, but it is incorrect in some cases (does not match with actual buffer with). So creating/erasing with new polygons is not an option, since it will leave alot of trash polygons. Buffers is also dissolved in most cases.


texts: Adding symbol to multiple layers with auto refresh in QGIS 3

texts:
How to add a symbol to multiple layers at once and how to set up that the layers are refreshed once the symbol is changed?

I am aware of this thread Change style properties of several layers simultaneously?. Perhaps it is a duplicate, but the question was in 2012! Now it is 7 years later ... I thought that in the meantime the workflow may be changed, is not it?


texts:
On my desktop PC, I am using Ministry of Transportation Ontario road shields to label major highways on my maps.  I often need to work with a laptop when I'm visiting a client at their location, so I copied the appropriate style files to the same folder on my laptop.  When I try to apply the roads shield markers on the laptop, it gives me weird symbols and text instead of the appropriate shield.  I had no problem doing this on my desktop, so I must be missing something.

On the left (in the image below), I have an example of the road shield marker style that displays on maps on my PC.  On the right (column highlighted in yellow) are the marker symbol options I get in the style manager menu when I try to access these shields when using the laptop.  I assume this is a simple fix and I'm overlooking something.


texts: Create layer from ArcGIS REST URL in ArcObjects using C#

texts:
I am trying to create a layer and add it into a MapControl from ArcGIS REST URL. For example I use the following URL http://sampleserver6.arcgisonline.com/arcgis/rest/services/EmergencyFacilities/MapServer

I can do this with the help of this post

However, the code in the post seems to derive a SOAP url from the REST url. So I have few questions regarding this.

Can I create the layer without deriving a SOAP url from the REST url.
Does it always guaranteed that the SOAP url can be derived in this way. For example If the service URL uses the "URL mapping", how will this method work? For instance how can I work with this URL https://maps.stats.govt.nz/wss/service/arcgis1/guest


texts: How accurate is Google Maps when entering lat and long from a survey drawing

texts:
How accurate is inputting the latitude and longitude from a survey pin into google maps to look at the site location?

The problem I am trying to understand is I am about to purchase a property in Palawan / Philippines. The beachfront 2 hectare site was surveyed, by chance I tried to put the Latitude and Longitude from pin 1 into Google Maps to see where it would show the location. Much to my confusion it shows the pin to be 900 metres south of the site I visited and was shown by the owner.

I am from Adelaide Australia, took my info to a surveyor here in Adelaide and he plotted the site and came to the same conclusion. He did suggest a surveyor from the Philippines check this, which I am in the process of doing. As I look at other properties on the same island my simple question is if I put the lat and long from a sketch survey into Google should it be reasonably close?

latitude   10.836883 N
longitude 119.244194 E


texts:
I have data about the basic municipality areas of Mexico in a shapefile.
I also have data about which metro areas they correspond to (csv).  I would like to dissolve them -or so I understand from this article and others.

This picture gives an idea of the number of polygons.  But the white edges are NOT the ones I want to dissolve (they're there to give an idea).

I managed to get the .csv table into the shapefile via other procedures.  And I keep the remaining attributes in this table.

Overall there are about 3000 municipalities, but most of them are linked to no metro areas, and have NULL values.  Only 300 are linked to metro areas, of which there are about 80.  The metro areas with the larges number of municipalities has 75 (that is Mexico City, as it is very big).

I run Vector -> Geoprocessing tools -> Dissolve, and choose the zonas_metr = metro areas field as dissolver.

The problem that I have is that it goes very slowly, and it freezes after more than an hour taking up all RAM.  So I wonder if there's something I can do better or maybe leave it working more time.  I also hope that the NULLed entries are not being merged.  Are they?


texts: Automatically switch scalebar units at certain scale

texts:
I have a scalebar with a element-width of 70 mm. It has three segments and their width is set to be between 50 and 60 mm. Scalebar units are set to [m]. This works fine when it comes to automated adaption to different scales.
But when my scales get smaller (1:30,000, 1:50,000 and similar) the labels of the scalebar change from hundreds of meters to thousands, which reduces readability.
Is there a way to automatically change the scalebar units e.g. from [m] to [km] when crossing a certain threshold of scale?


texts:
How do I reorder the rows of the shapefile in the file? I know how to order them in the attribute table view, but I need to actually switch up the order (by descending id) in the file itself.

It's a little weird but the reason being eventually it will be turned into a geojson and drawn up in a webmapping project but the way the webmapping API works, the draw order is naturally in the order that is in the file, and when icons get very clustered, I'd like items with the lower id to take priority (e.g. be on top / drawn last).

There's probably a way around that in the webmap API but I think it'd just be much simple to reorder the table.


texts:
I got a case where WMS-T is desired, preferably using GeoServer. The input data contains two date fields with start and end date, which describes a timespan.

However, in Geoserver it seems that one can only specify one time field for dimensions. Is there a way to get around this ?

I am thinking that one possible scenario would be to use PostGis for serving data, and create a view/PL-pgSQL that does some of the logic.

Any suggestions are appreciated.


texts: Navbar on ArcGIS JS API causing locator to behave buggy

texts:
I have added a nav bar with combo boxes to extend an ArcGIS JS API application further.  When a location is selected on the map the point does not appear and I receive this error.

Failed to get text for stylesheet 59: No style sheet with given id found

Additionally, when the map is zoomed in, the point will appear, albeit at another intersection or place close to where the user initially clicks.

    div id="divMainContainer" dojotype="dijit.layout.BorderContainer" design="headline"
    gutters="false" style="width: 100%; height: 100%; position: relative">
    <div id="map" class="map" region="center" dojotype="dijit.layout.ContentPane"

     < class='navbar' style="height:450px;">

        <form class="navbar-search pull-left">

            <col>

                                            <select style="width: 80px"> <select id="Select1" dojotype="dijit.form.ComboBox" autocomplete="true" hasdownarrow="true"
                                                style="width: 250px; background-color:Black; font:black;  color: white;" onchange="ValidateRequestType">
                                            </select>

edit: http://jsbin.com/bogosuqari/1/edit?html,css,js,console,output


texts: What values are needed for a SimpleMarkerSymbol?

texts:
My graphicGemetry object contains X and Y values (-9212442.86532579,4421916.145327024) but my image isn't appearing.  I suspect there are other data values I need in order to plot my point but I do not know what they are.

Here is what I have and no errors are being reported:

cppr.App.map.graphics.clear();
var symbol = new esri.symbol.SimpleMarkerSymbol(
esri.symbol.SimpleMarkerSymbol.STYLE_SQUARE,
18, // symbol size
new esri.symbol.SimpleLineSymbol(
    esri.symbol.SimpleLineSymbol.STYLE_SOLID,
    new dojo.Color([0,0,0]),
    1),
    new dojo.Color([255,255,0,0.80])
);
cppr.App.map.graphics.add(new esri.Graphic(graphicGeometry, symbol));

I should also add that I've used this code above when I queried the map service.  What I'm doing now is querying a MSSQL Database to pull back the name.  the database also contains the x & y coordindates so I thought I'd be able to reuse the same code as shown above but apparently there is something missing.


texts: How to calculate and export tasseled caps from MODIS collection in GEE?

texts:
I would like to calculate and export the brightness, greenness and wetness tasseled caps for the MODIS MCD43A4 V006 collection based on the Lobser coefficients, but I have no idea how to do it.

var imageCollection = ee.ImageCollection("MODIS/006/MCD43A4")
    .filter(ee.Filter.calendarRange(5,9,'month'));
   var modbands = imageCollection.select("Nadir_Reflectance_Band1","Nadir_Reflectance_Band2","Nadir_Reflectance_Band3",
   "Nadir_Reflectance_Band4","Nadir_Reflectance_Band5","Nadir_Reflectance_Band6","Nadir_Reflectance_Band7","BRDF_Albedo_Band_Mandatory_Quality_Band1",
   "BRDF_Albedo_Band_Mandatory_Quality_Band2","BRDF_Albedo_Band_Mandatory_Quality_Band3","BRDF_Albedo_Band_Mandatory_Quality_Band4","BRDF_Albedo_Band_Mandatory_Quality_Band5",
   "BRDF_Albedo_Band_Mandatory_Quality_Band6","BRDF_Albedo_Band_Mandatory_Quality_Band7");

var qas = function(image){
  var mask1 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band1").eq(0);
  var mask2 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band2").eq(0);
  var mask3 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band3").eq(0);
  var mask4 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band4").eq(0);
  var mask5 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band5").eq(0);
  var mask6 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band6").eq(0);
  var mask7 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band7").eq(0);

  return image.updateMask(mask1).updateMask(mask2).updateMask(mask3).updateMask(mask4).updateMask(mask5).updateMask(mask6).updateMask(mask7);
};
var masked= modbands.map(qas);

And from this point, I have no idea how to do it. I think I need something like the following, but I t feel I am not any close to the solution.

var tascap= function (image){
  var bands = image.select("Nadir_Reflectance_Band1", "Nadir_Reflectance_Band2", "Nadir_Reflectance_Band3", "Nadir_Reflectance_Band4" ,"Nadir_Reflectance_Band5","Nadir_Reflectance_Band6","Nadir_Reflectance_Band7").multiply(0.0001);

     var bright= ee.Image([0.4395,0.5945,0.2460,0.3918,0.3506,0.2136,0.2678])
      var green= ee.Image([-0.4064,0.5129,-0.2744,-0.2893,0.4882,-0.0036,-0.4169]);
      var wet= ee.Image([0.1147,0.2489,0.2408,0.3132,-0.3122,-0.6416,-0.5087]);

    var br = image.expression(
            '(Band * BRIGHT)',
            {
                'Band':bands,
                'BRIGHT': bright
                }
            );
  var gr = image.expression(
    '(Band * GREEN)',
            {
                'Band':bands,
                'GREEN': green
                }
            );
  var wet = image.expression(
    '(Band * WET)',
            {
                'Band':bands,
                'WET': wet
                }
            );

};

Any help?


texts: Selecting features using expression with 'maximum' aggregates function in QGIS python console

texts:
I am trying to accomplish what should be a simple task: selecting features using an expression in the python console. I am using QGIS 2.18.14 and the Alaska Airports dataset from the QGIS sample data.
I want to select the feature with the highest value in the “ELEV” field.
My code snippet is:

layer = iface.activeLayer()
expr = QgsExpression('"ELEV" = maximum"ELEV"')
selection = layer.getFeatures(QgsFeatureRequest(expr))
ids = [s.id() for s in selection]
layer.setSelectedFeatures(ids)

However, the result is that no features are selected (with no error messages).

If I use this expression : "ELEV"  =  maximum(  "ELEV"  ) in the Select by expression dialog from the GUI, all works as expected, selecting the feature with the maximum elevation value.

Using other expressions in the python console also work successfully, for example:

layer = iface.activeLayer()
expr = QgsExpression('"ELEV" <100')
selection = layer.getFeatures(QgsFeatureRequest(expr))
ids = [s.id() for s in selection]
layer.setSelectedFeatures(ids)

The above code snippet works well, selecting all features with elevation values less than 100.

As I am very much a beginner with Python, I am assuming the problem is with the syntax of the expression within my python code.
To this end I have tried typing the expression in different ways including:

expr = QgsExpression('"ELEV" = maximum"ELEV"')

expr = QgsExpression('"ELEV" = maximum\"ELEV\"')

expr = QgsExpression('"ELEV" = maximum\( "ELEV" \)')

expr = QgsExpression('"ELEV" = maximum\("ELEV"\)')

Which all give the same result- no error message but no features selected.

How do I use the ‘maximum’ aggregates function as an expression with PyQGIS code?


texts: How to get a result close to the QGIS Relief algorithm using only GDAL based commands

texts:
I'm trying to figure a way to output some nice vector features out of a raster DEM in order to have a pleasant visualization on a vector map. I was first trying to play with some algebraic manipulations and combinations between the DEM itself, the hillshade, aspect, slope, TPI, TRI, and many more kernel based approaches...

TPI (Topographic Position Index)

Slope

Aspect

Hillshade

Actual underlying DEM
But I quickly noticed that the possibilities are almost infinite in the space of derived raster that you can compute from a single DEM.
Based on this observation, I was wondering if there are any 'good practices' to render nice elevation maps in a vector form from a raster DEM. But this question is unfortunately too broad and may be opinion based, even if I am intimately convinced that it can be the subject of genuine scientific research, with all its rigor, at the intersection between the disciplines of GIS and imagery visual perception.
Among my researches, I found this nice QGIS algorithm which I never seen before -- this may be because it is quite well hidden; it does not exist in any of the Raster sub-menus --, called "Relief" and which output looks like:

QGIS Relief
Its documentation is here at the time of writing: https://docs.qgis.org/3.16/en/docs/user_manual/processing_algs/qgis/rasterterrainanalysis.html
I usually like to go under the hood and watch what's going on, and I was expecting a bunch of GDAL commands here, because one can see on the resulting image that at least three algorithms have been used (the computation of some contours, a TPI, and a hillshade). So I looked at the console output, but it is surprisingly quiet:
Loading resulting layers
Algorithm 'Relief' finished
QGIS version: 3.16.3-Hannover
QGIS code revision: 94ac9f21b8
Qt version: 5.9.5
GDAL version: 2.4.0
GEOS version: 3.7.1-CAPI-1.11.1 27a5e771
PROJ version: 493
Processing algorithm…
Algorithm 'Relief' starting…
Input parameters:
{
  'AUTO_COLORS' : True,
  'COLORS' : '',
  'INPUT' : '/media/data/gis/DEM.tif',
  'OUTPUT' : 'TEMPORARY_OUTPUT',
  'Z_FACTOR' : 5
}

Execution completed in 1.05 seconds
Results:
{
  'FREQUENCY_DISTRIBUTION': '',
  'OUTPUT': '/tmp/processing_1g3wgU/63301c1a5/OUTPUT.tif'
}

Loading resulting layers
Algorithm 'Relief' finished

I now want to tweak, experiment and configure this algorithm on a headless server where QGIS is not installed (but where GDAL is), hoping I can get some results for my desired vector features out of it.
Hence the first question that came to my mind was; is the QGIS Relief algorithm only based on GDAL commands? Because it looks so to me. And as kindly commented by @user2856 it seems that it's a native QGIS algorithm: https://github.com/qgis/QGIS/blob/master/src/analysis/raster/qgsrelief.cpp
Question
So, how close to the Relief result could we go using a combination of GDAL based tools only?
Other interesting links I came across:
https://www.nature.com/articles/sdata201840.pdf
https://blog.datawrapper.de/shaded-relief-with-gdal-python/
https://grass.osgeo.org/grass78/manuals/addons/r.geomorphon.html


texts: Get Point at distance on a LineString using NetTopologySuite

texts:
I am quite sure there is a function in NTS that actually gets us a point at a distance along a line starting at first vertex. I am using NetTopologySuite version 1.12.1.

In Shapely, the python version of JTS, it is quite simple to use

>>> LineString([(0, 0), (0, 1), (1, 1)]).interpolate(0.75, normalized=True).wkt
'POINT (0.5000000000000000 1.0000000000000000)'

I am trying to do the same thing using NTS but i hightly doubt they could be using a jargon word for this. I tried LinearReferencing namespace, but could not find a suitable method.

Does anyone have idea on this?


texts:
I am unable to load a PostGIS raster layer into QGIS. I have looked through previous posts Loading raster from PostGIS DB into QGIS
and Cannot add the PostGIS raster layer with QGIS but was not able to solve the problem.

I tried both on QGIS Lisboa 1.8.0 and QGIS Dufour 2.0.1. In the Lisboa version I installed the plugin "Load PostGIS Raster to QGIS" 0.5.4 but I get the following error

An error has occurred while executing Python code:

Traceback (most recent call last):
  File "C:/Users/Esteban/.qgis//python/plugins\wktraster\wktRasterPlugin.py", line 90, in callAddLayer
    dlg = DlgAddRasterLayer()
  File "C:/Users/Esteban/.qgis//python/plugins\wktraster\DlgAddRasterLayer.py", line 40, in __init__
    self.updateUIMode()
  File "C:/Users/Esteban/.qgis//python/plugins\wktraster\DlgAddRasterLayer.py", line 109, in updateUIMode
    self.listTables()
  File "C:/Users/Esteban/.qgis//python/plugins\wktraster\DlgAddRasterLayer.py", line 57, in listTables
    tables=conn.listTables(self,connstring) #returns a list of pairs (index, value)
  File "C:/Users/Esteban/.qgis//python/plugins\wktraster\conn.py", line 76, in listTables
    db = GeoDB(host=parmlist[2].split("=")[1],dbname=parmlist[1].split("=")[1],user=parmlist[3].split("=")[1],passwd=parmlist[4].split("=")[1],port=int(parmlist[5].split("=")[1]))
  File "C:/Users/Esteban/.qgis//python/plugins\wktraster\postgis_utils.py", line 134, in __init__
    raise DbError(e)
  File "C:/Users/Esteban/.qgis//python/plugins\wktraster\postgis_utils.py", line 83, in __init__
    self.msg = unicode( error.args[0], 'utf-8')
UnicodeDecodeError: 'utf8' codec can't decode byte 0xe9 in position 23: invalid continuation byte

Python version:
2.7.2 (default, Jun 12 2011, 15:08:59) [MSC v.1500 32 bit (Intel)]

QGIS version:
1.8.0-Lisboa Lisboa, 6416f38

Python path: ['C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python', 'C:/Users/Esteban/.qgis//python', 'C:/Users/Esteban/.qgis//python/plugins', 'C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python/plugins', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\bin\\python27.zip', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\DLLs', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\plat-win', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\lib-tk', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\qgis\\bin', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages\\PIL', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages\\win32', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages\\win32\\lib', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages\\Pythonwin', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\qgis\\python\\plugins\\fTools\\tools']

How to solve the above problem in Lisboa version ?

Is it possible to load a PostGIS raster layer into the Dufour version ?


texts:
I have hundreds of raster files of some agricultural fields. I am using a PyQgis code that takes each raster from a specific folder then calculate the raster's minimum pixel values and maximum pixel values then assign a newly created color ramp to the raster. Here is the codes

from PyQt4.QtCore import *
from PyQt4.QtGui import *
import glob, os

# collect each raster file location
rast_path = r"E:\NDVI"
rasters = glob.glob(os.path.join(rast_path, "*.tif"))

for fileName in rasters:

    #extraction of Information of the raster
    fileInfo = QFileInfo(fileName)
    baseName = fileInfo.baseName()
    layer = QgsRasterLayer(fileName, baseName)
    QgsMapLayerRegistry.instance().addMapLayer(layer)

    #rendering details
    renderer = layer.renderer()
    provider = layer.dataProvider()
    extent = layer.extent()

    #max and min value calculation
    stats = provider.bandStatistics(1, QgsRasterBandStats.All,extent, 0)
    min= stats.minimumValue
    max = stats.maximumValue
    range = max - min
    add = range/2
    interval = min + add
    valueList =[min, interval, max]

    #defining color ramp and assining values
    colDic = {'red':'#ff0000', 'yellow':'#ffff00','green':'#36fc09'}
    lst = [ QgsColorRampShader.ColorRampItem(valueList[0], QColor(colDic['red'])), QgsColorRampShader.ColorRampItem(valueList[1], QColor(colDic['yellow'])), QgsColorRampShader.ColorRampItem(valueList[2], QColor(colDic['green']))]

    #applying color ramp
    myRasterShader = QgsRasterShader()
    myColorRamp = QgsColorRampShader()
    myColorRamp.setColorRampItemList(lst)
    myColorRamp.setColorRampType(QgsColorRampShader.INTERPOLATED)
    myRasterShader.setRasterShaderFunction(myColorRamp)
    myPseudoRenderer = QgsSingleBandPseudoColorRenderer(layer.dataProvider(), layer.type(),  myRasterShader)

    #refreshing layers
    layer.setRenderer(myPseudoRenderer)
    layer.triggerRepaint()

The code automatically imports each raster, calculate max and min values, finally applying the color ramp I would like to have. The raster layers are now visiable on the display and present on the layer panels (as shown below).

I need to export or save the rasters with the exisiting color ramp and the exact same name (what it has now on the layer panel) to a directory ("E:\NDVI\save\"). Is there anyway to do that?

I have tried the following codes by placing inside the loop after the previous codes -

-----Continuation of previous Codes -----------
#refreshing layers
layer.setRenderer(myPseudoRenderer)
layer.triggerRepaint()
------------------new codes ------------------
pa_name, file_name = os.path.split(fileName)
save_folder = "E:/NDVI/save"
save_raster = os.path.join(save_folder, file_name)
width, height = layer.width(), layer.height()
crs = layer.crs().toWkt()

pipe = QgsRasterPipe()
pipe.set(provider.clone())
pipe.set(renderer.clone())
file_writer = QgsRasterFileWriter(save_raster)
file_writer.writeRaster(pipe, width, height, extent, layer.crs())

It is providing the following errors -

Seems like I need different codes for saving the raster.


texts: Filter Features depending on User with GeoServer

texts:
I now use viewparams to filter the features by using sql views in GeoServer which is connected to our postgis-db. This means the parameter comes from the request url.

I want to use keycloak for authentication. And i wonder if its possible to set the filter by user attribute of keycloak. If it is not working with keycloak, is there some other authentification method for this requirement?

Why?
Now the request may filter anything. I want to define the filter for every user.
For example:
UserA have the right to request buildings where city='Berlin' and UserB have the right to request buildings where city='Paris'.
The datasource should be the same and also the layer in GeoServer, because i don't want to have for every city a new layer.

(As Endpoint we use wmts with vector tiles)


texts: How to convert Garmin Waypoints to Google Mercator?

texts:
I have georef some topocards with QGIS and the Google Hybrid Layer.
So, this cards are connected with the Google Mercator System.

Now I have some waypoints taken with Garmin GPSmap62 (in WGS84).

When I try to load both datasets in QGIS, I get a mismatch of around 2000 km (don't know, the topocards are in bulgaria and the waypoints somewhere in the southern atlantic).

I tried to convert the topos to WGS84 and/or the waypoints to Google Mercator.
Either qgis crashed down or nothing happens.

I'm working with Win7, and QGIS 1.7.4

\edit:
First the Topocard. And I try to load up the gpx File, but I'm not sure how gg.
Ah okay, it isn't allowed me to upload pictures, cause I'm a new user. sry


texts: Using custom marker name in Leaflet from external GeoJSON file

texts:
I have the following code (which works great):

<!DOCTYPE html>

<link rel="stylesheet" href="https://unpkg.com/leaflet@1.5.1/dist/leaflet.css" integrity="sha512-xwE/Az9zrjBIphAcBb3F6JVqxf46+CDLwfLMHloNu6KEQCAWi6HcDUbeOfBIptF7tcCzusKFjFw2yuvEpDL9wQ==" crossorigin="" />

However, I wish each marker on the map to use a custom marker name that's specified in the GeoJSON file.

My points.json file contains:

{
"type": "FeatureCollection",
"features": [

    {
        "type": "Feature",
        "geometry": { "type": "Point", "coordinates": [-3.216660, 54.456070] },
        "properties": {
            "Icon": "Icon1",
            "Street": "street one",
            "City": " city here",
            "State": " state here",
            "Zip": " ZIP 1",
            "Misc": " some misc 1"
        }

    },

    {
        "type": "Feature",
        "geometry": { "type": "Point", "coordinates": [-0.551140, 51.202747] },
        "properties": {
            "Icon": "Icon2",
            "Street": "street TWO",
            "City": " city here",
            "State": " state here",
            "Zip": " ZIP 2",
            "Misc": " some misc 2"
        }

    }

]
}

Note the Icon property.

In my HTML code, if I change this:

                    return L.marker(latlng, { icon: Icon4 })

To:

                    return L.marker(latlng, { icon: feature.properties.Icon })

I get the following error in the browser console and no marker is displayed:

TypeError: t.icon.createIcon is not a function

What am I doing wrong here and how do I correct this?


texts: Unable to install ‘sf’ R package on Ubuntu 16.04 LTS

texts:
I'm not able to install sf R package on Ubuntu 16.04 LTS. I tried

install.packages("sf")

as well as

devtools::install_github("r-spatial/sf")

after running

sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable
sudo apt-get update
sudo apt-get install libudunits2-dev libgdal-dev libgeos-dev libproj-dev

Getting the following error message:

  configure: error: gdal-config not found or not executable.

  ERROR: configuration failed for package ‘sf’

Wonder what is the way to install sf?


texts: ArcGis--how to lower a DEM's Z value along a line

texts:
I'm a raw newbie, so please excuse the ignorance.

I'm trying to do the following in ArcGis 10: reduce the Z value of a DEM along a line.  My goal is to "etch" a very shallow line into the DEM surface, so I can then hillshade it and see the etched track in 3d.  I've been messing around with the Interpolate Shape and Add Surface Information (and other) tools, but can't get the right result.

As I see it, the necessary steps probably include these:

Transfer Z information from raster DEM to line.  (use Interpolate Shape or Add Surface Info for this?)
Subtract from the line's Z value to effect the etching or depression. (Adjust 3D Z for this?)
Subtract the difference between DEM raster and line Z-values from the DEM raster.  This is, I think, where I'm getting stuck.  Not sure which tools are appropriate: Raster Calculator, Raster Math (Minus), Intersect, etc.?
Now with "etched" DEM (where Z values have been subtracted along the line), hillshade it.

So basically I'm trying to define a shallow (around 1 meterx1 meter) 3d line and subtract it from the DEM.  I'm sure I'm not going about this the right way--any help would be much appreciated.  Thanks.


texts: Given a shapefile with polygons, how do I convert the polygons to individual lines instead? I know how to do this in QGIS, but I'm looking for a Shapely function which does the same thing.


texts: Selecting polygons with a minimum internal straight line value - ArcGIS 10

texts:
I have a map of polygons representing forest stands. I am trying to sort through them to find those polygons that have a minimum internal distance of 200 m in a straight line (my field sampling will consist of ground truthing the stands, and will require a 200 m line to fall within the stand).

I have tried to find if there is an automated way to do this, which there doesn't appear to be. There was something made thats no longer applicable for arcgis 10 ( http://www.jennessent.com/arcview/longest_lines.htm ).

If it comes to it, I don't mind eyeballing the distance, and drawing a line myself. But I can't seem to figure out a way to draw it so that I can pull geographic coordinates off of it. As in, when I draw two points I somehow convert the drawn points to a feature set of X and Y coordinates. All I can figure out is just doing it by hand, and viewing the values at the bottom right of the screen and recording them by hand.

Thanks in advance, and let me know if you need anything cleared up, not great at explaining my problems (probably why I can't find anything online too).

Ayden


texts:
I've inherited a large somewhat messy database and am attempting to clean it up a bit.
Does anyone know if it is possible to search a database to see what arcGIS map documents use a feature class? I'd like to delete old versions of my feature classes, but don't know if they are in use in a mxd somewhere.
This would also be useful to see what mxds need to be updated because they're using old versions of data.
I'm envisioning a search in ArcCatalog that specifies a feature class and returns a list of mxds and maybe relationships that use said feature class.

Any suggestions appreciated. Thanks in advance!


texts: Replace null pixels from one raster with values from another raster in QGIS

texts:
I have 2 rasters of the same area, one has some null pixels and the other one doesn't and I would like to fill in the nulls in raster 1 by values from raster 2.

How to perform in raster calculator?


texts: What is parentHWnd parameter used in some ArcObjects methods?

texts:
So the parameter parentHWnd appears in a number of ArcObjects methods, such as IExportOperation.ExportFeatureClass.

Does anyone know what the parameter is?

I can't seem to find any explanation. The only value I've ever passed to it is 0. Passing 1 threw an exception that I didn't recognise.

The reason I ask is because I'm trying to figure out if it's possible to hide the "exporting features" progress dialogue that pops up when you call this method, and I'm wondering if this has something to do with it.


texts:
In R, I need to calculate the nearest neighbour index between the relocations of a single animal, for several different individuals. I want to determine whether the relocations of certain individuals are more randomly or regularly distributed than those of other animals.

I have a dataset similar to this:

ID<-  c(a,a,a,a,a,b,b,b,b,b,c,c,c,c,c,c,d,d,d,d,d, etc)
E<-c(167685,167945,167685,153985,167685,158675,167645,167667, etc)
N<-c(9876548,9879248,9876838,9596548,9926548,9878578,9876548,9166548, etc)

In which each "ID" represents a different animal, and E/N points represent locations where the animals were observed. The data frame, consisting of 120 observations of 6 individuals with 20 observations each is called "cor".

I know I can use the "nni" function in the "spatialEco" package to obtain the NNI, the problem is I need to sequentially subset the original data frame in order to calculate NNI for each individual in the sample:

library(spatialEco)
a<-subset(cor, ID=="a")
coordinates(a) <- ~E+N
nni(a, win = "hull")

Being this the output:

 $NNI
 [1] 1.942733
 $z.score
 [1] 8.065562
 $p
 [1] 7.289959e-16
 $expected.mean.distance
 [1] 0.1731123
 $observed.mean.distance
 [1] 0.336311

I would like to create a loop of some kind which would do this set of commands for each individual (ID).
In the output table, each row would summarize the NNI and the rest of parameters fr each individual, like so:

ID          NNI     z.score     etc

a           0.87     2.34       -
b           1.45     3.12       -
c           2.13     4.05       -


texts: Putting features (IFeature) into ifeatureclass with ArcGIS Engine?

texts:
I am trying to get the selected features from axMapControl and put them in a new feature class. I can get the selected features but when i'm trying to put them in a IFeatureClass and then cast it as IFeatureLayer I get the null reference exception.

Can any one help?

List<IFeature> features = GetSelectedFeatures(axMapControl1.get_Layer(0) as IFeatureLayer);
        IPointArray pa = new PointArrayClass();
        for (int i = 0; i < features.Count; i++)
        {
            p.Add(features[i].ShapeCopy as IPoint);
            pa.Add(p[i]);
        }

        IPolyline pl = new PolylineClass();
        pl.FromPoint = pa.Element[0];
        pl.ToPoint = pa.Element[1];
        pl.ToPoint = pa.Element[2];

        pl.SpatialReference = p[0].SpatialReference;

        IWorkspaceFactory workspaceFactory = new ShapefileWorkspaceFactoryClass();
        IWorkspace workspace = workspaceFactory.OpenFromFile(@"C:\", 0);
        IFeatureWorkspace featureWorkspace = workspace as IFeatureWorkspace;

        IFeatureClass FC = CreateFeatureClass("FC", featureWorkspace, pa.get_Element(0).SpatialReference);

        IFeatureCursor fCur = FC.Insert(false);
        IFeatureBuffer buff = FC.CreateFeatureBuffer();

        buff.Shape = pl;
        fCur.InsertFeature(buff);

        fCur.Flush();
        System.Runtime.InteropServices.Marshal.FinalReleaseComObject(fCur);
        ILayer layer =FC as ILayer;
        axMapControl1.AddLayer(layer);

        axMapControl1.Refresh();


texts: I'm trying to calculate and export the static of an image (error:Unknown element type provided)

texts:
var imagem =ee.ImageCollection("LANDSAT/LC08/C01/T1_RT_TOA")

imagem=imagem.filterBounds(region)
    .filterDate ( '2019-01-01' , '2019-12-31' )
    .map(maskClouds)

var mean=imagem.mean();

mean=mean.select(['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7']);

var statistic_combine= mean.reduceRegion({

  reducer: ee.Reducer.min()
                      .combine(ee.Reducer.max(), null, true)
                      .combine(ee.Reducer.mean(), null, true)
                      .combine(ee.Reducer.stdDev(), null, true)
                      .combine(ee.Reducer.median(), null, true)
                      .combine(ee.Reducer.variance(), null, true),

  geometry: region.geometry(),

  scale: 30,

  maxPixels: 1e9
});

Export.table.toDrive({
    collection: statistica_combine,
    description: 'statistica_table',
    folder: 'GOOGLE_GEE_DIR',
    fileFormat: 'CSV'
  });

{Error:Unknown element type provided: object. Expected:  ee.Image,
ee.ImageCollection, ee.FeatureCollection or ee.Element}


texts: Creating View with new geometry field in PostGIS

texts:
I have a table in my data base. This table was created with differents fields as 'x','y' and 'z'.

Now I need to make a view to generate a geometry column, but I found problems, because the created geometry field is not correctly named. I need to force my column "geom geometry" to "geom geometry(PointZ, 25830)".

I did this:

st_force_3dz(st_setsrid(st_makepoint(table.x, table.y, table.z), 25830))

Any solutions?


texts:
I am looking for district boundaries (what OSM puts at level 9) for Tokyo major Wards (level 7 in OSM).

I have found on OSM and affiliated websites some of the district boundaries, but hardly all.

Example : Shinjuku which is a ward , so level 7. OSM has almost no more granular boundaries for the districts contained inside this ward.

Where can I get this data?


texts: How to get accurate mile markers from DOT ShapeFile?

texts:
I'm having trouble getting accurate mile markers through route traversal using PostGIS.

Based on the discussion in Geocode from mile marker/milepost to latitude/longitude?, I downloaded Texas Department of Transportation's official ShapeFile (link, see TxDOT Roadways 2010 at right).

I first imported the ShapeFile into Postgres using shp2pgsql. Then, using a C# program, I used PostGIS's ST_Line_Interpolate_Point function as the meat of my code that traverses each line one mile at a time. Here's the main part of the code:

        public SortedDictionary<int, string> GenerateMileMarkers()
    {
        // make the sorted dictionary
        SortedDictionary<int, string> dict = new SortedDictionary<int,string>();

        // get length of route
        double length = (double)(new NpgsqlCommand(
            "SELECT ST_Length('" + geometry + "')",
            connection)).ExecuteScalar();

        // find each mile marker
        // note that "increment" is 1/69.17, the reciprocal of miles/degree
        for (int i = 0; i <= length/increment; i++)
        {
            double segmentLength = (i * increment)/length;

            string mileMarkerPosition = (string)(new NpgsqlCommand(
                // have to use ST_Linemerge to get from MULTILINESTRING to LINESTRING
                // since ST_Line_Interpolate_Point requires LINESTRING
                "SELECT ST_AsEWKT(ST_Line_Interpolate_Point(ST_Linemerge('" + geometry + "'), " + segmentLength + ")) As foo;",
                connection)).ExecuteScalar();

            dict.Add(i, mileMarkerPosition);
        }

        return dict;
    }

I believe the SHP's reference measurement is a degree of latitude, so I did come calculations and figured that a degree has 69.17 miles, although thanks to Nicklas Avén's comment below, I now realize this is not correct since longitudinal lines are closer together the further north one goes.

Long story short, I generated serially-numbered points for each route that represent mile markers. The problem is I am coming up with significant errors.

I found the error by comparing Google Maps-indicated mile markers with those at a sampling of Interstate highway (IH) intersections or IH intersections with state borders. I'm finding error between ~4% and ~22. Most errors are approx ~15% +/- 1%. Here's an example of the error: where I-30 hits the Texas/Arkansas border is officially mile marker 223 (per Google Maps), but through this method I got 260. That's an error of about 16.6%. Also, where I-30 hits I-635 in Mesquite is about mile marker 56, but I calculated 67. That's a 19.6% error.

This variable error is really concerning to me. Since the errors are inconsistent, even within the same highway, it's not like I can apply a scaling factor to everything.

This is a problem because I am trying to do analysis on a few million events, and these events are expressed as happening at milepost X on highway Y.

Below a user rightfully questioned the projection. It appears the file may in fact have a projection, as it has a PRJ file:

GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

Is it possible to make accurate measurements from this TxDOT-issued SHP?

(NOTE: I used "linear referencing" above, but I don't think this is true linear referencing because I am not correlating a point to a location along a route. The ST_Line_Interpolate_Point PostGIS function just happens to be in the Linear Referencing section of the PostGIS Reference.)


texts: Creating minumum polygon coverage from set of lines in ArcGIS Desktop?

texts:
Does anyone know of a tool or methodology I could use to create a polygon layer enveloping the end nodes of the line layer using ArcMap?

I need a coverage layer for the lines below, I was hoping to find a way to follow the edge of the line (like the illustration by the red line but over the entire dataset) to create a "minimum fit" polygon layer to the lines.


texts: ogr export gdb to PostGIS doesn't contain geometry for display in GIS

texts:
I'm currently trying to export a .gdb containing polygons to PostGIS which other users will need to display using any software like QGis or ArcGIS.
I'm using ogr fom osgeo to export the GDB to PostGIS using the following command:
ogr2ogr -f "PostgreSQL" "PG:[db-connection]" "[.gdb-path]" -nlt PROMOTE_TO_MULTI -overwrite -skipfailure -progress

This manages to upload everything to PostGIS, however the "wkb_geometry" column uses "bytea" as the dafault datatype and when I try to display the geometries using QGIS I'll get the feedback that my table doesn't contain any geometry-type.
So I guess my question would be: Why does this happen and how can I use osgeo (or other tools) to export a gdb to PostGIS as a "proper"  geometry type?
I'm glad for any help/pointers


texts: Attribute Table Corrupted/Altered when Selected SDC Data Exported from ArcMap?

texts:
I am attempting to export selected features from a streets.sdc layer using ArcMap 10.1. After using the 'select by attribute' tool to select the features I want, I then right-click the data layer, choose Data --> Export data, and export the selected features as a new shapefile.

After successfully exporting these features I find that the new attribute table has been altered. For example, many values in the speed limit attribute column have reverted to zero. In addition, other values have been shifted around into different, incorrect attribute columns. For example, the name of a highway may be cut in half, with half of the name remaining in the 'Highway name' attribute column, and the other half migrating to the next attribute column to the right.

How do I prevent the attribute table from being corrupted when I create a new shapefile using selected features?


texts:
I have installed SCP on QGIS 2.18.28 from the Plugins-manage and install plugins. When I tried opening it, it gave me an error of crashed dumped as shown below:

  minidump written to C:
  \users\KLED\AppData\Local\Temp\qgis-20191030-151848-13368-13856-6551c2a73e.dmp


texts: When using ArcGIS Georeferencing tool, why does fourth point show up as yellow crosshair?

texts:
When using the Georeferencing tool of ArcGIS Desktop, why does the fourth point show up as a yellow crosshair?

I have never had this happen before and am worried that it may cause issues.


texts: Using RGB values stored in field of feature class to apply colour component for ArcGIS Pro symbology?

texts:
Symbolizing using RGB values from Excel in ArcGIS Desktop? and Symbolizing vector features using RGB values stored in three fields using ArcPy? are questions that distil to whether it is possible to use RGB values in fields of an attribute table to apply colours direct to a feature layer in ArcMap/ArcPy.  It seems not to be possible so I advised voting for an existing ArcGIS Idea (for ArcGIS Desktop).

My searches of the GUI and web for the equivalent functionality in ArcGIS Pro have also proved fruitless but I am hoping that I am wrong.

Is it possible to use RGB values stored in a field of a feature class to apply the colour component to symbols used to display them in ArcGIS Pro?


texts:
I want to save georeferences readed from input data to output data.
I don't want use copy georeferences because input data is in .h5 format and georeferences was read from metadata.

I use GDAL and C but i accept program in C++ and python

Correct georeferences is EPSG 2180
https://epsg.io/2180
PROJ.4:

+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

In file i can read:

+proj=tmerc +lat_0=0.0 +lon_0=19.0 +k=0.9993 +x_0=500000.00 +y_0=-5300000.00 +ellps=WGS84

It is possible to use:

OSRSetWellKnownGeogCS( hSRS, "+proj=tmerc +lat_0=0.0 +lon_0=19.0 +k=0.9993 +x_0=500000.00 +y_0=-5300000.00 +ellps=WGS84" );

This system does not have UTM zone but if I don't use this it not working:

OSRSetUTM( hSRS, 11, TRUE );

I want to save this georeferences:

+proj=tmerc +lat_0=0.0 +lon_0=19.0 +k=0.9993 +x_0=500000.00 +y_0=-5300000.00 +ellps=WGS84

Convert to this:

PROJCS["ETRS89 / Poland CS92",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",19],PARAMETER["scale_factor",0.9993],PARAMETER["false_easting",500000],PARAMETER["false_northing",-5300000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","2180"]]

To file using:

GDALSetProjection( hDstDS, georeferences);

Is it possible or how to convert it?

GDAL Info:

bash-4.2$ gdalinfo 201804010000_acc0010_grs.h5
Driver: HDF5/Hierarchical Data Format Release 5
Files: 201804010000_acc0010_grs.h5
Size is 512, 512
Metadata:
  dataset00_data0_data_CLASS=IMAGE
  dataset00_data0_data_IMAGE_VERSION=1.2
  dataset00_data0_what_quantity=ACRR
  dataset00_quality0_data_CLASS=IMAGE
  dataset00_quality0_data_IMAGE_VERSION=1.2
  dataset00_quality0_what_quantity=QIND
  dataset00_what_prodpar=000
  how_software=grs
  how_source=GAU+ GAU_QI+ RAD+ RAD_QI+ SAT+ SAT_QI+
  how_source_gau=376
  how_source_rad=BRZ- GDA- LEG- PAS- POZ- RAM- RZE- SWI- czbrd- czska- skjav+ skkoj+ skkub+ sklaz+ UKLL+ gedrs- gepro- geros- ltlau- ltvln-
  how_sw_args=compr=ZLIB, int=OK, grs=QSP, grs_qi=AVG_W
  how_sw_version=2.0
  what_acclength=10
  what_date=20180401
  what_nodata=-999
  what_object=COMP
  what_product=SURF
  what_time=000000
  what_type=GRS
  what_version=RainDIM 1.1
  where_projdef=+proj=tmerc +lat_0=0.0 +lon_0=19.0 +k=0.9993 +x_0=500000.00 +y_0=-5300000.00 +ellps=WGS84
  where_xllcorner=50000
  where_xscale=1000
  where_xsize=900
  where_yllcorner=30000
  where_yscale=1000
  where_ysize=800
Subdatasets:
  SUBDATASET_1_NAME=HDF5:"201804010000_acc0010_grs.h5"://dataset00/data0/data
  SUBDATASET_1_DESC=[800x900] //dataset00/data0/data (64-bit floating-point)
  SUBDATASET_2_NAME=HDF5:"201804010000_acc0010_grs.h5"://dataset00/quality0/data
  SUBDATASET_2_DESC=[800x900] //dataset00/quality0/data (64-bit floating-point)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)


texts: Recommended approach to integrate GeoNetwork with SQL Server database

texts:
I am new to GeoNetwork. But now my organization wants to invest in it to see how we can make the best use of it for publishing, sharing and distributing metadata.

We currently have all our GIS data (mostly feature classes) stored in ArcSDE databases on SQL Server 2008 R2. Only brief information about the data has been captured, therefore no standardization has been pursued.

My question is if we get a GeoNetwork catalogue set up, what is the best way to integrate it with the existing SQL Server databases. Also how does each metadata record in GeoNetwork get associated with the actual dataset in SQL Server? Can we set up GeoNetwork for it to automatically extract simple metadata of our specified from the SQL Server and push it to the GeoNetwork catalogue?


texts: Mapnik: GeoTIFF flipped (mirrored vertically)

texts:
I generated a GeoTiff using the QGIS HeatMap plugin and want to render it on top of some OSM data. The problem is that the raster gets rendered upside-down in the output. Can someone point me in the right direction?

mapnik xml:

<?xml version="1.0" encoding="utf-8"?>
<Map background-color="#f2efe9" srs="+proj=latlong +datum=WGS84">
  <! ... Styles ... -->
  <Layer name="highways" status="on" srs="+proj=latlong +datum=WGS84">
    <StyleName>highways</StyleName>
    <Datasource>
      <Parameter name="type">osm</Parameter>
      <Parameter name="file">data.osm</Parameter>
    </Datasource>
  </Layer>
  <Layer name="heatmap" status="on" srs="+proj=longlat +datum=WGS84 +no_defs">
    <StyleName>heatmap</StyleName>
    <Datasource>
      <Parameter name="type">gdal</Parameter>
      <Parameter name="file">test.tif</Parameter>
      <Parameter name="format">tiff</Parameter>
    </Datasource>
  </Layer>
</Map>

$ gdalinfo test.tif

Driver: GTiff/GeoTIFF
Files: test.tif
Size is 709, 500
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (...)
Pixel Size = (0.000849347874585,0.000849347874585)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
...
Band 1 Block=709x2 Type=Float32, ColorInterp=Gray
  NoData Value=-9999

$ gdalsrsinfo test.tif

PROJ.4 : '+proj=longlat +datum=WGS84 +no_defs '

...

render.py

import mapnik
m = mapnik.Map(1024, 1024)
mapnik.load_map(m, "mapnik_style.xml")
bbox = m.layers[1].envelope()
m.zoom_to_box(bbox)
mapnik.render_to_file(m, "map.png", "png")

edit: I used gdal & numpy to reverse the pixels and it rendered correctly. But then that raster displays upside-down in qgis.


texts:
I'm trying to use getis ord on the new york taxi dataset to identify the hotspots based on the number of pickups.

I want to model the data into cells such that that score of each cell is calculated based on its queen neighborhood.

I found a package called spdep in R that has a func to calculate getis ord but I'm still not able to use it to solve my problem as I'm still new to R.

I would appreicate if you can give me some hints about creating the neighborhood or a guiding example.


texts: Partitioning value (100 tonnes/ha) of underlying shapefile to multiple features in same shapefile in areas that overlap using ArcGIS for Desktop?

texts:
I have tried searching for an answer but can't think of a way to generalize my question to other scenarios. Looking for answers using ArcMAP 10.1.

I currently have two shapefiles. One contains different values for the quantity of dry matter (DM) per hectare (e.g. one patch of 100 metric tonnes (t) of dm per hectare, another of 50, another 70 and so on...). The other shapefile has 1 km pixels scattered over different locations, and each of these pixels has its own time and date stamp attached to it.

If two pixels overlap I want to distribute the percentage of that 100 t per hectare that between the two (or more) overlapping areas proportionatly to the area.

For example If two 1 ha pixels (ha to make it simple) have a 50% overlap, both features would have 25 t for the overlapping area (in reality we would now have 4 pixels); and if their were four 4 ha pixels with 50% overlap in the same space it will be 12.5 t (with 8 pixels)etc. (example picture below)

I would like to apply the above workflow in a model in ArcGIS 10.1 to a big dataset. I am happy to convert the shapefiles to rasters etc and have full access to the different tools.


texts: Lot of roads are missing from pgRouting pgr_drivingDistance result

texts:
I started working with pgr_drivingDistance. Lot of edges (marked as yellow) are missing from result.

Those seem to be cases, where routing algorithm has approached to same edge from two directions.

Is such result normal and expected?

I can fix those missing edges with separate queries, but at first I would like to be sure, that I haven't messed something up myself.

EDIT: I have added some post processing now to fix those missing edges.

I have used different software versions (Win platform): latest versions (Postgres 11 + PostGIS 2.5 + pgRouting 2.6) and Postgres 9.6 + PostGIS 2.4 + pgRouting 2.3. Results are identical.


texts:
I am migrating a plugin from QGIS 2 to QGIS 3. The following code is QGIS 2:
    # Define a mapscale-dependent rectangle about the pick point.
    pntGeom = QgsGeometry.fromPoint(point)
    pntBuf = pntGeom.buffer(canvas.mapUnitsPerPixel()*4, 0)
    rect = pntBuf.boundingBox()

    # Convert the rectangle coordinates to the layer CRS.
    layerRect = canvas.mapRenderer().mapToLayerCoordinates(layer, rect)

    # Select the features that are in the rectangle and return their fids.
    layer.removeSelection()
    layer.select(layerRect, True)
    fids= layer.selectedFeaturesIds()

where canvas is a self.iface.mapCanvas() object and layer is a QgsVectorLayer object.
I need this to work in QGIS 3.  The mapRenderer() object is obsolete and it appears the mapToLayerCoordinates() function is a member of the QgsMapSettings and requires that the layer argument "must have type QgsMapSettings" so my replacement code:
    layerRect = QgsMapSettings.mapToLayerCoordinates(layer, rect)

fails with

TypeError: QgsMapSettings.mapToLayerCoordinates(): arguments did not
match any overloaded call:
overload 1: first argument of unbound method must have type 'QgsMapSettings'
overload 2: first argument of unbound method must have type 'QgsMapSettings'

The ultimate goal is to return any feature ids that are within the bounding box.
How do I resolve this conflict?


texts: How to create a sub dataframe from a dataframe using specific columns

texts:
I have a dataframe that has more fields/columns than I want, so I am trying to create a sub dataframe that only has certain fields/columns. The following throws a

  TypeError: '(['uniqueid'], ['PARCEL_ID'], ['LAND_USE_T'],
  ['PROP_IND_T'],['Bldg_Sq_Ft'],
  ['GEOID'], ['CensusPop'], ['HU_Pop'])' is an invalid key

feat_df = gpd.read_file(state_gdb, layer=join_feat)
new_df = feat_df[['uniqueid'], ['PARCEL_ID'], ['LAND_USE_T'], ['PROP_IND_T'], ['Bldg_Sq_Ft'], ['GEOID'], ['CensusPop'], ['HU_Pop']]

What is the easiest way to create a sub dataframe with only the specific columns I want?


texts: Can I use predefined list of values from Postgis when I draw lines in QGIS?

texts:
Lets say that we have 2 tables in Postgis:

Table 1 (PIPE):
-ID
-pipe_type_id
-geometry

Table 2 (PIPE_TYPE):
-ID
-type

foreign key PIPE.pipe_type_id - PIPE_TYPE.ID

Table PIPE is empty. Table PIPE_TYPE is filled with values.

I need to draw a line in QGIS directly to PIPE.geometry column. After the line is drawn I need a popup window to appear with an option to choose PIPE_TYPE.type value for that line. After submitting, correct value of PIPE_TYPE.ID should be stored in PIPE.pipe_type_id column.

I did this with Bentley Map using their Map Administrator and Oracle DB, but I want to try this with QGIS + Postgis. Is it possible?


texts:
I am using OSM/Google maps and WMS Layers. Now I want to know how my WMS will overlay on this map. I am using the following code for that:

function init()
             {
                map = new OpenLayers.Map('map');

                layer = new OpenLayers.Layer.OSM( "OpenStreet Map");
                map.addLayer(layer);
                layer = new OpenLayers.Layer.Google( "Google" );
                map.addLayer(layer);

                var jpl_wms = new OpenLayers.Layer.WMS("Administrative Layer",
                "http://localhost:8080/geoserver/test/wms",
                {layers: "maharashtra_administrative",transparent: true,visibility: false},{isBaseLayer:false});
                map.addLayer(jpl_wms);
                map.addControl(new OpenLayers.Control.LayerSwitcher());
                map.setCenter(new OpenLayers.LonLat(8461018,2120810),6.75);
                map.addControl(new OpenLayers.Control.MousePosition({}));
                map.addControl(new OpenLayers.Control.OverviewMap());
                map.addControl(new OpenLayers.Control.ScaleLine());
            }

Can anybody guide me for the same as this overlay does not appear in the map.


texts: Self-intersecting after exporting .shp from Personal GeoDatabase

texts:
I have a Feature Dataset in a Personal Geodatabase, with resolution = 0.00000001 and tolerance = 0.000000005.

I am exporting the Polygon Feature Classes to shapefiles, after setting Resolution = 0.0001 in the Environment Settings of the Export tool (generalising them).

The resulting shapefiles have too many self-intersecting polygons (reported by the Check Geometry Tool), while the FCs in the GeoDB were Ok!

Any suggestions?

I know these settings in the GeoDB look ridiculous, but it's not my choice...


texts: Grid from SAGA natural neighbor results in black grid

texts:
I'm new in QGIS and I'm following tutorials mainly on YouTube so I created a path in Google Earth then exported to kml which is converted to txt via GPS visualizer.
When added delimited layer it shows just fine but when used SAGA raster creation and Natural neighbor it run without errors but the grid that comes as result is black.

How do I convert the map to 3d contour to get terrain?
I tried to change options and settings but didn't help.


texts:
I'm using Leaflet right now and I have issues on ordering a list of FeatureGroups.
Each feature group is either an imageOverlay or a list of Features

I tried using bringToFront / bringToBack and it's working well with FeatureGroups containing a layer of imageOverlay but not with layers containing Features.

What am I missing ?

displayLayer(layer): void {
    const fg: L.FeatureGroup = L.geoJSON();

    if (layer.features !== null && layer.features.length > 0) {
      for (const sc of layer.features) {
        const feature = {};
        feature[`type`] = 'Feature';
        feature[`geometry`] = JSON.parse(sc.geometrie);

        const l: L.GeoJSON = L.geoJSON(JSON.parse(JSON.stringify(feature)));
        fg.addLayer(l);
      }
    } else {
      const bounds = JSON.parse(layer.bounds)
      const c1 = L.latLng(bounds.coordinates[0][1], bounds.coordinates[0][0])
      const c2 = L.latLng(bounds.coordinates[1][1], bounds.coordinates[1][0])
      const imgCalque = L.imageOverlay(layer.url, L.latLngBounds(c1, c2));

      fg.addLayer(imgCalque);
    }

    if (fg.hasLayer) {
      this.carte.addLayer(fg);
      this.visibleLayers.set(layer.id, fg);
    }
  }

bringToFront(layer): void {
    const fg = this.visibleLayers.get(layer.id);
    if (fg !== undefined) {
      fg.bringToFront();
    }
  }

  bringToBack(layer): void {
    const fg = this.visibleLayers.get(layer.id);
    if (fg !== undefined) {
      fg.bringToBack();
    }
  }


texts: How do you parse an HTML formatted attribute value to seperate columns in a database table?

texts:
All,

I have an attribute field formatted using HTML all ready to go for google earth markers.

The KML in question is here:

http://dl.dropbox.com/u/95231813/state_pts_MI.kml

Here is the value:

    Facility type:

    Geoid: 261450002001020

    County: Saginaw County

    2010 Correctional Population: 18

(i tried to paste the actual HTML but couldn't figure out how to display it as text rather than formatted html!!!)

I've love to have seperate attribute fields for the columns within that string:

Facility Type
County
2010 Correctional Population

How can this be done?

QGIS? LibreOffice? Python?

Thanks!!!

-mb


texts: Assign Custom Band Combination To Image Catalogs In ArcGIS

texts:
I have an image catalog composed of 4 band rasters.  I would like to use the same image catalog, and set custom renderers for a true color and false color combination and save each of those catalogs off as a layer file.  In addition I would like to turn off the default stretch for each catalog.  The code I am using is below (where lyr is a reference to a raster catalog):

                IRasterCatalogLayer rcl = lyr as IRasterCatalogLayer;

                IRasterRGBRenderer pRGB = new RasterRGBRendererClass();
                IRasterStretch pRasStretch = pRGB as IRasterStretch;
                rcl.Renderer = pRGB as IRasterRenderer;

                pRasStretch.StretchType = esriRasterStretchTypesEnum.esriRasterStretch_NONE;
                int r, g, b;
                pRGB.QueryBandIndices(out r, out g, out b);
                pRGB.UseRedBand = true;
                pRGB.RedBandIndex = 3;
                pRGB.UseGreenBand = true;
                pRGB.GreenBandIndex = 0;
                pRGB.UseBlueBand = true;
                pRGB.BlueBandIndex = 1;

                int r1, g1, b1;
                pRGB.QueryBandIndices(out r1, out g1, out b1);

                rcl.Renderer.Update();

The stretch is being applied without problem, but I cannot change the band combinations.  Thoughts?  Is there a better way to do this without resorting to ArcObjects to create layer files.


texts:
I made a tool in ArcGIS Python Add-Ins for my colleagues and right now wrote an update in main script. It would be easy to tell everyone to delete old add-in and install new one, but recently found the path to the initial script of a tool which is C:\Users\%username%\AppData\Local\ESRI\Desktop10.2\AssemblyCache\folder with a long name\file.py.
I tried to make a script that could make a backup of file.py to file_backup.py, take a new file and paste it in found folder.
Script works well, the new file with new updates is there, but when I launch ArcMap there are no any changes. I can see that the file is the same that I used before.

Is it technically possible to perform this operation using the change of files?

The script is:

import os
import shutil
loc = os.environ.get('LOCALAPPDATA')
filename = 'Pointer_addin.py'

for x in os.walk(loc):
    if filename in str(x) :
        print os.path.join(x[0],filename)
        path = x[0]
        filepath = os.path.join(x[0],filename)

os.rename(filepath,filepath.rstrip('.py') + '_backup.py')
newp=r'D:\Pointer\Install\Pointer_addin.py'
shutil.copy(newp,path)


texts: Dissolve polygons in Spatialite - New Table has no Geometry

texts:
I have followed some examples for dissolving polygons in Spatialite and successfully created my simpler polygons:

SELECT
  "col1"
  sum("total") AS "total",
  st_union("GEOMETRY") AS "GEOMETRY"
FROM "table"
GROUP BY "col1";

I can preview these new shapes in Spatialite (Map Preview, Check Geometries), but they're not showing up as 'geometry' tables, with the little world icon, and I can't view them in QGIS.

I've seen examples that say I need to use something like this:

SELECT AddGeometryColumn('test_geom', 'the_geom',
  4326, 'POINT', 'XY');

Can I convert the existing "geometry" column without adding a "geometry1" column?

Also, some of merged shapes are polygons and some are multipolygons due to minor errors in the original file coverage. I hope to edit the errors in the dissolve in QGIS, what new geometry type is best use for the table?


texts: ogr2ogr PDF unable to open datasource but PDF is recognized fine in GDALInfo

texts:
I am attempting to use ogr2ogr to convert a PDF to a SHP through the OSGeo4W shell, part of the QGIS installation (3.6 i think) on windows 10.  I have used this before successfully but been a while so unsure why it isn't working.  I have tried using 2 PDFs (both georeferenced - i can see the latlong).  Both PDFs are recognized in GDALInfo and all the info is provided easily.
When run in ogr2ogr i use the simple following line:
ogr2ogr dest.shp source.pdf

but get this error
Unable to open datasource `source.pdf' with the following drivers.
  -> `JP2ECW'
  -> `OCI'
  -> `SOSI'
  -> `PCIDSK'
  -> `netCDF'
  -> `JP2OpenJPEG'
  -> `PDF'

It has the PDF driver (which works fine in GDALInfo) so why not here?
Thanks
Joe
To add to this:
I have just researched this a bit further and found that both PDFs were made in ArcMap.  When exporting PDFs I know the georeference option needs to be checked for this to work.  I also found that only the middle option "Export PDF Layers Only" works in ogr2ogr.  If either the  flattened version "None" and the "Export PDF Layers and Feature Attributes" are chosen then the ogr2ogr conversion does not work.  I cant see a way around this
Examples PDFs here if you want to see.  Only GeoPDF2.pdf seems to work in OGR2OGR.
https://we.tl/t-xBBRCJg9TH


texts: Deleting symbols in QGIS symbol selector?

texts:
How do I delete symbols in the symbol selector QGIS?

I am trying to clear out a few old symbols from the library and was wondering how id do it.


texts: Testing if two GDAL datasets are in same projection

texts:
I'm attemping to write a function to test if a list of datasets are in the same projection.  Naively, I thought I could test the text of WKTs, but that does not work in the following case:
'PROJCS["NAD_1983_UTM_Zone_14N",GEOGCS["NAD83",DATUM["North_American_Datum_1983", ...'
'PROJCS["NAD83 / UTM zone 14N",GEOGCS["NAD83",DATUM["North_American_Datum_1983", ...'
Even though the projections are identical, and the rest of the WKT is identical, the tags leading them are not which makes me wonder if there is a formal mechanism for testing the equality of projection systems.
I am programming with the GDAL extensions in Python.


texts: I had a raster layer in my map document that suddenly became invisible in the map even though it is still in my table of contents. There is now a red exclamation mark next to it in the table of contents though. Even when I re-add the layer from my catalog the layer is invisible and keeps that red exclamation mark. I can't find any explanation of this or how to fix it. What does it mean and what can I do about it?


texts: Calculating distance between multiple points using ArcMap?

texts:
The shapefile I created contains 5 points. Their FID's are (0,1,2,3,4).
I need to calculate the distance between all of them. The output table should contain the following info -

points to compare      distance
        0-1           float value
        0-2           float value
        0-3           float value
        0-4           float value
        1-2           float value
        1-3           float value
        1-4           float value
        2-3           float value
        2-4           float value
        3-4           float value

I have found some ways of calculating the distance from one particular point to other multiple points. But trying to calculate the distance between all of them.
Is there any way to do that?

I am using ArcMap 10.6 and ArcPy 2.7.


texts: How to Create database using OGC compliant geometries and Custom type Geometries?

texts:
I am interested in storing geometry attributes in a PostGis database, so that it would be OGC compliant, along with OGC compliant geometry types like point,line, linestring, etc. I also wanted to store some custom type geometry types like ellipse,circle,arrow,some symbols (also combinations of these types) etc, these custom geometry types have different attributes(entities).

Is it possible to store all these geometry in a Generic OGC compliant DataBase?

I'm using pdf file "URL:http://portal.opengeospatial.org/files/?artifact_id=829"(99-049_OpenGIS_Simple_Features_Specification_For_SQL_Rev_1.1.pdf) as a reference, and as per page 2-20, it appears like I can represent any kind of geometry using GEOMETRY_COLUMNS, SPATIAL_REFERENCE_SYSTEMS, FeatureTables and GEOMETRY_COLUMNS. But a lot of the geometries I have are custom created - like thick curved arrows with control points to allow a user to adjust the curvature of the arrow, and a text label at the center of the arrow, polygons with an X mark at the center of every edge and control points at certain points of the polygon to allow the user to resize it etc.

like i have n no:of item(custom_geometries like,arrow,symbols etc)each item has unique attributes.
should i create multiple feature table for each item?


texts: Animating ArcGIS API for JavaScript Div Map Element?

texts:
I have an ArcGIS JS map with a time slider. I've created a Div element that shows the local time in two minute intervals. I want it to follow a shapefile as it moves across the screen diagonally, down and to the right.

The element is in the right place at first, but on the second time change it disappears. I'm not sure what I'm missing. I've defined the Div position in the css, and I want to change it with each time change. Here are the relevant sections of my code.

 #timePanel {
        position: absolute;
      top: 225px;
      left: 300px;
      z-index: 50;
  }

function timeMove() {
             var elem = document.getElementById("timePanel");
             var pos = 225;
             var id = setInterval(frame, 2000);
             function frame() {
                  pos++;
                  elem.style.top = pos + '25px';
                  elem.style.left = pos + '100px';
                }
             }

      timeSlider.on("time-extent-change", function(evt) {
        var startValString = evt.startTime.toLocaleTimeString(navigator.language, {hour: '2-digit', minute:'2-digit', timeZoneName: 'short'})
       ;
        dom.byId("daterange1").innerHTML = "<b>" + startValString + "<\/b>"
           dom.byId("daterange2").innerHTML = "<b>" + startValString + "<\/b>"

           timeMove();

        });


texts: Select different columns from each feature layer to be count values

texts:
How to count values (occurences of values) in each feature layer using a loop?
The problem is that the column names are not the same for all layers.
I have a table that has columns:
TheTmpQueryName tablename  query    duplicated
q1                 NTA      query1     0
q2                 IT       query2     1
q3                 INTS     query3     1

I want for each row that has in duplicated = 1 to count values for a specific column.
First do filtering of the rows that the count values will be applied
fc = r'C:\data.gdb\tableingdb'
field = 'duplicated'
with arcpy.da.SearchCursor(fc,field) as cursor:
     for row in cursor:
         if row[0]==1:
          #here how to make a new column that will count how many times occurs?

based on this?
dict1 = ({'NTA' :  ['AIN'] ,
'IT' : ['CIN'],
'INTS' : ['KAD'], 'IOT' : ['EK']

Where NTA, IT, INTS are the layers and AIN,CIN, EK the columns on which the count should be based?
Example:
NTA(layer)
AIN   COUNT
SSA   21
SSC   1
SDD   98

IT(layer)
CIN COUNT
SRS   23
SCV   14
SDE   12

and so on.


texts: How to display a map from a script in Grass Gis 7.6.1

texts:
I am working with Grass Gis and spyder (IDE) and I am trying to display a map on "Spyder console" from the command line of spyder. I have a raster map call "elevation" and i tried this script:

import grass.script as gscript
def main():

    gscript.run_command('d.rast', map='elevation')

if __name__ == '__main__':
    main()

but the CONSOLE return this:

ERROR: Neither MONITOR (managed by d.mon command) nor
       GRASS_RENDER_IMMEDIATE (used for direct rendering) defined

I dont know what is GRASS_RENDER_INMEDIATE. If I can display in WxGui and not in "spyder console" maybe I can work too.