texts:
I am trying to modify a dataset by replacing the geometry of each feature by its buffer. I try to do it directly, in update mode (not very secure I know). Here is what I tried:
ogr2ogr -update test.gpkg test.gpkg -dialect sqlite -sql "UPDATE test SET geom = ST_buffer(geom,100000)"

but it does not work. No error is returned but the dataset has then empty geometries.
Any idea about something I missed ?
The geometry column is geom and the CRS is in meters, so that the distance 100000 is realistic.


texts:
I need to calculate the quartile value (e.g. 75%) of the pixels of a raster layer within a polygon file. I try to do that using QGIS, and found that the zonal statistics might be able to do that. But I only find the below options with the zonal statistics: Count; Sum; Mean; Median; StDev; Min; Max; Range; Minority; Majority; Variety.

Any solution to calculate the quartile value in QGIS?


texts: I have a simple question: how do you rename your Assets in Google Earth Engine (I'm not talking about the files within the Assets)? For some reason I cannot rename it and I have a hard time creating a new one and delete the older folder!


texts:
I am building a map application for a northerly area and would love to bring in a feed of northern lights prediction/probability (like below).
http://www.gi.alaska.edu/AuroraForecast

Does anyone know of a service for space weather that I could digest into leaflet or something?


texts: US One Address locator style missing city, state and zip code fields?

texts:
I am using the file:
C:\Program Files\ArcGIS\Desktop10.0\Locators\US One Address.lot

To do address resolution on a parcel / polygon layer.  I'm trying to understand why the file doesn't contain city, state or zip code fields.  It only contains:

{ HouseNum, PreDir, PreType, StreetName, StreetType, SufDir }.

MWrenn suggested I try the locator named "US One Address with Zone" and map city to zone.  This technique works, but the results omit state and zip code.  How can I get city, state, and zip code as a "US One Address" locator?

ArcEngine 10, C#, VS2010

Also posted here:
http://forums.arcgis.com/threads/25632-Why-doesn-t-US-One-Address-locator-style-have-city-state-or-zip-code-fields?p=84515

Edit:
Sean made the point that

  "Having an address and ZIP code is
  sufficient to geocode the address"

which I don't take issue with, but maybe I didn't properly explain why I want City, State, Zip...

If the address is not found in the search and address candidates are found, I open a window which allows the user to select which address they want.  Unsurprisingly, the client wants a regular looking address with house number street name, city, state, zip.

Here is how I fetch addresses:

IArray resultsArray = addressCandidates.FindAddressCandidates(addressProperties);
IFields candidateFields = addressCandidates.CandidateFields;

for (int candidateIndex = 0; candidateIndex < resultsArray.Count; candidateIndex++)
{
ClosestAddressResults thisResult = new ClosestAddressResults();
candidatePropertySet = resultsArray.get_Element(candidateIndex) as IPropertySet2;

List<string> dataAvailableToDisplay = new List<string>();
for (int fieldIndex = 0; fieldIndex < candidateFields.FieldCount; fieldIndex++)
{
  addressField = candidateFields.get_Field(fieldIndex);
  dataAvailableToDisplay.add(addressField.AliasName + " " + candidatePropertySet.GetProperty(addressField.Name));
}

The list looks like this:

Shape System.__ComObject
Score 66
HouseNum 104
PreDir
PreType
StreetName 1ST
StreetType ST
SufDir
Zone 56441
Match_addr 104 1ST ST, 56441

In this case I mapped Zip to Zone.  If I had mapped City to Zone, I would get the city name back.  Because I am allowing the user to select an address in a window I want to display the addresses with all the data, but since only Zone comes back, I'm having a difficult time figuring out how to get it.

Right now my workaround is may City to Zone.  Then for each address match, get the IPoint, for each IPoint, run it through a find closest address function, with Zone mapped to Zip, grab the Zip code, and concatenate the results.  And of course I could do this for state too.  The problem with my workaround even though technically it works, is that its pretty slow.  So any ideas whatsoever are appreciated.


texts:
I have the latitude and the longitude of view of map ( latitude and longitude of north west and the latitude and longitude of southeast ) and a geojson files with MultiPolygon. (43mo)

In PHP i want to detect if a polygon is in view.

test they 2 point (NW & SE) if is in polygons is not a good way because polygons is imprevisible https://upload.wikimedia.org/wikipedia/commons/thumb/8/81/Rectangle_definition.svg/2000px-Rectangle_definition.svg.png in this case the result found is the corner NW and SE but the center and ne and sw not found

I need to know what are in view

example with this view:
I test NW and SE (of view) if is in polygon on center -> return false
but i see the polygon on center

how i can know all polygons in actual view ?


texts: Alternative to Farm Service Agency (FSA) CLU dataset?

texts: Is there an alternative dataset that can I find that is similar to the Farm Service Agency (FSA) Common Land Unit (CLU) ?


texts: ArcMap ModelBuilder Date query in line substition

texts:
I have successfully created the Minimum date and time variable as well as the maximum date and time variable(date datatype). I am only looking 30 seconds in either direction from the record I am querying. My problem is when I try to do an in line substation for the select layer by attribute it will not accept [Date_stamp] >= %#Min_Query#% it says there was a syntax error. The field I am querying is a date and time field and the variable I am using is a date time variable. Both appear to always have an associated time. I added a picture to show exact formats.

I am using ArcMap 10.6.


texts:
I have 4000+ raw .tif images. For each image, I have the four corner coordinates in Irish Grid System. I'd like to know how I can use GeoTools to georeference each .tif image using the coordinate values I have.
After following some search, I noted that GeoTools have functionalities for my specific need. But I'm still unable to find a concrete example.


texts: I'm using a dragbox to capture coordinates from an openlayers 3 map.  Is there a way to keep the dragbox visible once you release the mouse buttom and have it disappear only once you start a new dragbox?


texts: Is it scientifically valid to compare landcover/use generated from different types of spatial resolution products

texts:
This question is mainly for the validity of the comparison of the classification product of different types of spatial resolution.

Is it scientifically valid to compare very high-resolution image analysis product with that of moderate resolution. Say I have RapidEye image classification(8 class) of 2011 and forest loss pixels(just 3 class:forest loss, no forest loss, not applicable) of 2010 from landsat image products. Now I want to see what forest loss pixel of 2010 is converted to what in 2011.

Is it scientifically valid or not to compare product of different spatial resolution?


texts: ESRI Javascript - Dynamically Create Points from Non Spatialized JSON Using Lat/Longs

texts:
I am working on application that needs to make use of a non spatial data set. The data is maintained by the group responsible for the installation of the asset item and thus it makes sense to pull from them rather than keep track of a second data set for spatial solutions. I am able to send a request to an API that queries the database and returns an array of JSON. The objects look as follows:

[
    {
        attribute01:value01,
        attribute02:value02,
        ......
        Lattitude:lattitude,
        Longitude:longitude1
    },
    {
        attribute01:value01,
        attribute02:value02,
        ......
        Lattitude:lattitude,
        Longitude:longitude1
    }
]

What I would like to do is dynamically create the points from this Lat/Long data as well as maintain the other attributes in the response object. I do not require just the points mapped, but rather  I would like to be able to interact with the points in a same manner as I would a feature layer. Is this possible?


texts: QGIS DBManager error when selecting a feature

texts:
Possible Duplicate:
  Why does my PostGIS query fail with “current transaction is aborted”

having loaded a PostGIS layer on the QGIS canvas using the DBManager plugin, I am unable to select any ot its features using the Single Feature tool. I receive the following error:

1 cursor states lost.
SQL: CLOSE qgisf13
Result: 7 (ERROR:  current transaction is aborted, commands ignored until end of transaction block
)

The other selection tools give similar errors.
Any suggestions?
Thank you
George


texts: Varying size buffer along line with PostGIS

texts:
I want to create a buffer along a line where the buffer width corresponds to a fraction of the length of the line, something like that:

To do that my idea is to compute for each node of the line (except for the first and last node) the bisector angle between the previous and the next segment, example for the node, bisector angle in red #3:

At the moment I'm just trying to create one side (a line) of the buffer:

Here's my attempt:

WITH POINT_LINE AS
(
WITH SUB AS
(
                -- Extract all the segments and calculate for each segment the cumulative length
        WITH segments AS (
        SELECT gid, (pt).path[1]-1 as path, ST_MakeLine(lag((pt).geom, 1, NULL) OVER (PARTITION BY gid ORDER BY gid, (pt).path), (pt).geom) AS geom
          FROM (SELECT gid, ST_DumpPoints(geom) AS pt FROM MYLINE) as dumps
        )
        SELECT
            gid,
            sum(st_length(geom)) over (order by path asc rows between unbounded preceding and current row)-st_length(geom)  as l_start,
            sum(st_length(geom)) over (order by path asc rows between unbounded preceding and current row)          as l_end,
            ST_StartPoint(geom)                                             as g_start,
            ST_EndPoint(geom)                                               as g_end,
            path,
            geom
        FROM segments WHERE geom IS NOT NULL
)
SELECT
a.gid,
a.path,
a.l_end,
b.g_start,
-- HERE I try to calculate the bisector angle:
atan((ST_Y(a.g_start)-ST_Y(a.g_end))/(ST_X(a.g_start)-ST_X(a.g_end))) + asin((ST_X(a.g_end)-ST_X(a.g_start))/SQRT((ST_X(a.g_end)-ST_X(a.g_start))^2+(ST_Y(a.g_end)-ST_Y(a.g_start))^2)) as phi,
a.geom
FROM SUB a, SUB b
WHERE a.path = b.path-1
AND a.gid = b.gid
)
-- Intermediate node
SELECT gid, path+1 as path, ST_SetSRID(ST_MakePoint(ST_X(g_start)+cos(phi)*l_end,ST_Y(g_start)+sin(phi)*l_end),2056) as geom
FROM POINT_LINE
UNION
-- First node of the line
SELECT gid, path as path, ST_StartPoint(geom) as geom
FROM POINT_LINE
WHERE path = 1

Unfortunately, the node are not always in the right position, so I guess that my angle is sometimes wrong.

Do you have an idea on how to solve this problem?


texts: Find and destroy duplicate OSM multipolygon/polygons in TileMill/PostGIS

texts:
In OpenStreetMap there are sometimes multipolygons where attributes are present on both the polygon and the multipolygon relation, such as the Grampians National Park:

http://www.openstreetmap.org/browse/relation/2703380
http://www.openstreetmap.org/browse/way/201068605

These polygons get loaded by OSM2pgsql separately, so I get duplicate polygons in PostGIS. Both match a query like select * from planet_osm_polygon where leisure='nature_reserve'.

The problem is that when rendering such polygons using a transparent style (eg, over a terrain basemap) they appear darker than other single polygons.

I can't think of any way to avoid this happening, short of not rendering any multipolygons, or only multipolygons.

Any suggestions? Pre-processing in PostGIS is ok if needed.


texts: Can graduated values be added as an attribute in QGIS?

texts: With QGIS Layer Properties > Graduated, we can get values sorted by some criteria as Natural Breaks (Jenks). Is there a way to add this as attribute values of the data? For example, if a specific field is classified into 4 classes, I would like to add this to a field called class with values of 1, 2, 3, and 4.


texts:
I have a script to replace links for pictures from an old folder to a new folder.

  Truncates the link to '\env.gov.ab.ca\Parks\CityWorks\Referenced\Franchere_Bay_PRA', what I want is '\\env.gov.ab.ca\Parks\CityWorks\Referenced\Franchere_Bay_PRA'

 import arcpy
    arcpy.env.workspace = "C:/Users/vincent.law/Documents/CW/CW_Franchere_Bay_PRA.gdb"

    a2 = """def rep_field(in_fld, rep_value):
        targets = ['W:\OPERATIONS\CPI\Asset_Management\GNSS Projects\Parks Inventory\Franchere Bay PRA\FRANCHEREBAYPRA-20150804-CD~files','W:\OPERATIONS\CPI\Asset_Management\GNSS Projects\Parks Inventory\Franchere Bay PRA\FRANCHEREBAYPRA_AUG4_2015_JL~files','W:\OPERATIONS\CPI\Asset_Management\GNSS Projects\Parks Inventory\Franchere Bay PRA\FRANCHERE BAY PRA-04-08-15-TCN~files','W:\OPERATIONS\CPI\Asset_Management\GNSS Projects\Parks Inventory\Franchere Bay PRA\FRANCHERE BAY PRA_2015-08-04_CS~files']
        for targ in targets:
            in_fld = in_fld.replace(targ, rep_value)
        return in_fld"""

    a1 = """rep_field (!Picture!,'\\\env.gov.ab.ca\Parks\CityWorks\Referenced\Franchere_Bay_PRA')"""

    arcpy.CalculateField_management(in_table="Bench_and_Table", field="Picture", expression = a1, expression_type="PYTHON_9.3", code_block = a2)


texts: When to use spatial regression for the proximity analysis

texts:
I am not sure about whether I should use spatial regression for the proximity analysis I am trying to do.
I’ll just describe the scenario and goal:
People claim that convenience stores lower nearby home property values. I want to actually prove this. So I want to essentially prove: is there a statistically significant relationship between the value of your house and how close you are to a convenience store?
Here are the data layers I am working with for the city of focus:

vector polygons of houses in 2010 (contains property values of each house)
vector polygons of houses in 2015 (contains property values of each house)
vector polygons of convenience stores in 2015

Basically I am thinking I would need to first find out how much property values have changed between 2010 and 2015, by taking the 2015 houses layer and subtracting the 2010 houses layer. From this, I would produce a layer that would show the houses in either red (negative) or green (positive) with different shades. So in this new layer, if a house were red, it would mean the house went down in value and if it were green, it would mean the house went up in value. I would imagine these being in different shades, so pink would mean just went down a little bit in value, and dark red would mean went way down in value, and similarly for green, where light green means went up a little bit in value and dark green means went way up in value. Of course, I would first need to see if there is spatial autocorrelation associated with property values.
I would now have all of these red and green houses of various shades and convenience stores. My question is, how do I show if there is a relationship between proximity to convenience stores and the property values of the houses? Normally I would think this would just require regression, but since this is dealing with proximity, wouldn’t this require spatial regression or a geographically weighted regression? Does my approach sound reasonable? I am just trying to figure out what the output of this regression would look like and how to interpret it. Which tools and methods in QGIS would I need for this?


texts: What does coordinate system WG 21° mean?

texts:
I have a land survey document for a piece of land, it contains absolute coordinates that I would like to locate on google maps, but I don't know how to correctly interpret the data.

What does the "System: WG 21°" mean and how do I convert these coordinates to decimal notation for google maps?

The land is located in South Africa.

My Attempt:

I considered that the constants were somehow the baseline values, and tried adding the constants to each coordinate value, but the result doesn't match decimal notation. In any case, in South Africa, the latitude is negative and the longitude is positive. Both latitude and longitude are positive in this case, which makes me think that they're not latitudes or longitudes at all, but rather the result of a more complicated transform. The transform is probably indicated by the WG21 (hence the question), which is where I am currently stuck.


texts: Rebuilding a Polygon from the points within it, in PostGIS

texts:
This question is an update on the problem formulated here.

I wanted to redefine the border of a set of polygons, in order to "smooth it". My first approach (described in the former question) was to grab the points from the border and throw them at the Concave Hull function. However, the "pixelated" look of the border was not improved, since I am keeping the same set of points to define the border.

Thus my second approach, was to use the "raw" points within the polygon, that were used to generate it in the first place.

The idea would be to to something like this:

do a point in polygon query (ST_Contains) to find out the points
inside each polygon.
send each collection of points inside a polygon to the concave hull
constructor; this implies a group by query regarding the polygon
geometry.

Considering the raw points table (tweets_malaga) and polygon table (exploded), I figured the query would be something like this:

SELECT ST_ConcaveHull(ST_UNION(tweets_malaga.geom),0.8) from exploded,tweets_malaga WHERE ST_Contains(exploded.geom,tweets_malaga.geom) group by exploded.geom;

However I am getting a weird error, stating:

ERROR:  Operation on mixed SRID geometries
CONTEXT:  PL/pgSQL function st_concavehull(geometry,double precision,boolean) line 71 at assignment

This bugs me, as all the geometries are correctly set to 4326 :-(

UPDATE:

I was able to make the above query work, setting the target_percent parameter of ST_ConcaveHull to 0.99, which is very close from having a convex hull. This parameter controls the target percent of area of convex hull the PostGIS solution will try to approach before giving up or exiting.
If I set this parameter to anything less than 0.99 it will fail with the mixed SRID geometries error.

I was wondering if this had anything to do with having holes in the geometry, so I fiddled around with the allow_holes parameter of the ST_Concave without any success (in any case, the default is false).

As you can see bellow I was able to generate the concave hull, but I won't mark this question as solved as I don't understand why it is not possible to decrease the target_percent.


texts: QGIS startup error (symbol lookup error, undefined symbol) after system update

texts:
I'm using QGIS 2.4.0 on Linux Mint 16 MATE 64-bit. This afternoon I did a normal system update via apt-get update && apt-get upgrade followed by apt-get autoclean && autoremove. But then suddenly my QGIS stopped working. This is the error message I got

/usr/bin/qgis.bin: symbol lookup error: /usr/lib/x86_64-linux-gnu/libQtWebKit.so.4: undefined symbol: _ZN10QSslSocket16staticMetaObjectE

I tried to uninstall (apt-get purge) & reinstall QGIS but still got the same error message

Anybody knows how to fix this?


texts: Custom marker based on geojson attribute in Leaflet

texts:
I would like to customize a marker to add an id (number) on top of it, this id is coming from an attribute of the geojson ("id" in the following example)

{
"type":"FeatureCollection",
"features":[
    {
      "type":"Feature",
      "geometry":{
        "type":"Point",
        "id":1,
        "destination":"Police station",
        "coordinates":[
          "4.3520735404177699","50.846502697562201"
        ]
      },
    {
      "type":"Feature",
      "geometry":{
        "type":"Point",
        "id":2,
        "destination":"Museum",
        "coordinates":[
          "4.347903728485107","50.846658450225334"
        ]
      },
...

I'm stuck with L.ExtraMarkers :

var json = L.ExtraMarkers.icon({
            icon: 'fa-number',
            number: 'x',
            iconColor: 'white',
            markerColor: 'black',
            shape: 'circle',
            prefix: 'fa'
            });

    var geojsonLayer = new L.GeoJSON.AJAX("query_json.php?&do=DCT1",{pointToLayer: function (feature, latlng) {
                return L.marker(latlng, {icon: json});
    },

    onEachFeature: function (feature,layer) {
        var popup = [];
        for (var prop in feature.properties) {
            popup.push(prop + ": " + feature.properties[prop]);
        }
        layer.bindPopup(popup.join("<br />"));
    }}).addTo(map);

number: 'x', should be the variable "id" coming from the geojson.

I'm also stuck with L.circleMarker :

var stylegeojson = {color:'white', fillColor:'blue', radius:10, opacity:1.0, fillOpacity:1.0, weight:2};

var geojsonLayer = new L.GeoJSON.AJAX("query_json.php?&do=DCT1", {pointToLayer: function (feature, latlng) {
                return L.circleMarker(latlng, stylegeojson);
                // what should I add here to print "id" value on top of it ?
    },
    }).addTo(map);

I think I have to add an

  onEachFeature: function(feature, id) {

but I don't know to write it.

Ideally I'd also like to have a different marker color per "destination"

I am wondering if someone already did this before or has an idea how to do this?


texts: Writing mathematical expression in Google Earth Engine

texts:
I am trying to run the expression below in Google Earth Engine to calculate Index-base Built-up Index(IBI)
var ibi = l8compositeMasked.expression(
    '(2*SWIR/(SWIR+NIR)-((NIR/(NIR+RED))+(GREEN/(GREEN+SWIR)))) / (2*SWIR/(SWIR+NIR)+((NIR/(NIR+RED))+(GREEN/(GREEN+SWIR))))', {
      'RED':l8compositeMasked.select('B3'),
      'GREEN':l8compositeMasked.select('B4'),
      'NIR':l8compositeMasked.select('B5'),
      'SWIR':l8compositeMasked.select('B6')
    });

Map.addLayer(ibi, {min:-1, max:1}, 'ibi', false);

but I keep getting the error below.
Layer shown property is not a boolean: ee.Image({
  "type": "Invocation",
  "arguments": {
    "DEFAULT_EXPRESSION_IMAGE": {

Below is the link to the script
https://code.earthengine.google.com/bdb6397411fe912861900cba04bb9172
How can I resolve this?


texts: Save fails on raster object created from NumpyArrayToRaster?

texts:
I am trying to save an ESRI grid built from a numpy array and am getting an error I've never seen before:

  Unhandled Exception: System.IO.FileNotFoundException: Could not load file or assembly
  'GpMetadataFunctions, Version=10.2.0.0, Culture=neutral,
  PublicKeyToken=8fc3cc631e44ad86' or one of its dependencies. The
  system cannot find the file specified. File name:
  'GpMetadataFunctions, Version=10.2.0.0, Culture=neutral,
  PublicKeyToken=8fc3cc631e44ad86'    at
  System.Reflection.Assembly._nLoad(AssemblyName fileName, String
  codeBase, Evidence assemblySecurity, Assembly locationHint,
  StackCrawlMark& stackMark, Boolean throwOnFileNotFound, Boolean
  forIntrospection)    at
  System.Reflection.Assembly.InternalLoad(AssemblyName assemblyRef,
  Evidence assemblySecurity, StackCrawlMark& stackMark, Boolean
  forIntrospection)    at System.Reflection.Assembly.InternalLoad(String
  assemblyString, Evidence assemblySecurity, StackCrawlMark& stackMark,
  Boolean forIntrospection)    at System.Reflection.Assembly.Load(String
  assemblyString)    at GetManagedType(Char* _assembly, Char* _path,
  Char* _className)

  WRN: Assembly binding logging is turned OFF. To enable assembly bind
  failure logging, set the registry value
  [HKLM\Software\Microsoft\Fusion!EnableLog] (DWORD) to 1. Note: There
  is some performance penalty associated with assembly bind failure
  logging. To turn this feature off, remove the registry value
  [HKLM\Software\Microsoft\Fusion!EnableLog].

I'm using an alternative workflow for now, but I'd still like to try to solve the problem with saving a raster from a numpy array.

Here is what I was trying to do:

read in an ESRI grid,
create an numpy array,
do some manual reclassification, and
save the processed array to a new ESRI grid

The process proceeds as expected through step 3 and fails upon step 4 AFTER the array has been converted to a raster object (see code below). It is when the save method of the raster is called that the error is encountered.

I got the same error when trying to create a numpy array from geo-tiff rasters (failed on calling arcpy.RasterToNumPyArray); I was able to get around this by copying the original .tifs to ESRI grid format first.

Here is my code:

import arcpy
import numpy as np

arcpy.env.overwriteOutput = True
arcpy.env.workspace = 'SwbLandUseParameters.gdb'

# List tables to use in the reclassification
tables = arcpy.ListTables('SwbLandUseParameters*')

# List land use rasters
lu_rasters = arcpy.ListRasters('nlcd*', 'GRID')

for lu_raster in lu_rasters:

    print 'Processing raster', lu_raster

    # Get raster properties
    r = arcpy.Raster(lu_raster)
    ll = r.extent.lowerLeft
    x_cell_size = r.meanCellWidth
    y_cell_size = r.meanCellHeight
    ncol = r.width
    nrow = r.height
    year = lu_raster[-4:]

    # Convert the raster to a numpy array
    lu = arcpy.RasterToNumPyArray(r, ll, ncol, nrow)

    # Create a suite of reclassified rasters using parameter tables
    for tbl in tables:
        temp = np.copy(lu)
        with arcpy.da.SearchCursor(tbl, '*') as cur:
            for row in cur:
                print row
                old_val = row[1]
                new_val = row[2]
                temp[temp==old_val] = new_val
        new_raster = arcpy.NumPyArrayToRaster(temp, ll, x_cell_size, y_cell_size, 0)
        new_raster_name = tbl + '_' + year
        print '  ', new_raster_name
        new_raster.save(new_raster_name)


texts: ArcGIS 10, Converting from date to number

texts:
I am trying to transform a date field like "00:05:40" into an integer field (i.e. 340 sec).
Does anyone know how I could do that conversion please?


texts: Button in ArcGIS add-in extension shows as "Missing"?

texts: I do not understand why my add-in toolbar and button shows up as "Missing" when I open ArcMap. I am not sure whether this is an issue with ArcMap 10.3 or whether I did anything wrong. Have anyone ever encountered this before?


texts:
I´m new to GIS development. I have read about Openlayers, maps servers (I choose Geoserver), geodatabase (in my case i choose PostGIS) and theory of GIS in general. I´m developing a Java EE application that needs GIS functionality, but i´m having some problems defining the high level architecture for my requirements. Before i explain what alternatives i´m considering, i think it´s best if I give the main requirements that the present application must fulfill:

It has to show a map on the browser, with several layers. Among the features that it has to show, are points that represents cities, or users, etc. In some of this features the user is allowed to click and a pop up should appear with related information, and effects like that.
For some of the layers, the users may see only some information (i.e, some of the features) but not all, and this is based on some custom domain logic that its implemented server side with java EE.
The users can modify some map information. They may create certain new features or modify existing ones, but all of these is based on custom logic implemented server side. They can also move through the map, and other users should be capable of seeing this updates.

So, i have a lot of questions regarding how i should design something like this. My first question is:

Should i use WMS or WFS services? As the users will be modifying things on the map, i first think of using WFS services. But, remember that all the operations that modify the map require specific domain logic in order to validate the update, before hitting in the data store. So perhaps i can access directly into the database from java and let the user retrieve the map from a WMS service. At the same time, i´m not so sure about this, because then i´ll be hitting the data base from two sides (Geoserver and jdbc or jpa or whatever) and i don´t know what kinds of problems could arise.
Lets suppose i pick one of the previous services. How could be client-side logic be implemented? I think that with Ajax i can request the application to see if i can validate the operation on the map, and if the validation is succesful then make a call with Openlayers to the service in Geoserver (i don´t know if this is good, since the client ends up making two remote calls, one to the java EE app and other to Geoserver in order to retrieve the ¨results¨ (map)). Or how can i accomplish this? Thinking about this i´m not even sure about why I should be using a map server at all, since i have to do some geographical validations in the map before serving the geographical information.. perhaps i can just query PostGIS and serve the contents directly to openlayers through the ajax call if the previous validations pass? or is it possible to pass the validation logic to geoserver somehow?

As you can see, these are some fundamental questions and decisions I have to made before doing anything else, and I was hoping for advice and other alternatives as well.


texts: How to batch-process georeferencing of rasters (prelude to animation)?

texts: I have many (10k) raster files, all with identical size and pixel dimensions, which need to be georeferenced. I would like to georeference one of them manually, using either QGIS or GRASS, and then apply the same spatial reference to all the others. This seems like it should be a simple task. Is there some command-line utility that I could script to loop over all the files, or is there something in GRASS that will do this?


texts: How can I snap irregular GPS point data to regular points of a route?

texts:
I have irregular GPS points in several shapefiles with meteorological parameters of a mobile measurement campaign.
In order to compare the different shapefiles I would like to snap the irregular GPS points to predefined waypoints (e.g. regular with 1m distance) of the route.
I was able to prepare the regular points with 1m distance. But how can I assign the values of the irregular GPS points to my artificially created points? The number of GPS points is larger compared to the regular points.
I thought about buffering the regular points and take all GPS points within the buffer to calculate the mean of the measured values. What do you think about this idea?


texts: QDockWidget losing its objectName after getting added to QGIS Interface in custom Plugin

texts:
I'm creating a custom QGIS (3.10) plugin, using Python 3, PyQt5 and PyQGIS. I'm facing an issue when adding a QDockWidget to QGIS Interface using iface.addDockWidget() functionality.
Consider this case:
class MainClass:
    def __init__(...):
        button = QPushButton()
        button.clicked.connect(self.add_dock)

    def add_dock(self):
        dock = QDockWidget('Custom Dock Widget', iface.mainWindow())
        dock.setObjectName('CustomDock')
        DockWidget(dock)    #Focus Line

class DockWidget:
    def __init__(self, dock):
        self.dock_widget = dock
        iface.addDockWidget(Qt.RightDockWidgetArea, self.dock_widget)

In such case, a New Dock Widget with Title: 'Custom Dock Widget' is added to QGIS Right Dock area. It is also added to View > Panels with a Checkable Button. But it loses it's objectName. Because if I use the Python Console in QGIS and run the following code:
#Testing Code
docks = iface.mainWindow().findChildren(QDockWidget, 'CustomDock')
for i, d in enumerate(docks):
    print(i, d.objectName())

I do not get any output. When I change the Focus Line to:
dock_dialog = DockWidget(dock)

Still no output. BUT if I change the Focus Line to:
self.dock_dialog = DockWidget(dock)

Then I get the output as:
> 0 CustomDock

Dock Widget gets successfully added to QGIS Right Dock area in all three cases, but it loses it's objectName in 2 of them. I am not keen of creating an instance variable when there is no need to create any.
Can there be any better workaround this? And can anyone help me understand why this is happening? Something related to free object getting cleared out of memory, or something else entirely?


texts: Open source alternative to FME from Safe Software for Spatial ETL?

texts:
Just wondering if there is any open source solution which comes close to safe FME?
I would like to integrate it into my workflow, but it´s just another few thousand which my employer has to give out for software.

Does anyone know of anything similar?


texts: Plotting cross section using python, gdal, tiff?

texts:
How to plot a cross section between two pairs of lat/long using a tiff file in python?

At first glance, it seems i have to calculate points position between these two desired points and then looking for nearest points in tiff to get their elevation.

Until now i could:

# Read tiff data
ds = gdal.Open('cliped.tif')
# Read it as Array
a=ds.ReadAsArray()

but the problem is data point has no projection so i can not search for lat/long pairs!


texts: NameError: global name 'QgsDataSourceURI' is not defined

texts:
Following on from my previous question I cannot seem to get my function to run using script runner in QGIS. However, if i type then out in the Python Console then they work and load the layer in question.

Is it something to do with how i! am using the run_script function at the bottom of my code?


texts: Python batch script to delete Attributes

texts:
I currently have a script that will select and delete all the values less than 0.1 in a column. The script below works, however it specifies a specific .shp, I want to batch this so it runs through all the shapefiles in a geodatabase

import arcpy
arcpy.SelectLayerByAttribute_management("Sr90_contour","NEW_SELECTION",'"Contour" < 0.1')
arcpy.DeleteFeatures_management("Sr90_contour")


texts: Using subtype property of the esriaddinx for a AddIn button

texts:
I am trying to utilize the subtype property of an AddIn toolbar button in C#. In older versions of ArcObjects this could be accomplished using ICommandSubtype and extending the BaseCommand. However, I cannot see how to do it in the current ArcGIS 10.x AddIn model.
I would like to implement this functionality, because I have multiple buttons that will simply write out a single field value to a selected feature. I do not want to have to create multiple classes for each value. Being able to use the subtype capability of a single class would be ideal. I noticed that the "subtype" is available in the intellisense of the toolbar button in the .esriaddinx, so can this be done?
Is there a way to extend the ESRI.ArcGIS.Desktop.AddIns.Button and use the "subtype" property?
I am developing for ArcGIS 10.4.


texts: Handling 2s in cost path analysis output?

texts:
I am trying to compute shortest distance to stream from various point locations on the landscape. I am using the cost path tool in spatial analyst to do this . I am using a DEM as a cost distance raster and flow direction raster as a cost back link raster. The tool runs fine. What I see is where there are points close to each other multiple paths merge( as it is supposed to according to the algorithm).

"When multiple paths merge and follow the remaining distance back to a source on the same route, the segment where the two paths travel together is assigned the value 2. The merged portion of the path cannot be assigned the value of one of the paths, since the merged portion belongs to both routes."

I am later converting this cost path raster layer to a polyline to compute geographic distances of each path. But the problem is with these lines coded as "2". They are not attached to the parent line. How do I covert cost path to polylines( one line for each point from source to destination). I tried the unsplit lines tool but it still puts the "2s" as separate line segment.

I'm using ArcMap 10.1.


texts: How to Configure sextante in Qgis 1.9 ubuntu

texts: I am using Ubuntu 13.04 and the most up to date version of Qgis 2.1.0-Master. I am trying to get sextante to work correctly and have a hard time. Specifically SAGA 2.1.0 and GRASS 6.4.2 will not execute correctly. When I try to run any of the grass tools qgis hangs then crashes. When I try to run SAGA tools it tells me saga is not installed or needs to be configured. The latest version of Sextante seems to have removed the SAGA folder input field which I am assuming needs to be updated. I am not sure why grass is crashing at all. I have another similar question open HERE. Any help would be greatly appreciated.


texts: Creating QGIS multi-layer Topo Map of Lorraine Region of France?

texts:
I have a project I'm working on where I need a QGIS multi-layer Topo Map go the Lorraine Region of France between Metz and Nancy as a starting point, but I don't know where to start.  What I need is a georefenced map so that I can utilize the roads, raliroads and river layers as well as the elevation contour lines.

I'm a newbie, so this is a lot to learn.


texts: Getting geojson data of relations by overpass in Python?

texts:
I am trying to use overpass api for python to get osm data in the format of geojson. I have successfully run this api to get way in the format of geojson.
The code is like:

api.Get('%s(%d)' % ('way',319141887))

But when I try it for relation data like:

api.Get('%s(%d)' % ('relation',4431641))

, it returns no data. (I'm sure that there is a relation with this id because the way 319141887 is one part of this relation.)

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

and I am wondering how to get them.


texts: How to move a point shapefile to the centroid of a polygon?

texts: I have this point shapefiles on top of a bounded polygon. I would like to move each point to the centroid of each bounded polygon. Is there an automatic way in ArcGis 10.1 of solving this? Like in the image shown below.


texts: planet quick-search api not recognizing authentication?

texts:
I can successfully reproduce the first two API calls from this walkthrough: https://www.planet.com/docs/api-quickstart-examples/step-1-search/ (as I have an API key). However, the third call (the POST to 'https://api.planet.com/data/v1/quick-search') gives me:

{
  "message": "Please enter your API key, or email and password.",
  "errors": []
}

I have copied and pasted everything from this example.

Is this expected because I am a free user?

EDIT - Responding to Richard Law - Yes, I have a mac and echo $PLANET_API_KEY shows me my key. As would be expected given the fact that the first two calls worked..

EDIT #2 - Responding to Bosth - I can show you my code, but it really is not that interesting because it's the exact same code as the example. Here it is. I'm running the two requests in one script to make sure it's the same environment..

import os
import requests
from requests.auth import HTTPBasicAuth

geo_json_geometry = {
  "type": "Polygon",
  "coordinates": [
    [
      [
        -122.52227783203125,
        40.660847697284815
      ],
      [
        -122.52227783203125,
        40.987154933797335
      ],
      [
        -122.01690673828124,
        40.987154933797335
      ],
      [
        -122.01690673828124,
        40.660847697284815
      ],
      [
        -122.52227783203125,
        40.660847697284815
      ]
    ]
  ]
}

# filter for items the overlap with our chosen geometry
geometry_filter = {
  "type": "GeometryFilter",
  "field_name": "geometry",
  "config": geo_json_geometry
}

# filter images acquired in a certain date range
date_range_filter = {
  "type": "DateRangeFilter",
  "field_name": "acquired",
  "config": {
    "gte": "2016-07-01T00:00:00.000Z",
    "lte": "2016-08-01T00:00:00.000Z"
  }
}

# filter any images which are more than 50% clouds
cloud_cover_filter = {
  "type": "RangeFilter",
  "field_name": "cloud_cover",
  "config": {
    "lte": 0.5
  }
}

# create a filter that combines our geo and date filters
# could also use an "OrFilter"
redding_reservoir = {
  "type": "AndFilter",
  "config": [geometry_filter, date_range_filter, cloud_cover_filter]
}

stats_endpoint_request = {
  "interval": "day",
  "item_types": ["REOrthoTile"],
  "filter": redding_reservoir
}

# fire off the POST request
result = \
  requests.post(
    'https://api.planet.com/data/v1/stats',
    auth=HTTPBasicAuth(os.environ['PLANET_API_KEY'], ''),
    json=stats_endpoint_request)

print result.text
print

# Search API request object
search_endpoint_request = {
  "item_types": ["REOrthoTile"],
  "filter": redding_reservoir
}

result = \
  requests.post(
    'https://api.planet.com/data/v1/quick-search',
    auth=HTTPBasicAuth(os.environ['PLANET_API_KEY'], ''),
    json=search_endpoint_request)

print result.text

Result from running code:


texts: One-to-many spatial join with results in one row

texts:
I have two feature Datasets. Both are Polygons. Dataset A has very big Areas and Dataset B just small ones. Now, B still overlaps in 2 or 3 Zones from Dataset A.

I want to perform a spatial join and put the result comma seperated in just one row like this:

Zone_ID_A  |  Zone_ID_B
   1       |   2,3
   2       |   2,4
   3       |   5
   4       |   1,6

I know this is possible with ArcMap but i could not figure out a way to do this in QGIS, any ideas?


texts: Add a scale or category to a rule which is only applied to points within original subset in QGIS

texts:
I have subset a shapefile of points to only include those having a value of (for example) 33 in the variable "Bout" using a rule based assignation:

Now I would like to apply a sub-category to that rule which colours the points by date like so:

However this uses the range of the original data rather than the range within the subset. Is there a way to only spread the colour values across the range of the subset data without manually creating each category? It is doable by simply removing the values from the sub-category manually and recalculating the value-based colours but I was hoping for a more elegant coded solution as I am a python amateur.


texts: Creating road center line from Grid-Polygon Shapefile using QGIS

texts:
I have the shapefiles for the polygon part which are colored as yellow but there is no any specific shapefile for the white parts( roads).
Note: I have seen the same question here as below but it didn't work for me using Qgis 3.16.0 with grass 7.8.4.

Finding centrelines from polygons in QGIS?

and

Creating centerline of river in QGIS?


texts: How to install PostGIS correctly with OpenGeo Suite 4.0 community edition on Windows 7?

texts:
I have installed the OpenGeo Suite community edition, version 3, onto Windows 7 multiple times in the past, with virtually no problems.
With the latest version of OpenGeo Suite (4.0.1) (which includes PostGIS 2.1), I am able to install almost all of the components successfully, except PostGIS. With earlier editions, PostGIS was installed as part of the full Opengeo installation process, but now, although PostGreSQL 9.3 is successfully installed, there are no PostGIS or opengeo database templates installed, nor are the spatial functions installed. Nevertheless PostGIS 2.1 is installed.

It seems that OpenGeo Suite 4 intalls all the component but no longer creates a default database or user account in Postgresql/PostGIS.

What steps are necessary to configure OpenGeosuite (windows installation) in order to be able to fix this issue?
thanks


texts:
I'm trying to convert IMG data from USGS's 1/9 arc second DEM to something that Pathloss 5 can read.

Is there a way to do so?


texts: Extraction of stream gradient for branched river networks

texts:
At the moment I am working on the extraction of stream gradients for a large basin, as a consequence, rivers are strongly branched.So far, I read through the relevant posts on this site, but so far I am not able to apply the proposed methodologies to my problem.

I started with the common procedure of extracting the river network from the 30-by-30 m AsterDEM.  What I want to do now, is to measure the gradient of the resulting streams in fixed intervals. The interval for a river reach depends on its mean gradient, which I determined using HEC GeoHMS (the steeper the mean slope, the more often I want to sample). Thus I divided the streams in four categories.

Now, I need to divide the resulting streams in equal intervals. I wanted to make use the technique as proposed here: 1. Thus converting the stream network into a route (using the “linear referencing toolbox”) and then use route events to divide the route in reaches of equal lengths. Unfortunately I have no idea how to prepare an event table in order to divide the routes in equal intervals (e.g. one event every 100 m), especially in a branched network.

I would highly appreciate if you could guide me through the workflow, as I have no experience with the linear referencing tools!

Rafael


texts:
Our map (not publicly visible, sorry) shows owner names and property boundaries.  Sometimes the names are truncated about a vertical line, as in these examples.  The first example also shows one name spans the vertical line.

I'm wondering if the text is only laid over a tile if a small circle around the centre point happens to fall on the tile.  In the case of Kauri Glen Limited in the first example, perhaps the circle falls onto two tiles.  Pure guesswork.

Can anyone explain definitively what's going on, and provide a way to prevent the truncation?


texts: Simple raster calculation in batch mode (for folders)

texts:
I have thousands of raster organized in several folders. All of their cell values have to be divided by 10. Seems simple, but I didn't find any solution to do it in batch. QGIS offers batch mode, but all raster calculators (SAGA, GRASS, GDAL) are way to complex for simply filling the batch parameters.

I'm sure, there must be a simple command to do that calculation:  [Raster Folder] / 10.

(working on ubuntu)


texts: Control custom panes for Leaflet Geojson SVG icons

texts:
In Leaflet v1.03 I have n layers of SVG icons which should be rendered on custom panes. I try to use the panes to control the zindex of items.

Within the L.geoJson options I set the pane as well as the SVG renderer with a pane option itself. Neither works. The SVG render order seems to be random.

Is the renderer option actually supported in Leaflet v1.03? The option itself is not listed in the current documentation.

Would I need to create a layer group in order to control the panes of SVG icons created with the pointToLayer method?

map.createPane('400');
map.getPane('400').style.zIndex = 400;
map.getPane('400').style.pointerEvents = 'none';

new L.geoJson(dots400, {
    pointToLayer: function (feature, latlng) {
        return L.marker(latlng, {
            icon: L.icon({
                iconUrl: 'blue_dot.svg',
                iconSize: [10, 10]
            }),
            interactive: false
        });
    },
    renderer: L.svg({pane: '400'}),
    pane: {pane: '400'},
    onEachFeature: null
}).addTo(map);

map.createPane('500');
map.getPane('500').style.zIndex = 500;
map.getPane('500').style.pointerEvents = 'none';

new L.geoJson(dots500, {
    pointToLayer: function (feature, latlng) {
        return L.marker(latlng, {
            icon: L.icon({
                iconUrl: 'red_dot.svg',
                iconSize: [10, 10]
            }),
            interactive: false
        });
    },
    renderer: L.svg({pane: '500'}),
    pane: {pane: '500'},
    onEachFeature: null
}).addTo(map);

Here is a jsfiddle. https://jsfiddle.net/goldrydigital/abrurqwn/

I want the red dot to be drawn on top.


texts:
I've followed the tutorial available here. The installation worked alright, but opening the tile in a browser leads to a 404 error. The only error from the server is "Failed to read cmd on fd 10".

After a long while, the tile is rendered (here, 1600 seconds):

renderd[18586]: DEBUG: DONE TILE default 13 4264-4271 2984-2991 in 1594.968 seconds
debug: Creating and writing a metatile to /var/lib/mod_tile/default/13/0/16/11/170/136.meta

How can I solve this problem?


texts: Why does data imported from AutoCAD appear different in ArcMap 10.1?

texts:
I am importing information from AutoCAD into Arcmap 10.1.  I found after some issues that I need to give the incoming information geographic location, which works fine.  The problem is that all the data appears the same.  I have two maps with different points of collections, but when I bring over both maps, they give me the same data.

Pic1: is the data that I am pulling over from AutoCad

Pic2: is what appears in Arcmap 10.1.

If you look at the bottom, you will notice that they are both suppose to be the same, but they are not.  The title is different, and also when I bring it over to ArcMap 10.1, it has alot more points of data in the different areas.
I have bringing over alot of different maps, is it possible that my maps are reading data from a different map that I brought over before I was looking at these?


texts:
I need help splitting a linestring into smaller linestrings, using PostGIS.

I need to route to closest point on the linestring, so I need to have more points on the line

Sample Linestring = LINESTRING(11.959813031048 55.3269904906477,11.9613429171704 55.3257458375313,11.9656473425314 55.325719907258,11.9675402524794 55.3251753715195)

I have no luck with the following

select st_astext(st_segmentize(st_linefromtext('LINESTRING(11.959813031048 55.3269904906477,11.9613429171704 55.3257458375313,11.9656473425314 55.325719907258,11.9675402524794 55.3251753715195)'
),2) ) as geom

But something is missing? I would like to have my line split into lines max 10 meters between points.

I hope somebody could get me back on track.

I then dump the line into points

select st_astext((st_dumpPoints(st_segmentize(st_linefromtext('LINESTRING(11.959813031048 55.3269904906477,11.9613429171704 55.3257458375313,11.9656473425314 55.325719907258,11.9675402524794 55.3251753715195)'
),10 ))).geom) as geom


texts: Executing script in Python window of ArcMap?

texts:
I've been writing Python in my text editor so I have a Python script, script.py. I'm still developing it so don't want to put it in my toolbox yet, and at the moment this means copying and pasting it from my text editor to the Python window in ArcGIS Desktop.

Is there a command that lets me call it from the ArcGIS Python window?

Something along the lines of:

run script.py

It seems logical that it should exist, but perhaps not.


texts: In ARCGIS 10.5 mxd I need to put definition query for one layer as Current Year in Year fields , What i need to set the query

texts:
I have a feature class and related info in another table, the table having fields as Year, Year having values like - 2016,2017,2018,2019... Presently I Put the definition query for table Year = 2018 after that I have joined table to feature and published map service and the results perfectly fine. for next year also I need to change year=2019. So that without changing the code every time, I required based on the system date need to change the year. Please provide the code for that how to get the date in select by attribute

Year field is Integer

The Table is SQL Table which is not registered as versioned


texts: How do I construct (or select) an appropriate local Cartesian coordinate system, for an arbitrary small area?

texts:
I'm relatively inexperienced with coordinate reference systems, so please bear with me if I'm not putting this clearly.

I have existing software which works with an abstract (un-referenced) 3D coordinate system.  That is, everything has X,Y,Z coordinates in meters, but those coordinates are not defined relative to any absolute location.

I want users to be able to "attach" this coordinate system to the real world by saying "the (local) origin corresponds to point FOO in the real world."  Here FOO is specified in any well-known coordinate reference system (probably EPSG:3785, but I'll allow anything Geotools can recognize), and the orientation is fixed.  I then want to have well-defined translations between a relative coordinate (X,Y,Z) and the corresponding absolute location (e.g. Y meters North of FOO, X meters East of FOO, Z meters above FOO).

What I think I want, in order to make this happen, is to (a) pick or create an appropriate Cartesian coordinate system for the area around FOO, and then (b) do the simple affine transformation between my internal coordinate system and the geo-referenced one.

All the advice I've seen so far says "use a local projection for small areas, or something like WGS for large areas," but I don't know what the right local projection is in advance, because FOO could be anywhere.  Can I either reasonably (1) look-up a standard, accepted projection for that location, or (2) define a good new projection for that point?

I'm working in Geotools, but I'm perfectly happy with a general conceptual answer.  Thanks!


texts: How to set time field when generating a GPX track with Python OGR?

texts:
I want to generate a GPX file with Python OGR. In particular I want to generate a tracks (layer of features of OGR type wkbMultiLineString)

My script is below. I get the error ERROR 1: No such field: 'time'
Also it won't let me create a field time because it is a default field in GPX.

Any ideas how to solve that?

import ogr, osr

in_data_source = ogr.Open('test.geojson')
in_layer = in_data_source.GetLayer()
feature_count = in_layer.GetFeatureCount()

driver = ogr.GetDriverByName("GPX")
out_data_source = driver.CreateDataSource("t1.gpx")
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
out_layer = out_data_source.CreateLayer("t1", srs, ogr.wkbMultiLineString)

line = ogr.Geometry(ogr.wkbLineString)

c = 0
while c in range(feature_count-1):
    print c
    f1 = in_layer.GetFeature(c)
    f2 = in_layer.GetFeature(c+1)
    g1 = f1.GetGeometryRef()
    g2 = f2.GetGeometryRef()

    line.AddPoint(g1.GetX(), g1.GetY())
    line.AddPoint(g2.GetX(), g2.GetY())
    time = f1.GetField("datetime")

    feature = ogr.Feature(out_layer.GetLayerDefn())
    # feature.SetField("time", time)
    feature.SetGeometry(line)
    out_layer.CreateFeature(feature)

    c += 1

This problem can also be solved by using the gpxpy library, though I am interested in solving it with OGR.


texts: QGIS Composer value based on http request param

texts:
is it possible to set value in QGIS Print Composer based on HTTP request param?

For example i want to add table in composer with something like "value = $requestedvalue", and $requested value should be in http request: serveraddress?map=...&SERVICE=WMS&REQUEST=GetPrint&...&$requestedvalue=something


texts: How to find the location of parking areas in OSM with QGIS?

texts:
I am making a small size map with QGIS.
I need to mark the parking areas in three streets in Philadelphia: 54th street, Residential street and Berks street.

I was wondering, is there a way to get the locations of these parking areas?

Maybe it can be download in a similar way to the one that is used to download street and buildings layer from
OSM?

I could not find this data in the OSM plugins that I am using in QGIS: OSMDownloader and osmSearch.

If it is of any help, here is the map area and the streets are marked with arrows.


texts: L.Control.Search add layers to map when added to map

texts:
I implemented a search functionality in my leaflet map using L.Control.Search. As the doc and the examples say, you have to pass for the Control the layerGroup you want to perform the search on. My layers are alreday clustered and added to the map but when I add the L.Control.Search the layers are added too and a mess is displayed on my map.
Here is my code of how I'm doing things :
var searchLayer = L.layerGroup([sibe, parc, daya, faune, biotope])
        //... adding data in searchLayer ...

    var searchControl = new L.Control.Search({
        layer: searchLayer,
        propertyName: 'titre',
        circleLocation: false,
        moveToLocation: function(latlng, title, map) {

            map.setView(latlng, 14); // access the zoom
        }
    })

    mymap.addControl(searchControl);

   var sibesMarkers = L.markerClusterGroup();
   var parcsMarkers = L.markerClusterGroup();
   var dayaMarkers = L.markerClusterGroup();
   var fauneMarkers = L.markerClusterGroup();
   var biotopMarkers = L.markerClusterGroup();

   sibesMarkers.addLayer(sibe);
   parcsMarkers.addLayer(parc);
   dayaMarkers.addLayer(daya);
   fauneMarkers.addLayer(faune);
   biotopMarkers.addLayer(biotope);

Is there any way to prevent the L.Control.Search from adding the layers to the map.


texts: I have geoTIFF raster with cells size 5x5 m in WGS84. I would like to make a point layer with values from raster. I use ArcMap extension: raster to point, but it makes a point net with size 5x~9 m, not each cells centroids. Do you have some idea why?


texts: Adding two fields in attribute table and assigning result to third field if some cells contain null values using ArcGIS Pro Field Calculator

texts:
I am using ArcGIS Pro, field calculator to add the values of two fields (Result_A and Result_B)  and assign the value to new Field Combined_Result. The type for all three fields is "double". Both fields for Result_A and Result_B have some cells with null values. This is generating errors in their addition. In the Calculate Field window (for field Combined_Result) For the calculate field, I am using the following expresssion:
=!Result_A! + !Result_B!

I get the following error:

WARNING 002858: Certain rows set to NULL due to error while evaluating
python expression: TypeError: unsupported operand type(s) for +:
'float' and 'NoneType'

I have tried changing the fields to float but get the same error.
Should I convert null values to zero before running the expression?
I have beginner's level python skills and can set up a code to convert null values if you point me in the right direction.


texts: How to prevent ArcGIS SOC from locking Oracle tables?

texts:
I have a python geoprocessing that performs CRUD operations on feature classes and oracle tables.

Sometimes, when handling Oracle tables, ArcGIS SOC gets hung and locks the Oracle table. I need then to restart ArcGIS SOM service to unlock the table.

Is there a way to prevent ArcGIS SOC from locking tables and getting hung?


texts: Unable to Generate Popup with Leaflet From Geoserver

texts:
I've been learning to use GeoServer. I found a tutorial: https://www.e-education.psu.edu/geog585/node/765. I followed all of the steps, however the Leaflet map does not generate a popup even when I go to localhost:8080.

Here is the code where the popup is generated:

          function Identify(e)
          {
            // set parameters needed for GetFeatureInfo WMS request
            var sw = map.options.crs.project(map.getBounds().getSouthWest());
            var ne = map.options.crs.project(map.getBounds().getNorthEast());
            var BBOX = sw.x + "," + sw.y + "," + ne.x + "," + ne.y;
            var WIDTH = map.getSize().x;
            var HEIGHT = map.getSize().y;

            var X = Math.trunc(map.layerPointToContainerPoint(e.layerPoint).x);
            var Y = Math.trunc(map.layerPointToContainerPoint(e.layerPoint).y);

            // compose the URL for the request
            var URL = 'http://localhost:8080/geoserver/geog585/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&LAYERS=philadelphia:FarmersMarkets&QUERY_LAYERS=philadelphia:FarmersMarkets&BBOX='+BBOX+'&FEATURE_COUNT=1&HEIGHT='+HEIGHT+'&WIDTH='+WIDTH+'&INFO_FORMAT=application%2Fjson&TILED=false&CRS=EPSG%3A3857&I='+X+'&J='+Y;

            //send GetFeatureInfo as asynchronous HTTP request using jQuery $.ajax

            $.ajax({
               url: URL,
               dataType: "json",
               type: "GET",
               success: function(data)
               {
                 if(data.features.length !== 0) {  // at least one feature returned in response
                   var returnedFeature = data.features[0]; // first feature from response

                   // Set up popup for clicked feature and open it
                   var popup = new L.Popup({
                     maxWidth: 300
                   });

                   popup.setContent("" + returnedFeature.properties.NAME + "
" + returnedFeature.properties.ADDRESS);
                   popup.setLatLng(e.latlng);
                   map.openPopup(popup);
                }
              }
             });
           }

           map.addEventListener('click', Identify);

         }

Does anyone know why this may not be generating anything? I also tried it with the OpenLayers code and it worked perfectly.

I'm very new to this.


texts: GeoServer 2.7.4 admin/geoserver password remains the same

texts:
This solution presented here: How do I change geoserver admin password? doesn't work.

I have GeoServer 2.7.4 up and running, I went into the admin pane/ users, changed the default password, still no joy. Restarted the server, the old password still works. Interesting enough, the new password works too.

I went restarting the whole machine, with no luck. Then I deleted the admin user after creating a replacement one (admin2 with a different password), so now the "admin" account is gone, still I can use it to login. I know this is not the master root account, but why I can't get rid of it? I think all the other GeoServer instances of the same version have the same vulnerability.

Any ideas?


texts: Toggle visibility of MVT features in OpenLayers 6

texts:
Following code add layer with MVT features:

class Mvt {
    constructor() {
        this.format = new MVT({
            idProperty: 'id',
            geometryName: 'geom',
            featureClass: Feature,
        }),
        this.source = new VectorTileSource({
            url: 'http://{a-d}.myhost/mvt/{z}/{x}/{y}',
            format: this.format,
            tileGrid: new TileGrid({...}),
        })
        this.layer = new VectorTileLayer({
            declutter: true,
            source: this.source,
            style: feature => someStyleFunction(feature)
        })
    }
}

const mvt = new Mvt()
map.addLayer(mvt.layer)

Each feature has layer attribute with correspond layer name,
so it's possible to set layer list for reading of features by specifying that list in format property (or using setLayers method) as mentioned here

Knowing about this, it's tempting to implement toggling visibility of "layers" where "layer" is collection of features with the same layer attribute.

I've tried mvt.format.setLayers([...]), but changes take effect after subsequent request only(i.e zoom/pan interaction).

It would be perfect to redraw the tiles immediately or to refresh source of layer at least.

Does anyone know a trick to do it?


texts:
Running the command:

gdaldem hillshade in.tif out.tif -z 30

I'm getting artifacts like:

and at -z 1:

How might I correct this? Here is the gdalinfo for the input GeoTiff. All dimensions should be in meters:

Driver: GTiff/GeoTIFF
Files: in.tif
Size is 58273, 58273
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (-20037508.339999999850988,20037508.339999999850988)
Pixel Size = (687.711576201671392,-687.711576201671392)
Metadata:
  AREA_OR_POINT=Area
  lat#actual_range={-90,90}
  lat#axis=Y
  lat#long_name=latitude
  lat#standard_name=latitude
  lat#units=degrees_north
  lon#actual_range={-180,180}
  lon#axis=X
  lon#long_name=longitude
  lon#standard_name=longitude
  lon#units=degrees_east
  NC_GLOBAL#Conventions=CF-1.7
  NC_GLOBAL#GMT_version=6.0.0 [64-bit]
  NC_GLOBAL#history=grdmath SRTM15_land_masked.nc debug/polished.unmasked.grd AND = SRTM15.nc
  NC_GLOBAL#node_offset=1
  NC_GLOBAL#title=Produced by grdmath
  z#actual_range={-10997.708984375,8627}
  z#long_name=z
  z#_FillValue=nan
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-20037508.340,20037508.340) (180d 0' 0.00"W, 85d 3' 4.06"N)
Lower Left  (-20037508.340,-20037508.340) (180d 0' 0.00"W, 85d 3' 4.06"S)
Upper Right (20037508.340,20037508.340) (180d 0' 0.00"E, 85d 3' 4.06"N)
Lower Right (20037508.340,-20037508.340) (180d 0' 0.00"E, 85d 3' 4.06"S)
Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
Band 1 Block=58273x1 Type=Float32, ColorInterp=Gray
  NoData Value=nan
  Metadata:
    actual_range={-10997.708984375,8627}
    long_name=z
    NETCDF_VARNAME=z
    _FillValue=nan


texts: Can not export some annotaton layers using QuickExport in ArcObjects.NET

texts:
I develops add-in for ArcMap that exports different layers to .mif format using ESRI.ArcGIS.DataInteroperabilityTools.QuickExport

Geoprocessor gp = new Geoprocessor();
gp.OverwriteOutput = true;
QuickExport qe = new QuickExport();
qe.Input = layerPath;
qe.Output = string.Format("MIF, {0}", exportPath);
gp.Execute(qe, null); // error is here

Common geometry layers (points, lines, polygons) are exports without any problems.

But I get Error HRESULT E_FAIL has been returned from a call to a COM component when tries export some annotation layers, and some annotation layers exports successfully. Inner exception is null.
Haw can I debug this error?
What difference could be between annotations?
And what could be reason of such behaviour?

Update:

Also I tried such way:

IVariantArray parameters = new VarArrayClass();
parameters.Add(layerPath);
parameters.Add(string.Format("MIF, {0}", exportPath));
RunTool("QuickExport_interop", parameters, null, false);

public static bool RunTool(string toolName, IVariantArray parameters, ITrackCancel TC, bool showResultDialog)
{
    Geoprocessor gp = new Geoprocessor();
    IGeoProcessorResult result = null;
    // Execute the tool
    try
    {
        result = (IGeoProcessorResult)gp.Execute(toolName, parameters, TC); // error is here
        string re = result.GetOutput(0).GetAsText();
        if (showResultDialog)
            ReturnMessages(gp);
        if (result.MaxSeverity == 2)
            return false;
        else
            return true;
    }
    catch (COMException err)
    {
        string message = err.Message;
        Exception e = err.InnerException;
        while (e != null)
        {
            message += "\n" + e.Message;
            e = e.InnerException;
        }
        MessageBox.Show(message + " in RunTool");
        ReturnMessages(gp);
    }
    catch (Exception err)
    {
        MessageBox.Show(err.Message + " in RunTool");
        ReturnMessages(gp);
    }
    return false;
}

private static void ReturnMessages(Geoprocessor gp)
{
    if (gp.MessageCount > 0)
    {
        for (int Count = 0; Count <= gp.MessageCount - 1; Count++)
        {
            Console.WriteLine(gp.GetMessage(Count));
        }
    }

}

But result was same as in first way.

QuickExport tool works only in ArcPy for me both for annotation layers and common geometry layers:

arcpy.QuickExport_interop(layer_path, 'MIF, {}'.format(export_path))


texts: Is it possible to calculate NDVI using a landsat 8 .img file?

texts:
I have a mosaic image with this file extension (.img) When I use image analysis in arcmap.10.4 to calculate NDVI I get values the following values

High: 3
Low : -5

NDVI values are supposed to be between -1.0 and 1.0 so this is wrong. Just wondering if it's even possible to calculate NDVI with an Erdas Imagine Image in arcmap.


texts:
I am trying to edit bus routes in QGIS.
Multiple bus routes can run on the same street.

When I want to change one bus route, I select one of them, then I enable editing mode. With the Node tool, I then edit or delete a node.

The problem is that this change is not applied to the bus route I had selected originally.

How can I make sure that I edit a specific bus route if multiple routes use the same street?


texts: Modify Attribute Table while opening using PyQGIS

texts:
I have a problem related to this question and the script given as the answer. I can add some functionalities to Attribute Table (AT) by the script in the post. But according to the post, I have to open the AT using showAttributeTable method, because the method returns the reference of the AT created/opened.

Of course, I can get references of all opened ATs using the following line.

tables = [w for w in qApp.allWidgets() if w.objectName() == 'AttributeTable']

Let's say, I want to add a button to all ATs to be opened or sort by a column for a specific layer. Somehow, I think I need to catch "the AT's opening event" or "a widget/child was added to iface.mainWindow()" etc.

I have looked at Qt5 and QGIS API(QgsApplication) documentation, but I couldn't find anything helpful or I missed something.

A pseudo code for a possible solution:

def something_opened(something):
    if something is an attribute_table and active_layer is foo_bar:
        do something

main_window.addedsomething.connect(something_opened)

Note: @Ben's answer is pretty nice. But, it still requires me to focus on the AT.


texts: How to compute surface type of raster or any point using ArcGIS?

texts:
I have a raster image whose values are a digital elevation model of terrain.

I need to know the surface type (topographic characterization) of particular points along transects of this DEM.

How to find it?

(In my application I have two points: a transmitter end and receiver end.
I find the highest point between the transmitter end and receiver end: see the figure below.  I would like to know whether the highest point is smooth or is a knife edge; that is, to within the resolution of the image it appears to be non-differentiable.)


texts: Why an extra iteration on the first feature of a GeoJSON layer?

texts:
I'm using Python and GDAL to iterate over the features in a GeoJSON layer and I'm seeing this strange behavior where the first iteration is repeated a second time before moving onto the next feature.

I only observe this when:
- opening the vector with gdal.GA_Update option
- including a layer.SetFeature() within the loop
- gdal version 2.4.1 (not observed in 2.2.4)

This example shows there is only 1 feature, but the print statement within the loop is called twice and prints FID 0 both times. With a GeoJSON containing more features, I observed that only the first feature iteration is repeated twice.

from osgeo import gdal

vector_path = 'test.geojson'

with open(vector_path, 'w') as file:
    file.write('{"type": "FeatureCollection", "name": "test", "features": [{ "type": "Feature", "properties": { "foo": 1 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -1, -1 ], [ -1, 1 ], [ 1, 1 ], [ 1, -1 ], [ -1, -1 ] ] ] } }]}')

vector = gdal.OpenEx(vector_path, gdal.OF_VECTOR | gdal.GA_Update)
layer = vector.GetLayer()

print(layer.GetFeatureCount())

for feature in layer:
    print(feature.GetFID())
    layer.SetFeature(feature)

layer = None
vector = None

Is there something wrong with my GeoJSON file? Or with the GeoJSON driver?


texts: How to get a visualized overview over huge PostGIS table?

texts:
within a python application a uge PostGIS geometry table (Polygon) is polled for certain processing. Yesterday we have set our testing system back to an earlier snapshot and failed to notice that this application still used the table on our testing system. So we have lost some million records of much more million records.
No we want to visualize the remaining 26 million records to evaluate the 'damage'. We tried PostGIS Layer in QGis (no chance), export to shape or GML (no chance, 28GB GML file), export to GTiff (aborted), we tried to union to reduce number of geometries (query is still running after 2 days) ...
Is there a method out there to 'quick' visualize such an amount of data? Our testing system is running on a server so we can add additional processors and RAM if needed.

Any help is welcome.


texts: ArcGIS, detect adjacent polygons with different values

texts:
I am developing an ArcMap extension and I need to validate a particular case. I need to detect adjacents polygons that have one value equal, and the others not equal. I you take a look at the image, the yellow polygons are invalid and I need to inform the user. The sector is "Sector 2" on both polygons but, the next value is not the same (100 vs 200). If the sectors were different, the polygons would be valid. The violet and green polygons are okay, because for the violet one, the values are all the same and for the green one, the polygons are not touching.

The reason why I am doing this is because the user must use only one polygon for a same sector if the polygons are touching. So if all the values are the same, I can dissolve the polygons otherwise the dissolve won't work.

Is there a way to do this?


texts: Global Mapper (improving KMZ quality):"raster_export-KMZ" batch script creates a poor quality .kmz compared to the "raster export - KMZ" manual comm

texts:
I'm sure it's something I'm missing...
In Global Mapper (GM) I've got my workspace loaded.. and I run a script that does nothing more than "EXPORT_RASTER" to KMZ. The orgin rasters loaded in GM will be cropped and saved in different separate .kmz(s) using a grid (polygon shapefile).
This is the script:

GLOBAL_MAPPER_SCRIPT VERSION=1.00
SET_LOG_FILE FILENAME="H:\temp\GlobalMapperScriptLogCapoluoghi.txt"
EXPORT_RASTER FILENAME="H:\Lavoro\_dati\Regioni\_Ortofoto_capoluoghi regione\_kmz\.kmz" \
TYPE=KML \
KML_RASTER_FORMAT=JPG \
KML_SUPER_OVERLAY=YES \
KML_TILE_SIZE=1024 \
POLYGON_CROP_FILE="H:\Lavoro\_dati\QU_ritagli_capoluoghi.shp" \
POLYGON_CROP_USE_EACH=YES \
POLYGON_CROP_NAME_ATTR="<Feature Name>"

Everything works fine but.... the .kmz(s) created with the batch-script are very poor in quality compared to the same identical kmz(s) exported using the manual File-->Export-->Export raster/image format-->KML/KMZ command.
i.e. the final kmz of the exactly same cropped area is 40mb exported manually on the contrary it's 24Mb using the script.
Moreover, the .kml(s) inside the two different .kmz(s) of the exactly same cropped area are very different (long and "rich" info-text in the manual export, very short and "poor" info-text in the batch script export - see image here attached).

I also tried to use these tags in the script without success:

KML_TILE_SIZE=1024 \
PIXEL_SIZE=1024x1024 \
SPATIAL_RES=1.79663056823901e-006,1.79663056823901e-006 \

What am I missing? What's wrong in my script?
I'm using GlobalMapper 19.


texts: How to measure similarity of SpatialLines objects

texts:
I created two SpatialLines objects in R: .

These objects were created this way:

library(sp)
xy <- cbind(x,y)
xy.sp = sp::SpatialPoints(xy)
spl1 <- sp::SpatialLines(list(Lines(Line(xy.sp), ID="a")))

Now I want to somehow conclude that this is the same line rotated and flipped, and that the difference between them is equal to 0 (i.e. shape is equal).

To do that, one can use maptools package and rotate line #1, e.g.:

spl180 <- maptools::elide(spl1, rotate=180)

Each rotated line must then be checked versus line #2 using rgeos package, e.g.:

hdist <- rgeos::gDistance(spl180, spl2, byid=FALSE, hausdorff=TRUE)

However, this is so computationally expensive way to match SpatialLines objects, especially if the number of objects is around 1000.

Is there any clever way to do this job?

P.S. Moreover, the above described approach does not guarantee all possible rotations and flips.

P.S2. If line #1 is zoomed out with respect to line #2, the the difference between line #1 and #2 must still be equal to 0.

UPDATE:


texts: Make lots of maps at once using QGIS print composer

texts:
I have around 1250 maps of species distributions.
Each species is a different shapefile projected on one base-map.

I can use the TomBio addon to export everything as a single map (PNG) without too much trouble, but I don't want seperate images. I want 8 small maps on a A4 format. Making this in the QGIS print composer is simple, but it takes a lot of time to produce all the maps. I was wondering if there is a faster way to automatically generate all my maps using QGIS.

So: Is there a way to use the template in my print composer (8 small maps of 8 different species on A4 format) to generate this for all of my 1250 species without having to do this all manually.


texts: Parallel Computing error using GP tool in C#

texts:
I have been using ArcGIS Geoprocessor(GP) tool whenever possible even though it is like a black box in my ArcObject development.
Most of my operations are involved with data processing, such as dataset projection, format conversion.
In order to take the advantage of my multi-core system, I decide to use parallel computing/multithreading in C# since it supports since framework 4.0:
http://msdn.microsoft.com/en-us/library/dd460720(v=vs.110).aspx
So I basically modified my loop which calls GP tool.
But it works for only a few GP operations, most of others failed with error handle.
So I want to know how to resolve this issue, is there anyway I can continue to use GP in parallel computing?


texts: Snapping to added guides in QGIS Layout does not actually align objects

texts:
I added some horizontal guides to QGIS layout with the idea of aligning my main map frame to the legend. I also activated the "Snap to Guides" option under View>Snap to Guides. When I snap the map frame with the legend, they are both snapped to the guide but appear not to be aligned:

I am using QGIS 3.2.3. Why does this occur? If the guide is not supposed to align objects on the map layout then what on earth is its purpose?


texts: GeoPandas: return lat and long of a centroid point

texts:
I would like to return the lat and lon or x and y from the point object in the following series as I plan to link them with an API I built that links OSM and NetworkX. The centroids will be added as new nodes for network analysis.

gp.GeoSeries(zones.centroid).x, and gp.GeoSeries(zones.centroid).y as outlined in docs raise AttributeError: 'GeoSeries' object has no attribute 'x' error.

Modifying things a bit and printing list(gp.GeoSeries(zones.centroid)) return thousands of shapely points of the following format:

[... <shapely.geometry.point.Point object at 0x0000000024035940>,
<shapely.geometry.point.Point object at 0x0000000024035978>,
<shapely.geometry.point.Point object at 0x00000000240359B0>,
<shapely.geometry.point.Point object at 0x00000000240359E8>,
<shapely.geometry.point.Point object at 0x0000000024035A20>,
<shapely.geometry.point.Point object at 0x0000000024035A58>,
<shapely.geometry.point.Point object at 0x0000000024035A90>,
<shapely.geometry.point.Point object at 0x0000000024035AC8>]

The code I'm using is the following:

import geopandas as gp

zones = gp.GeoDataFrame.from_file(shp_file)

for index, row in zones.iterrows():
    print index, gp.GeoSeries(zones.centroid)

# result:
# 9700022.00    POINT (-122.8196050489696 54.00617624128658)
# 9700023.00    POINT (-122.7474362519174 53.99998921974029)
# 9700100.00    POINT (-121.4904983300892 53.98447191612864)
# 9700101.00    POINT (-122.5513619751679 53.73999791511078)
# 9700102.00    POINT (-123.0624037191615 53.62317549646422)
# 9700103.00    POINT (-123.0848175548173 54.05921695782788)

How can I return the x and y from the GeoPandas POINT object?


texts: Installing PCL module for Python 3.6 in Anaconda

texts:
I am using Python 3.6 in Anaconda3 for Windows10(64) and would like to install PCL module. I tried conda install -c sirokujira python-pcl from (https://anaconda.org/sirokujira/python-pcl) and successfully installed
pcl: 1.8.1-hd76163c_1 conda-forge
python-pcl: 0.3.0-py36h0881630_0 sirokujira
However, import pcl in python gives error ImportError: DLL load failed: The specified procedure could not be found.


texts: QGIS geoprocessing script, ## hand-over processing-specific information error

texts:
I have a QGIS python geoprocessing script starting like this:

## Networking = group
## lIn = vector
## lOut = output vector

from qgis.core import *

v = processing.getObject(lIn)

#[...the rest...]

Executing it make a wonderful gui

but clicking 'Run' gives me an error

What am I probably missing?!


texts: Convert a .tif raster in .img format in QGIS

texts:
I have a raster file in .tif format and I want to convert it to .img with QGIS.

Does anyone know how I could do that?


texts:
I am trying to save my raster layer as follows:
file_writer = QgsRasterFileWriter('file_to_write.tif')
pipe = QgsRasterPipe()
rlayer = lu.getLayerByName('raster_name')
provider = rlayer.dataProvider()

if not pipe.set(provider.clone()):
    return
file_writer.writeRaster(
    pipe,
    provider.xSize(),
    provider.ySize(),
    provider.extent(),
    provider.crs())

The problem is that the output is a totally black image where I can barely see the content.
The script used to load the raster can be found below:
def load_raster(path, iface, opt=0):
   fileInfo = QFileInfo(path)
   baseName = fileInfo.baseName()
   layer = QgsRasterLayer(path, baseName)
   crs = layer.crs()
   crs.createFromId(2100)
   layer.setCrs(crs)
   raster_trans = layer.renderer().rasterTransparency()
   ltr = QgsRasterTransparency.TransparentSingleValuePixel()
   tr_list = []
   ltr.min = 150  # Or another value
   ltr.max = 255  # Or another value
   ltr.percentTransparent = 100  # Or another value
   tr_list.append(ltr)
   raster_trans.setTransparentSingleValuePixelList(tr_list)
   provider = layer.dataProvider()
   provider.setNoDataValue(1, 0)
   QgsProject.instance().addMapLayer(layer)

For some rasters, the above work like a charm but most of the times the output is total black. What could it be wrong?


texts: Is there a way to create a network based on road elevations , to build a sewer/Strom drain network . The idea is to draw a polyline from higher elevation to lower elevation automatically.


texts:
I have a point layer which symbolize trees. This trees has fields "Specie" and "Diameter" among others. I would like to symbolize them by size depending on the "Diameter" but also (at the same time) by colours depending on the "Specie".
So at the end I would have trees on different sizes but also colours. How can I achieve this?


texts: Display SVG icons in polygon layer using PyQGIS3

texts:
I have this code to set SVG icons only for selected features. It works, but only with point geometry, while I have a layer with polygon geometry. How can I modify this code to works also with polygon geometry?

vlayer = iface.activeLayer()

svgStyle = {}
svgStyle['name'] = path
svgStyle['outline'] = '#000000'
svgStyle['size'] = '15'
symbolLayer = QgsSvgMarkerSymbolLayer.create(svgStyle)

symbol = QgsSymbol.defaultSymbol(vlayer.geometryType())
symbol.changeSymbolLayer(0, symbolLayer)

renderer = QgsRuleBasedRenderer(symbol)
root_rule = renderer.rootRule()
rule = root_rule.children()[0].clone()
rule.setLabel('test')
rule.setFilterExpression('is_selected()')
root_rule.appendChild(rule)

root_rule.removeChildAt(0)
vlayer.setRenderer(renderer)

I thought that probably it's not possible to apply a SvgMarkerSymbol to a polygon layer and that I have to use before QgsCentroidFillSymbolLayer, but I don't know how.


texts: How to set values from list to row using arcpy

texts:
I have list with values and I want to put them in rows for specific column. Here is my code:

den = kol() # list with values
with arcpy.da.UpdateCursor('Wojewodztwa', 'gest_zal') as cursor:
    for row in cursor:
        for i in range(len(den)):
            row.setValue(den[i])

During execution I've got error: 'list' object has no attribute 'setValue'. I also try to iterate throught list and set row[0] = den[i] but all I got was that all values in rows are the same last number. How can I make it work?


texts: is it possible to use hillshade function in ArcMap to calculate shadow areas of buildings or simulate shadows of buildings, using DSM of the buildings and the sun angle and azimuth?


texts:
I am struggling to plot points using WGS84 co-ordinates on a map of the UK using the OSGB 1936 / British National Grid CRS. I'm using the rgdal package in R.

Here is my reproducible example, using the co-ordinates of Buckingham Palace...

library(rgdal)
foo <- data.frame(lat = 51.500829, long = -0.142994)
bar <- SpatialPoints(foo)
proj4string(bar) <- CRS("+init=epsg:4326")
bar <- spTransform(bar, CRS("+init=epsg:27700"))
coordinates(bar)

This outputs the following transformed co-ordinates:

         lat     long
[1,] 7464126 -5554500

Where I should be getting latitude of 179622 and longitude of 528993.

Can anyone tell me what I'm doing wrong? The only similar problem I could find was here and I'm now pretty confident I am choosing the right CRSs. I'm trying to do something similar to this post but again, I don't think the solutions there can help me.

P.S. I'm using this site to check what co-ordinates I should be getting.


texts: gdal2tiles doesn't use resampling method correct?

texts:
I'm using gdal2tiles to create tiles from MODIS satellite imagery (250 m resolution) and I have noticed that I get the same bad result if I use -r near or -r bilinear (or lanczos for that matter). See attached image.

If I do the same generating with Mapnik and generate_tiles.py the images looks so much better.

The source file has a smaller resolution than the output tile, but I guess that's what the resampling bilinear (or lanczos) should "take care of"? Now it seems like gdal2tiles only scales the image without adding any resampling method?

I've tried to generate the tiles from geotiffs and jpegs (with jgw) but it's all the same.

The reason I want to use gdal2tiles rather than generate_tiles.py with Mapnik is that gdal2tiles are so much faster. Look here: The fastest way of generating tiles from raster?


texts:
I'm relatively new at working with the CMV (and .css).

What's the easiest way to change the Home button's position?

To clarify, I'm attempting to move the position of the button itself down about 100px.


texts: QTableWidget Resize with Dialogue Python QGIS

texts:
I am developing python qgis plugin with QTableWidget inside a tabWidget. How can i set my QTableWidget size and QTabwidget with the Dialogue box. So that when someone expand the dialogue box of the plugin from QGIS window,it won't show a small window of the Attribute table. First image is showing original window.But after changing dialogue size it becomes like second image

second image:


texts:
In the past I have used OL extensively in order to make WFS requests to Geoserver and add layer on a map.
Now I am trying to learn LeafletJS and I want to do something similar.
I have made a WMS request successfully but I can not find documentation of how to do WFS request.
This is the WMS request to geoserver:

var layer_name = L.tileLayer.wms("http://something:xxxx/geoserver/wms", {
        layers: 'workspace:geoserver_layer',
        format: 'image/png',
        transparent: true,
        version: '1.1.0',
    });
    layer_name.addTo(map);

Also, how do I solve the problem with the proxy. I remember in OL I used a cgi script.


texts: QGIS 3.0 point cluster renderer, is it possible to control rendering order (z-order)?

texts:
I've been having a play with the new point cluster renderer in QGIS 3.0.2 and it's great - a lot easier (and faster) than having to do the clustering yourself.

To show what I mean, here's an example of road accidents in Edinburgh using STATS19 data and a 30m radius.

  I can't work out how to set up the z-order so that the largest clusters are drawn in front of smaller clusters.

I'm aware of the @cluster_size variable - it's what I'm using to scale and colour the symbols. For this I'm using the size and colour assistants. I removed the font markers to remove the text from each cluster.

Size of markers is defined as follows:-

coalesce(scale_exp(@cluster_size, 0, 50, 1, 200, 0.57), 0)

Colour of markers is defined as follows:-

coalesce(ramp_color('custom ramp',scale_linear( @cluster_size, 0, 50, 0, 1), '#000000')

What I've tried so far:-

Symbol Levels - No luck; unable to select a variable name.
Layer rendering > Control Feature Rendering Order and choose @cluster_size as sort. Doesn't appear to have an effect.
Graduated style using @cluster_size as the expression doesn't work.

  Does anyone know if this is possible, or is this a limitation with this particular renderer?

I know there's a lot of other ways to do this (buffer/spatial join, plugins), but the new renderer should be a lot faster. Or have I overlooked something obvious? :)


texts: I have requirement where i have to use web map . I have have lots of feature layers which will loaded conditionally hence in my json i have provided visibility as false for all layers then i have modify json conditionally and set visibility to true for some layers. The problem is once we set visibility to true , all the features will be loaded instead i want to set Definition or pass some query task where i can load features based on need.


texts: Cropping a grid of polygons to a single polygon gives an odd result

texts:
I am trying to create a grid over a polygon from which I can choose a random subset of cells for sampling. The problem occurs when I crop the grid, which I am converting to polygons, by the polygon shapefile. For some reason this results in a polygon the shape of the cropping polygon but with lines across it. After opening the resulting shapefile in QGIS I can see that the lines are "edges" of the polygon, because when I get rid of the outline of the polygon, the lines go away.

Can you explain what error, if any, I am making, or do you have an alternate means of generating a grid over a polygon from which a subset of grid cells can be selected (and ideally turned in a shapefile)?

Here is some sample code:

(California boundary found at https://data.ca.gov/dataset/ca-geographic-boundaries)

require(rgdal)
require(sp)
require(raster)
require(rgeos)

# Load shapefile
ca <- readOGR(dsn = "...", layer = "CA_State_TIGER2016")
ca.utm <- spTransform(ca, CRSobj = "+proj=utm +zone=11 ellps=WGS84")

# Calculate number of grid cells necessary
x.cells <- ceiling((ca.utm@bbox[1,2]-ca.utm@bbox[1,1])/10000)
y.cells <- ceiling((ca.utm@bbox[2,2]-ca.utm@bbox[2,1])/10000)

# Generate grid
ca.grid <- raster(extent(ca.utm), nrow=y.cells,
                      ncol=x.cells, crs=ca.utm@proj4string)
ca.grid.poly <- as(as(ca.grid, "SpatialPixels"),"SpatialPolygons")
ca.grid.utm <- spTransform(ca.grid.poly, CRSobj = "+proj=utm +zone=11 ellps=WGS84")

# Crop grid to polygon
ca.grid.inter <- gIntersection(ca.grid.utm, ca.utm)
# or
ca.grid.crop <- crop(ca.grid.utm, ca.utm)

# both gIntersection and crop give same result as far as I can tell

Here is a screen shot:

plot(ca.grid.inter)


texts: Leaflet popup does not work with esri rest feature services

texts:
I am writing a code that basically loads an arcgis service layer using leaflet.The code does as below:

Loads POINT layer from esri rest feature server and Categorizes
points based on attribute named Zone.
Binds popup to each of these points using the contents of that
features service layer.
Loads another esri rest map service POLYGON layer and symbolize it based on attribute named Zone.

Now my problem is step 2 is not working!

Mycode:

  <meta charset=utf-8 />
  <title>Custom popups</title>
  <meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

  <link rel="stylesheet" href="https://unpkg.com/leaflet@1.0.2/dist/leaflet.css" />

    <link rel="stylesheet" href="https://unpkg.com/esri-leaflet-geocoder@2.2.3/dist/esri-leaflet-geocoder.css">

N.B. I followed http://leafletjs.com/examples/geojson/example.html


texts: GeoServer REST API: how to create an MBTiles datastore

texts:
I want to create an MBTiles datastore using the REST API of GeoServer.

I was able to create a datastore for a shapefile using:

curl -v -u admin:geoserver -XPUT -H "Content-type: text/plain"
  -d "file:///data/shapefiles/rivers/rivers.shp"
  http://localhost:8080/geoserver/rest/workspaces/acme/datastores/rivers/external.shp

I didn't find any example for MBTiles files. So, I tried the following (after adding the MBTiles extension to my GeoServer):

curl -v -u admin:geoserver -XPUT -H "Content-type: text/plain"
  -d "file:///data/path/to/mbtiles/file.mbtiles"
  http://localhost:8080/geoserver/rest/workspaces/acme/datastores/rivers/external.mbtiles

But got the following return code:

400 Unsupported format: mbtiles

Can any body help me with the right way or a turn around way.


texts: Using Calculate Geometry tool to calculate area when disabled?

texts:
I have added a field in my attribute table but I cannot calculate the area. The units field is greyed out and the options at the top say "area - disabled".

What can I do about this?


texts: How to edit attribute table with decimal numbers in QGIS

texts:
I have been adding points to a vector layer and I want to define their value in a specific column. When I try to edit the values, that are in decimal numbers, it only allows me to type integers.

How do I edit points in my attribute table in decimal numbers?

This is before(Z_COORDS) and after(ZCOORD) I tried to edit the column.


texts: Surface area of waterbodies for multiple location using Google Earth Engine

texts:
Using this code, I've extracted surface area of water bodies (for a period of 6 years) from sentinel 1 GRD dataset. But how can I compute water area extent of multiple location and produce chart/CSV for all of them at once? I have a shapefile of 50 locations and I have created a ROI for each but adding them one by one is tedious. any workaround for this problem?
Code link: https://code.earthengine.google.com/24e31dff3a145f38bbfd7deb59636da5
// Load region defined by polygon and add it to the map

Map.addLayer(roi, {}, 'ROI')
Map.centerObject(roi, 12)

//Load Sentinel-1 SAR collection and filter according to data collection type
var S1 = ee.ImageCollection('COPERNICUS/S1_GRD')
  .filterBounds(roi)
  .filterDate('2014-01-01','2020-07-29')
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))

//Add first image to map to get an idea of what a SAR image looks like
Map.addLayer(S1.first(),{bands: 'VV',min: -18, max: 0}, 'SAR image')

// Filter speckle noise
var filterSpeckles = function(img) {
  var vv = img.select('VV') //select the VV polarization band
  var vv_smoothed = vv.focal_median(100,'circle','meters').rename('VV_Filtered') //Apply a focal median filter
  return img.addBands(vv_smoothed) // Add filtered VV band to original image
}

// Map speckle noise filter across collection. Result is same collection, with smoothed VV band added to each image
S1 = S1.map(filterSpeckles)

//Add speckle filtered image to map to sompare with raw SAR image
Map.addLayer(S1.first(),{bands: 'VV_Filtered',min: -18, max: 0}, 'Filtered SAR image')

//Classify water pixels using a set threshhold
//Here we are using -16. This is only an approximation and will result in some errors. Try adjusting the
var classifyWater = function(img) {
  var vv = img.select('VV_Filtered')
  var water = vv.lt(-16).rename('Water')  //Identify all pixels below threshold and set them equal to 1. All other pixels set to 0
  water = water.updateMask(water) //Remove all pixels equal to 0
  return img.addBands(water)  //Return image with added classified water band
}

//Map classification across sentinel-1 collection and print to console to inspect
S1 = S1.map(classifyWater)
print(S1)

//Make time series of water pixels within region
var ClassChart = ui.Chart.image.series({
  imageCollection: S1.select('Water'),
  region: roi,
  reducer: ee.Reducer.sum(),
  scale: 100,
})
  .setOptions({
      title: 'Inundated Pixels',
      hAxis: {'title': 'Date'},
      vAxis: {'title': 'Number of Inundated Pixels'},
      lineWidth: 2
    })

//Set the postion of the chart and add it to the map
ClassChart.style().set({
    position: 'bottom-right',
    width: '500px',
    height: '300px'
  });

Map.add(ClassChart)

// Create a label on the map.
var label = ui.Label('Click a point on the chart to show the image for that date.');
Map.add(label);

//Create callbakc function that adds image to the map coresponding with clicked data point on chart
ClassChart.onClick(function(xValue, yValue, seriesName) {
    if (!xValue) return;  // Selection was cleared.

    // Show the image for the clicked date.
    var equalDate = ee.Filter.equals('system:time_start', xValue);
    //Find image coresponding with clicked data and clip water classification to roi
    var classification = ee.Image(S1.filter(equalDate).first()).clip(roi).select('Water');
    var SARimage = ee.Image(S1.filter(equalDate).first());
    //Make map layer based on SAR image, reset the map layers, and add this new layer
    var S1Layer = ui.Map.Layer(SARimage, {
      bands: ['VV'],
      max: 0,
      min: -20
    });
    Map.layers().reset([S1Layer]);
    var visParams = {
      min: 0,
      max: 1,
      palette: ['#FFFFFF','#0000FF']
    }
    //Add water classification on top of SAR image
    Map.addLayer(classification,visParams,'Water')

    // Show a label with the date on the map.
    label.setValue((new Date(xValue)).toUTCString());
  });