texts: Masking WMS (Esri Oceans Basemap) layer in ArcGIS for Desktop?

texts:
I would like to display ESRI's Oceans basemap web map service (WMS) in one of my maps, but not the whole layer, just the landmass depicting Canada as I find it a nice fit for the look I'm trying to achieve. Is there any way to mask the features around Canada's landmass in that layer? I tried using the Advanced Drawing feature, but the Oceans basemap layer is not listed, probably because it's a WMS. I am working with ArcGIS 10.2. Are there any suggestions as to how I could do this?

Otherwise, if anyone can suggest a symbolized layer of Canada's landmass (and only containing Canada's landmass) ressembling that depicted on the ESRI Oceans Basemap that I could download or connect to somewhere, I would be grateful.

I could make my own. I tried with a few DEM's and hillshades but didn't get good results.


texts: PostGIS: Approximating a house number from address range

texts:
I'm developing a reverse-geocoder for Canada. So far, given a lat/lng, I can find the nearest street (line segment), which includes line segment direction and address ranges for both sides of the street. I'm now trying to figure out the best way to programmatically approximate the nearest house number to the given lat/lng point.

Here's an example of a row containing the street data:

-[ RECORD 1 ]-
[...]
l_adddirfg | Same Direction
l_hnumf    | 3219
l_hnuml    | 3235
l_stname_c | Breen Road North-west
r_adddirfg | Same Direction
r_hnumf    | 3224
r_hnuml    | 3236
r_stname_c | Breen Road North-west
the_geom   | 0105000020E610000001000000010200000002000000B0F6990E78885CC088DF2B5F3C8C49400875B39A89885CC0A0BCA6AC4B8C4940

So, given a lat/lng coordinate that lies near the "the_geom" line segment, a person could tell visually which side of the street the point is on (left or right side), and how far along the segment it is -- thereby approximating a house number. For example, if the point lies on the right side, three-quarters down the street, I would use the fields r_hnumf (right side, first number) and r_hnuml (right side, last number)...  The street address is probably close to:

3232 Breen Road North-west

What I'm looking for is a best practice in either computing/approximating this in PostGIS (which I'm new to), or in the application layer once the row is fetched.

Thanks!


texts:
I can no longer search for tools, etc using the Search window in ArcMap i.e. Windows | Search doesn't work and nor does ctrl + f.

I thought updating my student license would do the trick but after doing so it's still the same. I have also tried to fix it by deleting Normal.mxt but that didn't work either.

I read somewhere that ArcMap can sometimes misbehave because of problems related to the pc user account.

I may accidentally have shut down my computer the wrong way last week and before I was properly logged off (by pulling the plug too early).

What could be the cause of this and how may I fix it?

My version of ArcMap is 10.5.0.6491.


texts: Are the images produced from the "Export mosaic dataset items" orthorectified?

texts: I am trying to find whether or not the orthorectification is done to these images like it is for the export data, I can't find anything about this.


texts: Import multiple PostGIS layers in QGIS map canvas via QTask QGIS 3

texts:
The vertex tool in red and the geometry displayed by the layer in yellow are not in sync

from time import sleep
from qgis.PyQt.QtSql import QSqlQuery, QSqlDatabase
serverName = 'hostname'
database = 'dbname'
port = 'port'
usr = 'username'
pw = 'password'
uri = QgsDataSourceUri()
uri.setConnection(serverName, port, database, usr, pw)

def stopped(task):
    pass
def completed(exception, result=None):
    pass
def run(task, wait_time):
    # main code
    wait_time = wait_time / 100

    Layers= ['120_pipes']

    for i, lyr in enumerate(Layers):
        task.setProgress(i)
        sleep(wait_time)
        uri.setDataSource('120', lyr, "wkb_geometry")
        vlayer = QgsVectorLayer(uri.uri(), lyr, "postgres")
        if vlayer.isValid():
            g = QgsProject.instance().layerTreeRoot().insertGroup(0, '120')
            QgsProject.instance().addMapLayer(vlayer, False)
            g.insertChildNode(0, QgsLayerTreeLayer(vlayer))

        if task.isCanceled():
            stopped(task)
            return None

task = QgsTask.fromFunction(u'Task', run,
                            on_finished=completed, wait_time=50)
QgsApplication.taskManager().addTask(task)


texts: QGIS - Merge a polygon line that is splitted

texts:
Verry much a beginner at this but i have some data that contains trails that are split into smaller sections, id like to be able to merge those small sections into longer trailes based on the close proximity between the parts, how to do this? In total i have about 4000 of those small lines that i need to merge together to make about a total of 100 longer lines/trails.

Above is an example where i have highlighted an example of this.
I have read something about using buffer at 0 width but that does not seme to work for me, and when i use a higher buffer width im getting buffered with doubble the amount of points.
I have also worked with dissolve but that gives me all of thracks joind together into one big entity, i still want them to be split apart, so in the image id like to end up with 4 different tracks.


texts:
I have been playing with the ST_Distance function available within PostGIS and to my surprise, I have encountered some rather strange behavior...

In the simplest setting, I have just two tables N and P. Table N contains point coordinates, while table P contains multipolygons representing administrative divisions. Now, I would like to calculate the distance between a particular node and given range of (multi)polygons. To this end, I just invoked following command:

demo=> SELECT N.node_id, P.iid, ST_Distance(P.outline::geography, N.location::geography, false)
       FROM N, P WHERE (N.node_id=9 AND P.iid=25);

  node_id |  iid  |   st_distance
 ---------+-------+------------------
        9 |    25 | 1732718.95934996
 (1 row)

which seems to produce the expected result.
However, if I slightly modify the WHERE condition:

demo=> SELECT N.node_id, P.iid, ST_Distance(P.outline::geography, N.location::geography, false)
       FROM N, P WHERE (N.node_id=9 AND P.iid>=20 AND P.iid<=25);
 node_id |  iid  |   st_distance
---------+-------+------------------
       9 |    20 | 1738836.97788419
       9 |    21 | 1723001.48856517
       9 |    22 | 1736600.89344563
       9 |    23 | 1718786.25166805
       9 |    24 | 1503840.93486047
       9 |    25 |                0
(6 rows)

the result for pair 9, 25 suddenly changes. Am I missing something obvious? I am using PostgreSQL 9.4.5 (Postgis 2.1.8) on CentOS 7.

EDIT:

It seems that the problem persists even if I test only the "suspicious" polygon versus several points, i.e., let's consider only following points

 node_id |       st_astext
---------+------------------------
       8 | POINT(-3.9018 42.8116)
       9 | POINT(-6.1947 41.5758)
      10 | POINT(-6.0858 41.4912)

and only polygon 25 (http://pastebin.com/7deJPC3u). Then I get:

demo=> SELECT N.node_id, P.iid, ST_Distance(P.outline::geography, N.location::geography, false)
       FROM N, P WHERE (N.node_id IN (9) AND P.iid=25);
 node_id |  iid  |   st_distance
---------+-------+------------------
       9 |    25 | 1732718.95934996
(1 row)

as well as:

demo=> SELECT N.node_id, P.iid, ST_Distance(P.outline::geography, N.location::geography, false)
       FROM N, P WHERE (N.node_id IN (8,9,10) AND P.iid=25);
 node_id |  iid  |   st_distance
---------+-------+------------------
       8 |    25 | 1731869.60412558
       9 |    25 |                0
      10 |    25 | 1729708.79737747
(3 rows)


texts: How to call the Manage Map Cache Tiles tool hosted by default in the system folder of arcgis server. I am getting ERROR 001428

texts:
When I am trying to call the GP tool Manage map service cache tiles that is hosted on the arcgis server by default in the system folder, I am geeting the followning error:
ERROR 001428 Failed to retrieve configuration for service .

Here is the message of the call I am trying to perform:
esriJobMessageTypeInformative: Executing (Manage Map Cache Tiles): ManageMapCacheTiles http://go-infra.eurodyn.com:6080/arcgis/rest/services/OSCAR/OSCAR_GO_NMI_Cached_service_program_All/MapServer 591657527,591555;295828763,795777;147914381,897889;73957190,948944;36978595,474472;18489297,737236;9244648,868618;4622324,434309;2311162,217155 3 RECREATE_ALL_TILES # # "Feature Set"

and here the complete error message:
esriJobMessageTypeError: ERROR 001428: Failed to retrieve configuration for service http://go-infra.eurodyn.com:6080/arcgis/rest/services/OSCAR/OSCAR_GO_NMI_Cached_service_program_All/MapServer.
esriJobMessageTypeError: Failed to execute (Manage Map Cache Tiles).
esriJobMessageTypeInformative: Failed at Wed May 27 19:34:16 2015 (Elapsed Time: 0,16 seconds)
esriJobMessageTypeError: Failed.


texts: Coloring only half raster polygon in QGIS?

texts:
I'm wondering whether it's possible to fill only half of a polygon with a color so I can add another layer in which I color the other half and double the amount of information presented.

For example:

For a certain area I have 2 datasets. One with plant species seen this year and one with plant species seen 10 years ago. Adding a grid layer and doing a unique point count vector analysis will give me the amount of species per grid square (from this year, 10 years ago and by subtracting I can see whether the amount increased or decreased).

Now to present this data: Ideally I want to duplicate the layer. In one layer I want to color only the lower-left diagonal half of each polygon, where I present the data from last year. In the other layer I want to fill the top-right diagonal half with a different color, representing the change over 10 years.

A bit like this

My grid squares are 100 by 100 meters. I tried to fill the polygons with only diagonal lines setting the line thickness to half the square as well as the distance between the lines, but that became messy.


texts: How do I add fields to a layer attribute table .created by CSV import into QGIS?

texts:
I have added a layer from a delimited text (csv) file.  I now wish to add the usual X, Y and Rotation fields to allow me to move and rotate the layer labels.  However, it doesn't seem possible to edit the attribute table as the edit icon is greyed out for this table.

Can I somehow overcome this problem?


texts:
Does anyone has idea how to deal with this error:

  Warning 1: organizePolygons() recieved an unexpected geometry. Either a polygon with interior rings or a polygon with less than 4 points or a non-Polygon geometry. Retrun arguments as collection.

What does it mean 'retrun arguments as collection?

I am using data in osm format. Error occurs for some polygons while using:

temp_path = path_to_folder
element = name_of_file
path = path_for_element

r'ogr2ogr {0} {1}.osm multipolygons -nln {2}_ply -lco ENCODING=UTF-8'.format(
        temp_path, os.path.join(path, element), element)


texts: Rotation field is missing in rule-based renderer in QGIS 2.6

texts:
I am trying to do rotation of point symbols in QGIS. I am working with rule-based renderer. In other types of classification there is a button advanced (for example in categorized renderer), where I can choose rotation field.
But in rule-based renderer this button is missing.

Do you know how to make this rotation of point symbols otherwise?

Example in categorized renderer

In rule-based renderer advanced(and rotation field) is missing


texts: What is the link between proj4js and proj

texts:
Question

I am simply wondering, what is the relationship between the web oriented proj4js and the PROJ C++ well known library (on which for example GDAL/OGR relies)?

Is there any link between these two libraries?

It's written on the proj4js web page that it is:

  Originally a port of PROJ.4 and GCTCP C it is a part of the MetaCRS group of projects.

But what is the current situation?
How could I check what version of PROJ does proj4js rely on?

Sources repositories of each of these projets:

https://github.com/proj4js/proj4js
https://github.com/OSGeo/PROJ


texts: Do I need to perform any corrections on Level-2 Landsat imagery?

texts:
My objective is to perform a land cover classification with multiple scenes over historical imagery (one year every decade).

Most of the imagery I have downloaded is level-2 imagery. From my understanding level-2 imagery is already atmospherically corrected. But since I am using multiple scenes, does that mean I still have to manually perform corrections even though it is a level-2 product? If so, what are all the types of correction that I need? Atmospheric, solar, geometric?

Edit: When I mosaic the images I get a huge visual difference, even though the images are only taken days apart and I have mosaicked them so their stretch hand display settings should be the same.


texts: Postgis bounding box intersect vs polygon intersect

texts: I am running a point in polygon test for close to one billion points. The size of my polygons are very small (at about 12 vertices each polygon). Therefore, I would be satisfied even with a bounding box intersect test. However, when I checked,point_geom && polygon_geom (bounding box test) and ST_Contains(point_geom, polygon_geom) (polygon contains test) yield the same amount of time. Rather, I expected the bounding box test to be much faster than the polygon contains test.


texts: Error reading large virtual raster with rasterio: "Cannot invert degenerate transform"

texts:
I've created a virtual raster using gdalbuildvrt, which covers a large extent (England) at a high resolution. The data is fairly sparse, focused around the coastlines with no files covering inland areas. I'm trying to read the data in Python using rasterio.

When I try to invert the affine transformation (to get world coordinates converted to pixel coordinates) I getting the following error:

>>> layer.affine
Affine(0.0002777778000000001, 0.0, -5.5000000044,
       0.0, -0.0002777778000000001, 55.0000000836)
>>> layer.affine.__invert__()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "c:\Python27\lib\site-packages\affine-1.0-py2.7.egg\affine\__init__.py", line 407, in __invert__
    "Cannot invert degenerate transform")
affine.TransformNotInvertibleError: Cannot invert degenerate transform

I don't get this error accessing the individual raster files. I have some idea of what it means, but not how to work around it. Perhaps it needs approaching in a different way? How do I access this data without having to worry about which file it's stored in?


texts: QGIS's Difference function inconsistently applied in complicated shapefile

texts:
In QGIS, I'm using Difference (under the Geoprocessing menu) to subtract a complicated shapefile with plenty of features from another complicated shapefile.

For 90% of the features, it works correctly. But for some, it ends up deleting the underlying feature altogether, instead of returning the difference -- even though the "Difference" feature didn't entirely cover up the underlying feature.

Here's an example:

The selected feature on the base layer.

The soon-to-be-Difference layer, in pink. Note that the selected feature is not wholly obscured by the difference layer, so it should show up in the final result.

The resulting shapefile. Note that the previously selected area is entirely blank (along with a good deal of the map that shouldn't be, as well.)


texts: How can I generate batch of coordinates from linestring to simulate driving on map?

texts:
I have a openlayers map and I want to simulate driving a route on it.
I would like to draw a line (using GPX tool) and create gps coordinates (preferred with advancing time stamp) so I can query them based on time and show on the map the advancement.

what is the easiest way doing this interpolation ? I have PostGIS and QGIS (or other tool !?) but I'm not familiar with python programming for QGIS.


texts:
I´m very new to QGIS and I would just look to have a "screenshot" of a map that I created.

When using snipping tool I get images of rather bad quality and when I use the "create ouput" tool my image is very large and it seems that this tool is mostly used to add other elements to the maps.

How do I get simple "screenshots" of a maps in decent quality?


texts: Iterating spatial join with multiple merge rules (and then creating raster) in ModelBuilder?

texts:
I'm using a Basic licence in ArcGIS Desktop 10.3 with Spatial & 3D Analyst extensions.

I have a dbf file which contains point data for several days.  I want my final output to be

  a) a point feature class for each day.

  b) A polygon feature class for each day with stats (e.g.: minimum, mean and maximum) for those points within a cell in a hexagon
  grid.

  c) A raster based on the mean value in each cell for those days.

I've got as far as having the individual point feature classes for each day (using a couple of python scripts within a model) but I am stuck on getting the iteration right for the Spatial Join and Rasters to be created.

I can't see how I can tell the model that I want the name(s) to be in a certain format.

For the Feature Classes that result from the Spatial Join I want to use the original input Feature Class name with "_SJHexGrid" appended to it.

Then for the rasters that are generated from those "_SJHexGrid" files I would want to append with "_MeanRast".

This is a diagram of the model as it currently stands (the "(wildcard)" in the Name input is just to show I am using a wildcard for illustration purposes here).


texts: I'm running QGIS 2.1.4 and Windows 7. For each layer, under properties there is a drop down box called 'data source encoding'. Mine is set to 'system'. I am wondering what this means, there seems to be a lot of options one can choose from.


texts: How to multiply a user-defined input parameter with a field value within the graphical modeler?

texts:
I have the same Question as asked in this post (How to connect input parameter to the field calculator algorithm within the graphical modeler?) and the answer from Joseph already helped a lot:

##Update field by number=name
##Layer=vector
##Fields=Field Layer
##Number=number 0

from qgis.core import QgsExpression

layer = processing.getObject(Layer)
layer.startEditing()

num = str(Number)
idx = layer.fieldNameIndex(Fields)
e = QgsExpression(num + '* 5')    # Change value to suit your needs
e.prepare(layer.pendingFields())

for f in layer.getFeatures():
    f[idx] = e.evaluate(f)
    layer.updateFeature(f)

layer.commitChanges()

Anyway, I'm looking for a way to multiply the number that the user enters in the interface with the value of the chosen field. In the script provided in the answer, the field value is ignored and the user value gets multiplied with the number defined in the script (default 5).

I have to admit that I'm not that much of a python programmer, so any suggestion would be very helpful.


texts: Creating a heatmap from TIFF file in QGIS

texts:
The task which I am currently working on, has us to create a kind of a 'heatmap' with KDE (or similar).

As a first step I vectorized the TIFF File and now it looks like this:

The end product should look something like a heatmap (ie. not chopped up like it currently is).
Is this possible in QGIS? Or is there a Python code/library I could use that does the trick?


texts:
I have many vehicles GPS data points, which consist of "timestamp, longitude, latitude, fuel consumption and distance" and the data is OSM(open street map) based. I want to match the data points to OSM, divide it by road type, so that I could get road type specific data in groups (e.g, highway group, local roadway group). Then I could calculate the fuel efficiency corresponding to different road type.

I'm wondering if R has any packages could handle this type of task?


texts: Using gdal_rasterize with a MSSQL Connection, possible?

texts:
I'm a bit new to gdal, I mostly use ogr2ogr for data imports to MS SQL Server -- so I'm a bit out of my tree on gdal_rasterize

My goal is to create png raster files 200x200 from a geometry column.  I tried something like with no avail.

gdal_rasterize -src_datasource "MSSQL:server={server};database=
     {database};uid={user};pwd={pwd};" -a geom -l layer1 -of png -burn
     255 -b geom -sql "SELECT top 1 geom from {myTable}" tract.png -of png -ts 200 200


texts: I have an issue with qfield. Is there any way to prevent creating duplicated data at each synchronisation with qgis. I would be very pleased if I get pertinent help.


texts: Dissolve with stats plugin does not appear in QGIS 3.0

texts:
So far we have tried to:

find it in the plugin repository from QGIS, with no luck.
install the plugin downloading it from the web QGIS Python Plugins
Repository: https://plugins.qgis.org/plugins/DissolveWithStats/
but it does not work.

The message that appears when trying to install the recently downloaded .zip is:

  No se pudo cargar el complemento 'DissolveWithStats' debido a un error
  al llamar a su método classFactory()
  ModuleNotFoundError: No module named 'dissolve_stats'
  Traceback (most recent call last):

  File "C:/PROGRA~1/QGIS3~1.0/apps/qgis/./python\qgis\utils.py", line
  336, in startPlugin
      plugins[packageName] = package.classFactory(iface)

  File "C:/Users/Rocio/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\DissolveWithStats__init__.py",
  line 25, in classFactory
      from dissolve_stats import DissolveWithStats

  File "C:/PROGRA~1/QGIS3~1.0/apps/qgis/./python\qgis\utils.py", line 664, in
  _import
      mod = _builtin_import(name, globals, locals, fromlist, level) ModuleNotFoundError: No module named 'dissolve_stats'

  Versión de Python: 3.6.0 (v3.6.0:41df79263a11, Dec 23 2016, 08:06:12)
  [MSC v.1900 64 bit (AMD64)]
  Versión de QGIS: 3.0.0-Girona Girona, 001c80b0c3

  Ruta de Python:

  C:/PROGRA~1/QGIS3~1.0/apps/qgis/./python
  C:/Users/Rocio/AppData/Roaming/QGIS/QGIS3\profiles\default/python
  C:/Users/Rocio/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins
  C:/PROGRA~1/QGIS3~1.0/apps/qgis/./python/plugins
  C:\PROGRA~1\QGIS3~1.0\bin\python36.zip
  C:\PROGRA~1\QGIS3~1.0\apps\Python36\DLLs
  C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib
  C:\PROGRA~1\QGIS3~1.0\bin
  C:\PROGRA~1\QGIS3~1.0\apps\Python36
  C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib\site-packages
  C:/Users/[...]/AppData/Roaming/QGIS/QGIS3\profiles\default/python
  Y:/Proyectos/[...]/Work/Graphics-Maps/Proyectos


texts: Automatically exporting time series of multiple polygons to .csv

texts:
I want to export data to a .csv file without having to display the chart. Similar to what is done here.

I want to skip the chart-display part because I'm working with long times series and if I try to display the chart I get the following error:"Error generating chart: Collection query aborted after accumulating over 5000 elements" as it is shown in the code below.

var regions = ee.FeatureCollection([
  ee.Feature(    // San Francisco1.
    ee.Geometry.Rectangle(-122.79, 37.82, -122.6, 37.6), {label: 'SF1'}),
  ee.Feature(  // San Francisco2.
    ee.Geometry.Rectangle(-122.6,  37, -122.48, 36.7), {label: 'SF2'}),
 ]);
Map.addLayer(regions);

var chl = ee.ImageCollection("NASA/OCEANDATA/MODIS-Aqua/L3SMI")
    .filterDate('2003-10-25', '2017-08-25')
    .select('chlor_a');

// Create a time series chart.
var chltimeseries = ui.Chart.image.seriesByRegion(
    chl, regions, ee.Reducer.mean(), 'chlor_a', 200, 'system:time_start', 'label')
        .setChartType('ScatterChart')
        .setOptions({
          title: 'Chlorophyl levels',
          vAxis: {title: 'mg/m3'},
          lineWidth: 1,
          pointSize: 4,
          series: {
            0: {color: 'FF0000'}, //SF1
            1: {color: '00FF00'}, //SF2
}});

// Display.
print(chltimeseries);


texts: My QGIS says "PyQGIS couldn't be loaded"

texts:
My QGIS (Version 3.2.3-Bonn 'Bonn', 9b17680) says:

  PyQGIS konnte nicht geladen werden. Die Python-Unterstützung wird
  abgeschaltet.

  Traceback (most recent call last):   File "", line 1, in
  ModuleNotFoundError: No module named 'qgis'

  Python-Version:
  3.6.6 (default, Sep 12 2018, 18:26:19)  [GCC 8.0.1 20180414 (experimental) [trunk revision 259383]]

  QGIS-Version:
  3.2.3-Bonn 'Bonn', 9b17680

  Python-Pfad: ['/usr/share/qgis/python',
  '/home/christoph/.local/share/QGIS/QGIS3/profiles/default/python',
  '/home/christoph/.local/share/QGIS/QGIS3/profiles/default/python/plugins',
  '/usr/share/qgis/python/plugins',
  '/home/christoph/.local/share/QGIS/QGIS3/profiles/default/python',
  '/home/christoph/anaconda3/lib/python36.zip',
  '/home/christoph/anaconda3/lib/python3.6',
  '/home/christoph/anaconda3/lib/python3.6/lib-dynload',
  '/home/christoph/anaconda3/lib/python3.6/site-packages']

Could it be, it's not finding the right Python?


texts: How to set up a Coverage Layer for Atlas Plugin?

texts: Has anybody utilized the Atlas plugin? I am starting to take a look and setting it up on a project. My question is how do you set up the layer that is utilized for the coverage layer. Is there a way to get a polygon the mimics the view of the map area shown in the Print Composer?


texts: Time Manager single point for time frame

texts: I'm using time manager for QGIS to create an animation. I am animating locations of buses. Their locations were sampled between every 2-5 minutes. So if I have a time window of 5 minutes, it will sometimes show two or more points for a single bus. I'm basically trying to figure out how to show only one point for a vehicle_id per time frame. Additionally when I interpolating points I lose the attributes so can't seem to style interpolated data.


texts: How to use 'memory:name' output as VectorLayer in PyQGIS?

texts:
After running processins in a loop, I try to export results on a single shapefile.
I've created a vector layer called "regroupee", I merge results with it, and I overwrite "regroupee" with the result.

As Kadir Şahbaz said in this post on the option 8, run() method with memory output make a dictionnary, with an instance of a layer (QgsVectorLayer).

The loop works, but when I try to export the final result with QgsVectorFileWriter.writeAsVectorFormat, an error say me that "regroupee" is a string...

My code, in the end of the for loop :

        alg_params = {
        'CRS': None,
        'LAYERS': [regroupee, calculatrice['OUTPUT']],
        'OUTPUT': 'memory:merge'
        }

        merge = processing.run('native:mergevectorlayers', alg_params, context=context, is_child_algorithm=True)
        regroupee = merge['OUTPUT']

    QgsVectorFileWriter.writeAsVectorFormat(regroupee, "I:/moulin_z/DONNEES/temp_us235/test_iteration_pyqgis/regroupee.shp", "UTF-8", crs.createFromId(2154), "ESRI Shapefile")

And the error message :

Traceback (most recent call last):
File "<string>", line 231, in processAlgorithm
TypeError: QgsVectorFileWriter.writeAsVectorFormat(): arguments did not
match any overloaded call:
overload 1: argument 1 has unexpected type 'str'
overload 2: argument 1 has unexpected type 'str'
overload 3: argument 1 has unexpected type 'str'

How can I export the final result as a shapefile ?


texts: Using OGR to export 3d polygons from postgis to DXF

texts:
Today I tried to export a 3-dimensional geometry out of postgis into a dxf file. After some tests it came out, that I was able to export

a 3d-linestring which represents the geometry or
a 2d-polygon which is basically the 3d geometry wih dropped z-values.

but not a closed 3d-polygon.

My environment

PostgreSQL 9.6, Postgis 2.3.2
GDAL 2.0.2

Let's create a single polygon and two tables for the ogr-export:

--cleanup
DROP TABLE IF EXISTS polygon;
DROP TABLE IF EXISTS test1;
DROP TABLE IF EXISTS test2;

--create single polygonz
SELECT
    ST_GeomFromText('POLYGON Z ((0.0 0.0 0.0  ,  10.0 0.0 10.0  ,  0.0 10.0 10.0  ,  0.0 0.0 0.0))') as geom
INTO
    polygon;

--create ogr table with polygon z
SELECT
    1 as gid,
    'Layer 1'::text AS layer,
    'BRUSH(fc:#ff0066,bc:#ff0066,id:"ogr-brush-0")'::text AS ogr_style,
    geom AS geom
INTO
    test1
FROM
    polygon;

--create ogr table with closed linestring z
SELECT
    1 as gid,
    'Layer 1'::text AS layer,
    'BRUSH(fc:#ff0066,bc:#ff0066,id:"ogr-brush-0")'::text AS ogr_style,
    ST_ExteriorRing((ST_Dump(geom)).geom) AS geom
INTO
    test2
FROM
    polygon;

And use the following script to export the geometry into DXF:

ogr2ogr -f "DXF" test1.dxf PG:"dbname=test3d host=127.0.0.1 port=5432 user=postgres password=postgres" -sql "SELECT * FROM test1"
ogr2ogr -f "DXF" test2.dxf PG:"dbname=test3d host=127.0.0.1 port=5432 user=postgres password=postgres" -sql "SELECT * FROM test2"

Result of test1: 2d hatch

Result of test 2: closed 3d polyline

Setting

 DXF_WRITE_HATCH=NO

leads in both cases to a 3d polyline, as shown in the image above.

I also tried to force ogr to create polygonz, using

-nlt POLYGONZ

which didn't improve the result, that should be a 3d polygon. Any idea how to achive this?


texts: SNAP vs. GDAL vs. Orfeo - image rescale, histogram cut

texts:
I work with a Sentinel-2 .jp2 image (red band, 10950 x 10950 pixels). My aim is to reach the same result what SNAP does with a Python script.  See the SNAP method and parameters:

So this is my reference (result with SNAP), I want to reach this result (QGIS grayscale representation, cumulative cut - 2/98%):

So I tried to replicate it with GDAL:

import numpy as np
from osgeo import gdal, gdal_array

input = "d:/bitbucket/cnn-lcm/T33TWM_A012703_20171127T100339_B04.jp2"
output = "d:/bitbucket/cnn-lcm/T33TWM_A012703_20171127T100339_B04_gdal.tif"

dataset = gdal.Open(input)
array = dataset.ReadAsArray()

percentile_025 = np.percentile(array, 2.5) # 349.0
percentile_975 = np.percentile(array, 97.5) # 3735.0

command = 'gdal_translate -scale ' + str(percentile_025) + ' ' + str(percentile_975)+ ' 0 255 -of GTiff -ot Byte' + ' ' + input + ' ' + output

os.system(command)

The GDAL result is not the same, its a bit brighter, the white areas are bigger. The values are not the same on the layers panel (QGIS grayscale representation, cumulative cut - 2/98%):

The Orfeo code:

import otbApplication

Convert = otbApplication.Registry.CreateApplication("Convert")
Convert.SetParameterString("in", "d:/bitbucket/cnn-lcm/T33TWM_A012703_20171127T100339_B04.jp2")
Convert.SetParameterString("out", "d:/bitbucket/cnn-lcm/T33TWM_A012703_20171127T100339_B04_hcut.tif")
Convert.SetParameterString("type","linear")
Convert.SetParameterString("hcp.high","2.5")
Convert.SetParameterString("hcp.low","2.5")
Convert.ExecuteAndWriteOutput()

Rescale = otbApplication.Registry.CreateApplication("Rescale")
Rescale.SetParameterString("in", "d:/bitbucket/cnn-lcm/T33TWM_A012703_20171127T100339_B04_hcut.tif")
Rescale.SetParameterString("out", "d:/bitbucket/cnn-lcm/T33TWM_A012703_20171127T100339_B04_orfeo.tif")
Rescale.SetParameterOutputImagePixelType("out", 1)
Rescale.SetParameterFloat("outmin", 0)
Rescale.SetParameterFloat("outmax", 255)
Rescale.ExecuteAndWriteOutput()

The Orfeo result is very similar to GDAL (only 1-2 value differences in pixels). And there are big, problematic strips in the middle (QGIS grayscale representation, cumulative cut - 2/98%):

So finally my questions:

Is it possible to eliminate the differences? Is it possible to reach exactly the result of SNAP? And how?

Download link to data:
http://sentinel-s2-l1c.s3.amazonaws.com/tiles/33/T/WM/2017/11/27/0/B04.jp2


texts: Does an FME fmw file include the processed data in it?

texts:
If I share an FMW file does it include the actual data or is it just the steps in the model and links to the data sources but not actual data.

Some models have links to Microsoft SQL Server databases, etc (without passwords)

I want to make sure private data is not shared when getting advise on updating a model or making it available under c3 for others to use.

Looking at the FMW in a text editor it doesn't seem to include any data - so it should be fine but just want to confirm this.


texts: Erasing overlapping polygons from same layer using ModelBuilder?

texts:
I currently have a shapefile with 5 catagories of overlapping polygons ranked 0-4, 1 is the most recent polygon, four is the oldest and 0 is undefined.

I am looking to use ModelBuilderr to seperate this file into each of the five catagories and then merge the entire layer diplaying only the highest ranked polygon (no overlapping older/lower ranked polygons drawn) as the file has over a 1000 polygons and the goal is to speed up drawing and use the final layer as asa reference to quickly visually identify the most current file for an area.

New to ModelBuilder.


texts: QGIS (Mis)Detecting Escape Characters in SLD files

texts:
Okay, I am just running into some trouble with the string literal functionality with QGIS's rule based styling.

For example I have an attribute called LITH that equals 'BAS\b'. The rule where LITH = 'BAS\b' does not work however as the escape character is not being ignored.

How would I go about ignoring the escape character in QGIS? If I was using c# I could use @'BAS\b', if I was using XML I would use the ![] tags.

I am loading these rules from a SLD file, but the symbology is incorrect in places due to the escape character not being ignored.


texts: Performing one-to-many join with address and crime data using ArcMap

texts:
Goal: I want to show the police calls coming from a parcel on ArcMap. So if you click on the parcel, it will show all the calls for the year. I want to join the address attribute from the parcel data to the police file's call address (which was created by concatenating street number and address to match the parcel layer)

Issue: many addresses have more than one call from different dates; one property even had 100 calls. The parcel layer has only one polygon for each address, so I can't show multiple records for the join.

How do I join so that all 100 calls would be joined to that parcel?


texts: PyQGIS standalone app remove/swichoff layers in the legend

texts:
I try to swich off remove layers from a legend in a qgis standalone application (QGIS 2.18) and found that the one way is to manipulate the QgsLayerTreeGroup (GIS.SE 142016).

# remove two unwanted layers from the legend set
mapset.remove(lyrBathy.id())
mapset.remove(lyrIso.id())

# create the legend
cleg = tk.createLegend( page, mapset, LEG_LEFT, LEG_TOP, 'Legend')

If I render the stuff at least the program crashes:

To be sure that the map is not influenced I make a copy of the current tree only the expected layers:

# ----------------------------------------------
# Remove layers not in list by copy on write
# ----------------------------------------------
def updateLegendLayerTree(srcTree, layers):
    print '---- BEFORE -----'
    print srcTree.dump()

    print  '--- CLEAN UP ----'
    for child in srcTree.children():
        if child.layer().id() in layers:
            result.addChildNode(child)
            print "..added", child.layer().name()
        else:
            print "..removed", child.layer().name()
    print '---- AFTER -----'
    print result.dump()
    return(result)

Create the legend with the filtered tree:

# ----------------------------------------------
# Create the legend
# ----------------------------------------------
def createLegend (page, layers,  x, y, title, \
                  titleFont = QFont('Arial',12,2), \
                  groupFont = QFont('Arial',11,1),
                  itemFont  = QFont('Arial',10,1) ):
    #print "Erstelle Legende: %s %d %d " % (title, x, y)
    cleg = QgsComposerLegend(page)
    cleg.setTitle(title)
    cleg.setAutoUpdateModel(False)
    cleg.rstyle(QgsComposerLegendStyle.Symbol).setMargin(QgsComposerLegendStyle.Top, 1)
    cleg.setStyleFont(QgsComposerLegendStyle.Title, titleFont)
    cleg.setStyleFont(QgsComposerLegendStyle.Subgroup, groupFont)
    cleg.setStyleFont(QgsComposerLegendStyle.SymbolLabel, itemFont)
    cleg.setX(x)
    cleg.setY(y)
    # Update the list
    srcRoot = QgsProject.instance().layerTreeRoot()
    newTree = updateLegendLayerTree(srcRoot, layers)
    # Set the new tree
    cleg.modelV2().setRootGroup(newTree)
    return cleg

I got this debug messages and the fatal error:

---- BEFORE -----
GROUP:  visible=2 expanded=1
  LAYER: Wassertiefe visible=2 expanded=1 id=Wassertiefe20170922114344525
  LAYER: Schutzgebiete visible=2 expanded=1 id=Schutzgebiete20170922114344676
  LAYER: AWZ visible=2 expanded=1 id=AWZ20170922114344696
  LAYER: Windpark visible=2 expanded=1 id=Windpark20170922114344504
  LAYER: Tiefenlinie visible=2 expanded=1 id=Tiefenlinie20170922114344550
  LAYER: Eisente visible=2 expanded=1 id=Eisente20170922114344818

--- CLEAN UP ----
..removed Wassertiefe
..added Schutzgebiete
..added AWZ
..added Windpark
..removed Tiefenlinie
..added Eisente

---- AFTER -----
GROUP: Legend visible=2 expanded=1
  LAYER: Schutzgebiete visible=2 expanded=1 id=Schutzgebiete20170922114344676
  LAYER: AWZ visible=2 expanded=1 id=AWZ20170922114344696
  LAYER: Windpark visible=2 expanded=1 id=Windpark20170922114344504
  LAYER: Eisente visible=2 expanded=1 id=Eisente20170922114344818

...write PDF: clc.2013-01-01.2017-01-01.cla.map-points.pdf
Speicherzugriffsfehler # (german... memory access error)

I assume the legend render tries to get things from dependenies which are not covered by the model, created in:

   def updateLegendLayerTree(srcTree, layers):
         .....
         result  = QgsLayerTreeGroup('Legend', True).

Any ideas what is the right statement? I tried to manipulate the visibility as well with negative results.


texts: Are mixed geometry types allowed in shapefiles?

texts: Does the shapefile support mixed geometry types (ie. point and polygon all in the same file)?


texts: Assigning raster values to line and polygon by QGIS?

texts:
I have a raster data(.tif), and would like to assign its values to
line and polygon data by location.
Here, I'd like to take min, max, mean and median.

There is "Join attributes by location", but it only supports
joining between vector data and vector data.

Any tools or ideas ?


texts: How to replace an ImageCollection stored as variable by an ImageCollection stored as Asset in a GEE script?

texts:
I'm trying to clip and export climate data (annual max) from one ImageCollection to an already existing table with (annual max) NDVI values and clip them with conflict points stored in FeatureCollection.
The script works for the NDVI, which I have as an asset in GEE. But since I'm creating annual max images in the script without saving them as an Asset, I have problems to use it.
//Calc annual max for TerraClimate for Temp, soil and prec and clip with conflicts 25
var annual_max_temp =  ee.ImageCollection(ee.List.sequence(1981, 2015).map(function (year){
  var date_start = ee.Date.fromYMD(year, 1, 1);
  var date_end = date_start.advance(1, "year");

  return(ee.ImageCollection('IDAHO_EPSCOR/TERRACLIMATE')
        .select(['tmmx','soil', 'pr'])
        .filterDate(date_start, date_end)
          .max()
          .set({year: year, 'system:time_start':date_start.millis()})
          .set({year: year, 'system:time_end':date_end.millis()}));
}));

print(annual_max_temp,'at');
Map.addLayer(annual_max_temp,{color: '0036FF', pointRadius: 6}, 'at');

//Clipping the Climate IC with the Conflict FC and adding them to the NDVI_Slopes_Table
var GIMMSConfl25 = NDVI_Slopes
  .map(function(f) {
    var maxNDVI = annual_max_temp //when I use this it works: var maxNDVI = ee.ImageCollection('users/welp1992/GIMMS_Max_Annual_NDVI')
                      .filterDate('1982-01-01', '2015-12-31')
                      .map(function(img) {return img.clip(Conf25)})
                      (['tmmx','soil', 'pr'])
      .map(function(img) {
        var d = img.reduceRegion({
                    reducer: ee.Reducer.first(),
          geometry: f.geometry(),
          scale: 8000
        });
        return f.set(d).set('GIMMSdate', img.date().format('yyyyMMdd'));
      });

    return maxNDVI; //global variable, use later in script

  }).flatten();

// Get Info
//print(NDVI_Slopes.limit(650));
print('Annual_Max_Clipped',Annual_Max_Clipped.limit(650));
print('Annual_Max_Clipped',Annual_Max_Clipped.size())

Export.table.toDrive({
  collection: GIMMSConfl25,
  description: 'annual_max_temp_table',
  fileFormat: 'CSV'
});

I get this error.
annual_max_temp.filterDate(...).map(...) is not a function
How do I replace the ImageCollection in the code?
I think it's a simple problem in the code. But if needed I can share the code and the assets.


texts: Creating multiple polygon grid in ArcGIS for Desktop?

texts:
In ArcGIS 10.3 I'm able to create a polygon grid with "grid index feature" but the problem is that I would like create a very big polygon grid (the size like an entire continent) and with the cell size equal to a square with the side of about 250 metres. It's obvious that a normal computer is unable to manage this layer due to the big cell amount, we talk about billion cells!

My idea is to create "N" small polygon grid (always with a cell size about 250 metres) to cover my area and I don't want join attribute or other kind of operation. The problem is that when I create the first polygon grid, the second one needs to start (or to end) by overlap the last cell-line or cell-row of the precede polygon grid and so on.

I don't want an automatic process, I can  manually create each polygon grid in ArcGIS but which are the steps to create this N polygon grid with this "special case"?

An image explain what I want, in this image there is only the North overlap but I need to consider also the case when the overlap in South or East or West):


texts: problems setting up postgis and importing data

texts:
Using Postgres 9.2.3.1 on a new Mac.  I installed Postgres from the dmg file that I downloaded from the Postgres site.  I then installed PostGIS using stackbuilder. This all worked flawlessly so far as I can tell.

I have tweaked the pg_hba.conf file to trust connections from local and localhost and created a database "db".  Then I tried::

Mac-mini:pg postgres$ bin/createlang plpgsql db
createlang: language "plpgsql" is already installed in database "db"

so it would appear that postgres knows all about postgis.

When I tried:

Mac-mini:pg postgres$ bin/shp2pgsql -c ~rful011/tmp/weed_boundaries weed_plot_boundaries | psql db
Shapefile type: Arc
Postgis type: MULTILINESTRING[2]
SET
SET
BEGIN
NOTICE:  CREATE TABLE will create implicit sequence "weed_plot_boundaries_gid_seq" for serial column "weed_plot_boundaries.gid"
CREATE TABLE
NOTICE:  ALTER TABLE / ADD PRIMARY KEY will create implicit index "weed_plot_boundaries_pkey" for table "weed_plot_boundaries"
ALTER TABLE
ERROR:  function addgeometrycolumn(unknown, unknown, unknown, unknown, unknown, integer) does not exist
LINE 1: SELECT AddGeometryColumn('','weed_plot_boundaries','geom','0...
               ^

I must have missed something but what?

Thanks, Russell


texts: How can I get Sokkia total station .RAW file data into ArcGIS

texts: I have a .raw file from a Sokkia total station.  I want to view the data in ArcGIS and create contours from this data.  Can anyone help?  The survey analyst extension is longer supported in ArcGIS 10.1 so I'm not sure where to begin...


texts: BTS 2014 Highway Performance Monitoring System over 2gb dbf limit

texts:
http://www.rita.dot.gov/bts/sites/rita.dot.gov.bts/files/publications/national_transportation_atlas_database/2014/polyline

The dbf in question is part of the hpms.shp

The error I receive when trying to import it into postgis (either via shp2pgsql or ogr2ogr) is an fread(782) error. Basically the size of the unzipped dbf file is a bit over 2gb, which is known limitation.

Anyone know how to get around the limitation?

I've contacted the BTS to let them know about this issue...


texts:
I have some GEOJSON-Layers in my QGIS-project. The problem is, that when I zoom in or out, it takes very long to built up the window again. That's because it always loads it from the server. But I don't need that. It's enough, just to load it when I start the project and the handle it like a shapefile. Of course, I can save the WMS-Layer as a shapefile. But doing this everyday with every WMS-Layer would me to much work.

Is there an easy solution?


texts: Change cursor when hovering over features in WMS and ArcGIS tilelayers

texts: I'm using Mapbox with Leaflet to retrieve WMS and ArcGIS tiles.  I have added the ability to retrieve attributes from the WMS servers and ArcGIS servers with WMS capabilities turned on.  I would like to have the mouse cursor change to a pointer when hovering over a feature to show that it's clickable.  How do I do this?


texts: Removing part of polygon from another using ArcGIS for Desktop?

texts:
Now I have shp of a county, which includes lake part and land part.

I only have shp of the lake part but I need the land part.

How can i get the land part?

I think there should be some tools doing this, but I don't know its name. Clip is totally doing the reverse work.


texts:
Is it possible to select a polygon from a predefined Lat & Lon programmatically?

eg location point: x='1234' y='5678' intersects polygon '1234', therefore SELECT POLYGON '1234'

So far I have managed to get the lat lon of the drop point:

$( "#cow" ).draggable({ revert: true});
        $( "#map" ).droppable({
        drop: function( event, ui ) {
    var point = map.getLonLatFromPixel( map.events.getMousePosition(event, ui   ))
        OpenLayers.Util.getElement("pos").innerHTML = point

        }
        });

Now I just neeed to select the polygon that x and y intersect. Ideas?


texts: Delete photos with no location in EXIF information

texts:
I can extract the locations of my photos in QGIS using Photo2Shape or Geotag and import photos, but I want to be able to delete the photos that do not have location information with them (they are backed up elsewhere).

At the moment I have to look in the attribute table which photos get imported with a latitude and longitude value of 0. Then I have to find these files manually and delete them. This is very time consuming and it would be great if it was possible to do this automatically somehow.

I guess a another possible workaround would be for the photos that have locations to be copied to a sub folder.

I just wondered if anyone knew of a tool (GIS or otherwise) that could achieve this before I set about writing a script to do it for me.


texts: Format Coordinates Checker with QGIS or other Appilcations

texts:
I have thousand list coordinates from excel format (CSV file), its projection is  Geographic WGS 84 (EPSG 4326) most of all. All of the coordinates are inserted manually from our surveyor teams. But we know that sometimes we make mistakes by putting wrong format coordinates, for example the correct coordinates is "6.504302, 112.786543" but they insert it as "6504302, -112786543" because of human errors. It's quite tiring to find the mistakes among thousands list excel cell and sheet, it also takes times too.

Are there apps, tools, or methods to find and replace the mistakes with correct values?


texts:
All I need is the ability to create random points within a radius of another point.

so given self.fence.latitude and self.fence.longitude and self.fence.radius I should be able to generate random points with the following code (i'm writing in python/django and using numpy but all these calcs can be done using the standard math library)

for whatever reason my points are showing up in a line. I need random points.

# angle for this random point
angle = numpy.float128(2 * numpy.pi * numpy.random.random_sample())
# radius from fence center
radius = numpy.float128(self.fence.radius) * numpy.random.random_sample()

# diff from ceneter in radians
latitudinal_diff_rad = radius * numpy.sin(angle) / numpy.float128(settings.EARTH_RADIUS)
longitudinal_diff_rad = radius * numpy.cos(angle) / numpy.float128(settings.EARTH_RADIUS)

# radians to degrees
latitudinal_diff = latitudinal_diff_rad * 180 / numpy.pi
longitudinal_diff = longitudinal_diff_rad * 180 / numpy.pi

try :
     # create the item
    item = Item(
        latitude = numpy.float128(self.fence.latitude) + latitudinal_diff,
        longitude = numpy.float128(self.fence.longitude) + latitudinal_diff,
    )
    item.save()

except IntegrityError, e :
    continue
    print e

What am i doing wrong?

Edit

# this was the math I was using before I simplified it and switched to numpy lib
distance_from_point = float(self.fence.radius) * random.random()
unit_length = distance_from_point / settings.EARTH_RADIUS
longitudinal_diff_rad = ((unit_length) * (2 * random.random() - 1))
latitudinal_diff_rad = math.acos(math.cos(unit_length) / math.cos(longitudinal_diff_rad))


texts:
How can I render a vector layer into .tif file? I tried with below codes without luck.

PS: I am doing this on ubuntu server and running it via PyQGIS.

Python scripts

#!/usr/bin/python
from qgis.core import *
from qgis.utils import *
from PyQt4.QtCore import *

def main():
  qgis_prefix="/usr"
  QgsApplication.setPrefixPath(qgis_prefix, True)
  QgsApplication.initQgis()
  app = QgsApplication([], False)

  uri = "elevp.csv?delimiter=%s&xField=%s&yField=%s&elevField=%s" % (";","x","y","elev")
  vlayer = QgsVectorLayer(uri, "test", "delimitedtext")

  # create image
  img = QImage(QSize(800, 600), QImage.Format_ARGB32_Premultiplied)

  # set image's background color
  color = QColor(255, 255, 255)
  img.fill(color.rgb())

  # create painter
  p = QPainter()
  p.begin(img)
  p.setRenderHint(QPainter.Antialiasing)

  render = QgsMapRenderer()

  # set layer set
  lst = [vlayer.getLayerID()]
  render.setLayerSet(lst)

  # set extent
  rect = QgsRect(render.fullExtent())
  rect.scale(1.1)
  render.setExtent(rect)

  # set output size
  render.setOutputSize(img.size(), img.logicalDpiX())

  # do the rendering
  render.render(p)

  p.end()

  # save image
  img.save("render.tif","tif")

main()

I tried to capture the debug error with pdb and it show as below: -

-> /path/to/pyqgis/test.py(40)main()
-> rect = QgsRect(render.fullExtent())

Can I also ask if pdb is the correct debugging tool to use?


texts:
I have a carto-map that consists of four layers.
The users of my map should have the option to choose the map projection.
For one layer this works fine:

$('#selectProjection').on('click', function() {
            var srid = $(this).val();
            firstLayer.setSQL("SELECT cartodb_id, st_transform(the_geom, " + srid + ") as the_geom_webmercator FROM firsttable");
 });

But i have to project all four layers. If a add the same query just for another layer, there is no layer projected at all.
This code doesn't work:

$('#selectProjection').on('click', function() {
            var srid = $(this).val();
            firstLayer.setSQL("SELECT cartodb_id, st_transform(the_geom, " + srid + ") as the_geom_webmercator FROM firsttable");
            secondLayer.setSQL("SELECT cartodb_id, st_transform(the_geom, " + srid + ") as the_geom_webmercator FROM secondtable");
});

Does anoyone have an idea?

EDIT: a js fiddle with the whole code: https://jsfiddle.net/nmrs11kq/11/


texts: Array written by GDAL with Python script has numbers read as "nan" by QGIS?

texts:
I would like to write an array to a geotiff with Python script using the GDAL module.  When I try to do something simple like copying an existing GeoTiff via reading it to an array, then, writing to a new array, QGIS reads the integer value 1 as "nan".  I have verified that the actual values that QGIS presents as "nan" are the integer value 1 in the GeoTiff.  Is there some additional information that I need to write to the GeoTiff file so that QGIS will read my integer values of 1 as 1?  My code is copied below:

import os
import gdal
os.chdir(r'C:\Users\pipi\Documents\Rogaine\Tarlo\gpx')

in_ds = gdal.Open('0011aa0011aa.tif')
in_band = in_ds.GetRasterBand(1)
noData = in_band.GetNoDataValue()
if noData:
    pass
else:
    noData = 0

cols = in_band.XSize
rows = in_band.YSize

# Creates an output file with the same dimensions and data type.
out_ds = gdal.GetDriverByName('GTiff').Create(
'master.tif', cols, rows, 1, gdal.GDT_Byte)
out_ds.SetProjection(in_ds.GetProjection())
out_ds.SetGeoTransform(in_ds.GetGeoTransform())
out_band = out_ds.GetRasterBand(1)

data = in_band.ReadAsArray()
out_band.WriteArray(data)
out_band.SetNoDataValue(noData)

del out_ds


texts: how to upload big files to GEE as assets with the API?

texts:
I try to upload big files to Google Earth Engine (10 000 points on the surface of a medium country) as assets. If I go through the code editor and Assets -> New -> .csv file It work like a charm.
I am using the Python API and I would like to do the same without going to the code editor.
What I did is the following :
import ee
import geemap
import pandas as pd
import geopandas as gpd

# create a tmp gdf
df = pd.read_csv(load_df['pathname'], sep=None, engine='python')
gdf = gpd.GeoDataFrame(
    df,
    crs='EPSG:4326',
    geometry = gpd.points_from_xy(
        df[load_df['lng_column']],
        df[load_df['lat_column']]
    )
)

# convert it into geo-json
json_df = json.loads(gdf.to_json())

# create a gee object with geemap
ee_object = geemap.geojson_to_ee(json_df)

# upload this object to earthengine
asset = os.path.join(folder, asset_name)

#create and launch the task
task_config = {
    'collection': ee_object,
    'description':asset_name,
    'assetId': asset
}
task = ee.batch.Export.table.toAsset(**task_config)
task.start()

Problem being the created asset is too big to be exported. Is there another way to do it ?


texts: How do I physically order overlapping ArcGIS features?

texts: I have a feature class with multiple overlapping polygons. It's a highways fc and contains underpasses, overpasses and at grade roads & intersections. The problem I have encountered is that the items are not returned in the order I desire. Is there a way to correct this? From what I've read, it is related to the object ids so the only solution I can see would be to write some python code to loop through each feature based on a 'level' attribute and write those with the smallest value (say -1) to a new fc first, followed by levels 0 and 1. I'm using ArcGIS 10.0


texts: Spectral mixture analysis from GeoTIFF in ENVI

texts:
I have a single GeoTIFF created from some NGB aerial photography that has been stitched in Agisoft Photoscan to create a single image of an area of heath.

I would like to perform spectral mixture analysis on the GeoTIFF but am lost as to where to start. I have sample photographs for end member selection. Can anyone point me in the right direction please?


texts:
I am using QGIS 3.18.1-Zürich and I am trying to create a Python action. Although I have installed the required packages, I get an error as shown below.
Traceback (most recent call last):
  File "", line 2, in
  File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 799, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ModuleNotFoundError: No module named 'pandas'

Python path:
'C:/OSGEO4~1/apps/qgis/./python',
'C:/Users/User/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python',
'C:/Users/User/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python/plugins',
'C:/OSGEO4~1/apps/qgis/./python/plugins',
'C:\\OSGEO4~1\\apps\\Python37',
'C:\\OSGEO4~1\\apps\\Python37\\Scripts',
'C:\\OSGeo4W64\\bin\\python37.zip',
'C:\\OSGEO4~1\\apps\\Python37\\DLLs',
'C:\\OSGEO4~1\\apps\\Python37\\lib',
'C:\\OSGeo4W64\\bin',
'C:\\OSGEO4~1\\apps\\Python37\\lib\\site-packages',
'C:/Users/User/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python'

What can I do to solve this problem?
Edit:
I have already installed pandas in Python 3.9 but QGIS is using 3.7.
WARNING: Value for scheme.headers does not match. Please report this to <https://github.com/pypa/pip/issues/9617>
distutils: c:\users\user\appdata\local\programs\python\python39\Include\UNKNOWN
sysconfig: c:\users\user\appdata\local\programs\python\python39\Include
WARNING: Additional context:
user = False
home = None
root = None
prefix = None
Requirement already satisfied: pandas in c:\users\user\appdata\local\programs\python\python39\lib\site-packages (1.2.4)
Requirement already satisfied: python-dateutil>=2.7.3 in c:\users\user\appdata\local\programs\python\python39\lib\site-packages (from pandas) (2.8.1)
Requirement already satisfied: pytz>=2017.3 in c:\users\user\appdata\local\programs\python\python39\lib\site-packages (from pandas) (2021.1)
Requirement already satisfied: numpy>=1.16.5 in c:\users\user\appdata\local\programs\python\python39\lib\site-packages (from pandas) (1.20.2)
Requirement already satisfied: six>=1.5 in c:\users\user\appdata\local\programs\python\python39\lib\site-packages (from python-dateutil>=2.7.3->pandas) (1.15.0)
WARNING: Value for scheme.headers does not match. Please report this to <https://github.com/pypa/pip/issues/9617>
distutils: c:\users\user\appdata\local\programs\python\python39\Include\UNKNOWN
sysconfig: c:\users\user\appdata\local\programs\python\python39\Include
WARNING: Additional context:
user = False
home = None
root = None
prefix = None


texts: Trimble Net R9 Base Station Receiver File Structure

texts:
We just installed our new GPS Base Station Net R9 receiver that replaced our old 5700 receiver. The old 5700 used Trimble's GPS Base software, which the new Net R9 does not. My first problem that I can't seem to figure out is the file structure. With the old GPS Base software we were able to set the data logging system to write the files in folders that went Year -> Month -> Day -> hour. It was a nice clean structure. The new Net R9 logs all files in one folder where we have it set for every hour. This is going to be a massive folder after a while. Does anyone have any experience with the Net R9 file structure and how to get it to write like the old GPS Base software? Below are some pictures. Thanks in advance!!
Trimble 5700 GPS Base File Structure
Net R9 File Structure


texts: Using OpenLayers and JSTS to create buffer geometry

texts:
So far I have found similar questions like this one but none have quite expressed or fixed my exact issue.

I'm using Openlayers v4.6.5 and JSTS 1.4.0

The problem I'm having is with the buffer() method of the JSTS library. My conde snippets are as follows:

var feats = vectorSource.getFeatures();

// Takes in multiple features and combines them into one
var geom = joinGeoms(feats);

var parser = new jsts.io.OL3Parser();

var jstsGeom = parser.read(geom);
var meters = toMeters(500, 'ft'); // Result is 152.4 meters
var buffered = jstsGeom.buffer(meters);

// This is where the issue happens
var finalGeom = parser.write(buffered); <---

The geometry that is generated from the buffer() method and assigned to finalGeom is then used to select other features that intersect. I have captured this geometry and have shown it below along with measurements done with the my measure tool (which I have tested separately and know is accurate).

The buffer distance should be 500 feet or 152.4 meters but comes back as a different distance. I have to tried seeing if this is from a different projection or unit conversion but I cant find anything that would take my inputs and result in the numbers I am getting.

My projection is EPSG:3857, which I know isn't the best for these kinds of measurements but I have been unable to successfully convert where its accurate.

Any thoughts?


texts: Why do most GIS packages need a numeric id?

texts:
This is a simple yet possibly controversial question: why do most (if not all) GIS packages require that a determined layer have a unique not nullable numeric identifier?

Why there is the need for such a surrogate key instead of a natural one?

Examples:

ArcGIS enforces OBJECTID (or a GlobalID)
QGIS does not load layers when they don't have a numeric id.


texts: Convert .tiff data to bitmap just when rendering and at runtime by visualC++ code

texts:
I need to convert a .tiff image to a bitmap one without saving the bitmap. I know that I should use GDAL_translate but because I'm new, this function's interface is a little complicated for me. Could you please tell me how should I call the function in this special case?
And if you have any better suggestion than using GDAL_translate, please allow me be informed?

Note that I don't want to save the bitmap. It's just for displaying purposes in mfc and my superior has told me not to save the bitmap just convert data to bitmap at run-time.

well @user1240592,@Mapperz,@Alex Marcov,@Jamo,@Tonny,@johanvdw and @soheil I've seen the question you've answered here. my problem is I don't need a bitmap file to be created on the hard disk.
suppose that I have this MFC application  already I have written the code to print header information of a geotiff in the rich edit box now I need something to convert a 15931x16997 geotiff to bitmap at runtime and show it in a picture control ( I need a high performance way) and be able to do various image processing operations on this bitmap data at run-time (using opencv) and immediately show the result on the picture control and also be able to do these tasks when the image is viewed:
1-zoom in
2-zoom out
3-pan
4-fixed zoom in
5-fixed zoom out
6-go back to previous extent
7-go next extent
and then if needed I can convert the bitmap data created by image processing operations to a geotiff image file and convert it in the hard disk.
I don't have problems with image processing parts,opening and saving geotiff images, extracting raster data and write raster data.
my problem is converting data to bitmap at runtime and use a high performance way (like openGL) to display the image with high performance and do various navigation tasks on it.
Do you anyway or any library that have implemented these tasks Notice just need bitmap data on the buffer not saved on the hard disk?


texts: Open source GIS software to use for utility management?

texts:
I am looking for open source GIS software which can be used for utility management (like water, gas, electricity, road, etc.).

I have done this in ArcMap 10.1 but am now looking for open source GIS software.

Can anyone suggest to me any open source GIS software for the same?


texts: How to make contingency table for two satellite images in R language

texts:
I am working on forest change detection in an area. I have taken two satellite images 2001 and 2016 to depict the forest change. I have taken four classes while classifying the images i.e., forest, non-forest, snow and water bodies.

I want to know, how we can make contingency table of both the  two satellite images. The contingency table must show the change of pixel like how many pixels of class forest were present in 2001 and how many of those pixels transferred in 2016 in another class.


texts: Count overlapping features using Geopandas

texts:
Is there a way to easily count overlapping polygons using Geopandas, the same way as the ArcGIS Pro Count Overlapping Features works?

So far my approach was to do union overlay and then dissolve with aggfunc='count' but for some reason the results I get are not correct.
I have 3 overlapping polygons in a single geodataframe:

Then I do the overlay:
union = gpd.overlay(demo_pg, demo_pg, how='union')

As a result I get only 9 polygons, although I should get 10 (this is what union in QGIS or ArcGIS would return):

Is there anything wrong with my approach? What is the best way to count overlapping polygons in a single geodataframe?
EDIT: Full code is below. It returns 9 polygons. Based on my understanding on union/intersect operations, it should result in 10 polygons. The intersection of 3 polygons is counted only twice, not three times... The union operation in QGIS for the same set of polygons results in 10 polygons.
import pandas as pd
import matplotlib as plt
import geopandas as gpd
from shapely import wkt
data = {'name':  ['polygon A', 'polygon B', 'polygon C'],
        'id': [1, 2, 3],
         'geom': ['MULTIPOLYGON (((36.00000 11.00000, 36.00000 12.00000, 37.00000 12.00000, 37.00000 11.00000, 36.00000 11.00000)))', 'MULTIPOLYGON (((36.50000 11.50000, 37.50000 11.50000, 37.50000 11.00000, 36.50000 11.00000, 36.50000 11.50000)))', 'MULTIPOLYGON (((36.61799 10.80580, 36.61570 11.19321, 36.86327 11.29637, 37.34925 10.91813, 37.00540 10.71182, 36.61799 10.80580)))']
        }

df = pd.DataFrame (data, columns = ['name','id','geom'])
df['geom'] = df['geom'].apply(wkt.loads)
gdf = gpd.GeoDataFrame(df, geometry='geom')

demo_pg.plot(alpha=0.5, column='id')
union = gpd.overlay(demo_pg, demo_pg, how='union')
len(union)


texts:
I'm not sure of the correct terminology, so I'll describe what I'm looking for.

I've used it before, but can't remember where, but what is the tool/plugin to use when adding features to say a point feature and it retains the attribute data from the last point added, so that the new attribute data is pretty much already populated?

I am using QGIS 3.0.3.


texts: Which tool is recommended for a PostGIS design?

texts: I need to create a design for a spatial database. I've found this list http://wiki.postgresql.org/wiki/GUI_Database_Design_Tools, but it seems to me that this software is orientated towards a regular (not spatial) DB design.


texts:
How to extract streetnames with osmosis (Java)?
I have tried to process all Nodes. I put them in a Map<String, List> meaning Map<CityName, ListOfStreetnamesOfThatCity>:

@Override
public void process(EntityContainer entityContainer) {
    if (entityContainer instanceof NodeContainer) {
        final NodeContainer entityContainer1 = (NodeContainer) entityContainer;
        if(containsTagWithKeys(entityContainer1.getEntity().getTags(), "addr:city", "addr:street")) {
            final Map<String, String> tags = entityContainer1.getEntity().getTags()
                                                  .stream()
                                                  .collect(Collectors.toMap(Tag::getKey, Tag::getValue));
            if( cityStreetnames.get(tags.get("addr:city")) == null) {
                    cityStreetnames.put(tags.get("addr:city"), new ArrayList<>());
            }
            String streetname = tags.get("addr:street");
            cityStreetnames.get(tags.get("addr:city")).add(streetname);

        }
    }
}

This gives me a subset of all streets per city only. (Not sufficient)

I tried to process all Way nodes

if (entityContainer instanceof WayContainer) {
    Way myWay = ((WayContainer) entityContainer).getEntity();

    final Map<String, String> tags = myWay.getTags()
                                             .stream()
                                             .collect(Collectors.toMap(Tag::getKey, Tag::getValue));
    if(containsTagWithKeys(myWay.getTags(), "name", "highway") && tags.get("name").length() > 2) {
        names.add(tags.get("name"));
    }
}

This gives me a lot of streets but I don't know in which city these ways are located.

I'm really new to OSM data. Is there a way to map a way node to a city? Is there a better way to get all streetnames of a city?

I have used .osm.pbf for Germany from geofabrik.de.


texts: Overcoming error : classified: Layer error: Image.reduceConnectedComponents: Label band must be integer valued, found type 'Float'

texts:
with following code , I am getting error message
https://code.earthengine.google.com/15a22656f680c38ff0bf10afd2747c5f
classified: Layer error: Image.reduceConnectedComponents: Label band must be integer valued, found type 'Float'.
I adapted code from following link which is not clear is it OBIA , what ??
https://blog.csdn.net/qq_21567935/article/details/85564946


texts: QGIS vector field marker not respecting map rotation

texts:
I’m working on a map the shows point displacement between two dates. These are survey benchmarks that are surveyed periodically to observe slope movement. Symbolizing this is QGIS is straight forward with the vector field marker by supplying two fields with the displacement values (delta) in X and Y direction. The map and all the layers are in EPSG:3057 that uses meters. The displacement values are in meters as well. My problem is that I need to rotate the map 130° CCW for publication but the arrows do not follow the map rotation (like a north arrow for example). The arrows keep point in the same “paper direction” or “screen direction”. Am I missing something in the symbology setting for the vector field marker?
The layer is a shapefile and this is QGIS 3.14.15 running in Windows 10.


texts: How to print a vector layer over Google map? "Physical Printing"

texts:
I use Geoserver, Openlayers, Extjs and GeoExt to build an app.

I am trying to get vector layer on Google map printed with pop up a Print dialog, I tried :

var targetElement = Ext.getCmp('mapPanel');
                    var myWindow = window.open('', '', 'width=600,height=500');
                    myWindow.document.write('');
                    myWindow.document.write('<title>' + 'Title' + '</title>');
                    myWindow.document.write('');
                    myWindow.document.write(targetElement.body.dom.innerHTML);
                    myWindow.document.write('');
                    myWindow.print();

it works fine except that Vector layer were not printed. I got Google map only. How I can do that ?


texts: Use Google Earth Engine to read MODIS NDVI into python

texts:
I'm trying to adapt this example for my needs: https://mygeoblog.com/2017/10/06/from-gee-to-numpy-to-geotiff/

The difference is that I'm using the MODIS NDVI image collection.

# define the image
img = ee.ImageCollection('MODIS/MCD43A4_006_NDVI').filterDate(date_start, date_end)

median_data = img.median()

# get the lat lon and add the ndvi
latlon = ee.Image.pixelLonLat().addBands(median_data)

# apply reducer to list
latlon = latlon.reduceRegion(
  reducer=ee.Reducer.toList(),
  geometry=area,
  scale=10000);

# get data into three different arrays
data = np.array((ee.Array(latlon.get("NDVI")).getInfo()))
lats = np.array((ee.Array(latlon.get("latitude")).getInfo()))
lons = np.array((ee.Array(latlon.get("longitude")).getInfo()))

I thought I understood the example but my data, lats and lons arrays are coming out different sizes which is obviously incorrect.

In [2]: data.shape
Out[2]: (21005,)

In [3]: lats.shape
Out[3]: (21051,)

In [4]: lons.shape
Out[4]: (21051,)

Any ideas what I'm doing wrong?


texts: Finding clusters of one set of points around another set of points from different layer?

texts:
I want to see if there is clustering of a certain type of buildings (x) around another type of buildings (y).

The two point files are in different layers.

I cannot figure out which tool I would use to do this.


texts:
I'm embarking on a new web map project that seeks to display simple geometries (lines, points, polygons) and rasters/basemaps. Ideally, the map will also allow authorized users to add/remove/edit geometries and their attributes.

What are the available options for the storage of data (e.g. SQL Server Spatial)?

What are the available options for the presentation of data (e.g. ArcGIS Server)?

I'm new to web mapping and am attempting to build a solid understanding of the available options and their pros/cons.


texts: Custom forms with Python logic: Run error

texts:
I'm a beginner on Python, PyQt5 and QGIS, I'm develop an form to capture some information and when a try to define a custom Python to make some validations it doesn't run.
I really need fix this error somebody can help me to identify how can I solve this problem ?
For picture sample it is code

This is message error with an simple code, but is similar with any code.


texts: Enhanced vegetation index using ENVI 4.7 bandmath

texts:
I computed Enhanced Vegetation Index using band math in ENVI 4.7 software. The expression used was:

2.5*((float(b1)- float(b2))/(float(b1)+ (6*float(b2)) +(7.5* float(b3))))
BYTSCL(float(b1), min = -1.0, max = 1.0)

The output EVI file has values varying from 0 to 255 rather than -1 to +1. Can anyone suggest an edit to the expression or a method.


texts: Optimising a very large point in polygon query

texts:
I have a national dataset of address points (37 million) and a polygon dataset of flood outlines (2 million) of type MultiPolygonZ, some of the polygons are very complex, max ST_NPoints is around 200,000.  I am trying to identify using PostGIS (2.18) which address points are in a flood polygon and write these to a new table with address id and flood risk details.  I have tried from an address perspective (ST_Within) but then swapped this starting from the flood area perspective (ST_Contains), rationale being that there are large areas with no flood risk at all.  Both datasets have been reprojected to 4326 and both tables have a spatial index.  My query below has been running for 3 days now and shows no signs of finishing any time soon!

select a.id, f.risk_factor_1, f.risk_factor_2, f.risk_factor_3
into gb.addresses_with_flood_risk
from gb.flood_risk_areas f, gb.addresses a
where ST_Contains(f.the_geom, a.the_geom);

Is there a more optimal way to run this?  Also, for long running queries of this type what is the best way of monitoring progress other than looking at resource utilisation and pg_stat_activity?

My original query finished OK albeit for 3 days and I got sidetracked with other work so I never got to dedicate the time to trying out the solution.  However I have just re-visited this and working through the recommendations, so far so good.  I have used the following:

Created a 50km grid over the UK using the ST_FishNet solution suggested here
Set the SRID of the generated grid to British National Grid and built a spatial index on it
Clipped my flood data (MultiPolygon) using ST_Intersection and ST_Intersects (only gotcha here was I had to use ST_Force_2D on the geom as shape2pgsql added a Z index
Clipped my point data using the same grid
Created indexes on the row and col and spatial index on each of the tables

I am ready to run my script now, will iterate over the rows and columns populating results into a new table until I have covered the whole country.  But just checked my flood data and some of the very largest polygons appear to have been lost in translation!  This is my query:

SELECT g.row, g.col, f.gid, f.objectid, f.prob_4band, ST_Intersection(ST_Force_2D(f.geom), g.geom) AS geom
INTO rofrse.tmp_flood_risk_grid
FROM rofrse.raw_flood_risk f, rofrse.gb_grid g
WHERE (ST_Intersects(ST_Force_2D(f.geom), g.geom));

My original data looks like this:

However post clipping it looks like this:

This is an example of a "missing" polygon:


texts: Python script for Data Driven Pages to run a series of DDP at once?

texts:
I have a several different polygon shapefiles that I am using to create individual maps for each polygon within the shapefile.

The layout is the same for each map, the only thing that changes is the polygon and scale bar. The size of the polygons vary greatly (i.e. < 1 acre to > 5,000 acres). Because of the size variability, I've had scale bar issues with running all polygons in DDP with the the division value either being too large or small.

Currently, I am running a DDP on each shapefile with correct scale bar division value, one shapefile at a time. I have 12 shapefiles in this project.  Is there a way to automate this using python to run DDP in a series on the different shapefiles?

Any advice would be appreciated, Thanks.


texts:
I will try to be as detailed as I can with the purpose and walkthrough of the situation.

I have 2 Feature Classes (fc1 and fc2).

fc1 is a polygon shapefile with no overlapping polygons and an id for each polygon on field code (Int)

fc2 is a polygon shapefile with overlapping polygons and an id for each polygon on field pcode (String)

I want to find for each polygon in fc1 the pcode that covers it the most. If more than one pcodes covers it the most, then I pick the pcode that is sorted FIRST.

In order to find a way through this, I wrote a code with the following steps and it produces the result that I want.

from collections import OrderedDict
F = OrderedDict([])
arcpy.Union_analysis([[fc1],[fc2]],'Test1','ALL','','GAPS')
arcpy.MakeFeatureLayer_management('Test1','Test2','"code" <> -1')
arcpy.Dissolve_management('Test2','Test3',['code','pcode'],'','SINGLE_PART','DISSOLVE_LINES')
arcpy.AddField_management('Test3','area','DOUBLE')
with arcpy.da.UpdateCursor('Test3',['SHAPE@AREA','area']) as cursor:
    for row in cursor:
        row[1] = row[0]
        cursor.updateRow(row)
arcpy.Sort_management('Test3','Test4','code ASCENDING;area DESCENDING;pcode ASCENDING')
arcpy.Statistics_analysis('Test4','Test5',[['pcode','FIRST']],'code')
with arcpy.da.SearchCursor('Test5',['FIRST_pcode','code']) as cursor:
    for row in cursor:
        F.update({row[1]:row[0]})
with arcpy.da.UpdateCursor(fc1,['code','pcode']) as cursor:
    for row in cursor:
        row[1] = F[row[0]]

The previous version of my question was posted because I was looking for a faster way because the above code takes about 6 minutes to finish.

From my search here and in other websites, I didn't find anything that can alter Dissolve with sth else unless someone has sth to suggest (it is the most time consuming process here).

For the rest, I didn't manage to find anything, so I hope you can help me out with this.

Supplemental for @radouxju comment:

I used the Tabulated Intersection tool but for some reason it miscalculates the area and the results are not correct. The polygon with code "4" is completely covered by the polygons with pcode 45194229099040 and 45194229099039 but the tool understands that it is partly covered. See the image below:

As you can see there is a difference on the 3rd decimal and because of this the Summary Statistics selects the incorrect pcode (in this case 45194229099040 instead of 45194229099039).


texts: Is there any water quality toollbox which can be used in water quality analysis

texts: I just start my job in water quality analysis,I wonder that do we have any toolbox to do water quality analysis in GIS


texts: Geoserver - Oracle Plugin not Working....?

texts:
I am using Geoserver 2.2.5 & added 'geoserver-2.2.5-oracle-plugin.zip' at '..\geoserver 2.2.5\webapps\geoserver\WEB-INF\lib' to connect with oracle data.

after restarting Geoserver, No datastore like 'Oracle/Oracle NG' is found.

The two jar files(in 'lib' directory) are -
1) ojdbc.jar ;
2) gt_jdbc_oracle_8.7.jar.

kindly,note that in the same process I installed 'MySql'Plugin & it is working fine.

Thanks in advance


texts:
I was trying to start using qgis2web plugin to make webgis for the first time, but when i click "create web map", i got this error message:

WARNING Traceback (most recent call last):
              File "C:/Users/UAV/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis2web\qgis2web.py", line 60, in run
                self.dlg = MainDialog(self.iface)
              File "C:/Users/UAV/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis2web\maindialog.py", line 160, in __init__
                QgsProject.instance().cleared.connect(self.reject)
            AttributeError: 'QgsProject' object has no attribute 'cleared'

Could anybody here help me to solve because I'm very new to this.


texts:
I'm using Openlayers v6.1.1 and I'm trying to get URLs of all tiles in the current map extent. I was trying to do it in the following way:

var tileUrlFunction = tileSource.getTileUrlFunction(),
    grid = tileSource.getTileGrid(),
    extent = map.getView().calculateExtent(map.getSize()),
    zoom = map.getView().getZoom();

grid.forEachTileCoord(extent, zoom, function (tile) {
    console.log(tileUrlFunction(tile.getTileCoord(), 1, ol.proj.get('EPSG:3857')));
});

But forEachTileCoord loop is never triggered. What I'm doing wrong?

I also tested different approach to get the URLs like that:

tileSource.on('tileloadend', function (evt) {
   console.log(getTileUrlFunction(evt.tile.getTileCoord(), 1, ol.proj.get('EPSG:3857')));
});

But this also doesn't give me an URL's - just undefined's.


texts: How can I edit and save feature attributes of a WMS layer in OpenLayers?

texts:
I am creating a data editing application in OpenLayers where the user will need to be able to edit the attributes of features (which are stored on a GeoServer installation). I have seen documentation for editing features in a WFS layer, but I haven't found anything detailing how to edit the attributes of a feature in a WMS layer. Is this feasible to do using a WMS layer?

I am currently retrieving data from GeoServer using the GetFeatureInfo request. I would like to load this data into a div, modify it, then save it back to the server. I haven't seen any examples of OpenLayers maps doing this, but I feel like it should involve a simple POST request to the server.


texts:
I have installed python 27 in one directory c:\program files\python27 and arcgis 10.2 in another directory c:\program files\arcgis, and within this folder there is c:\program files\arcgis\desktop10.2\arcpy. After opening up the arcgis map interface I can activate python IDLE and start using python commands.

But what if I want to use sublime text for arcpy instead?


texts: ArcMap maintaining raster cell size when using math?

texts:
I have a raster created from importing xy data (it has lat, lon, and a value field [WGS84]). The raster resolution is 60 meters as that's my required accuracy and for my value field I'd like to maintain around 6 decimal places.

I need to change this to integer in order to change the system to New Zealand Map Grid 1949 to align with my other data. However, when I use times in order to maintain my value field accuracy the cell size gets huge (10 km or so). The same occurs when I use int.

How do i maintain my 60 meter resolution to obtain my polygon for later transformation, to then change back to raster to continue with my work?


texts: How to add support for derivedCRSType in GeoTools?

texts:
I am trying to parse some GeoJSON files using GeoTools (24.2). The file contains non-EPSG CRS references like urn:ogc:def:derivedCRSType:OGC:1.0:image which causes my parsing to fail. I am using the typical gt-epsg-hsql but it doesn't seem to include that sort of CRS. Is there a jar file I could add to add support to this?
On the other hand I don't really make use of that CRS, the file contains a point with geographical coordinates and another one with image coordinates. If there is a way to just ignore that property that's equally fine with me.
Edit:
Here is an example feature (redacted as it's not my property)
        {
        "type": "Feature",
        "geometry": {
            "type": "Point",
            "coordinates": [
                6.835789963766803,
                53.46926328556648
            ]
        },
        "properties": {
            "xxxxx": {
                "result": {
                    "yyyy": "xxx",
                    "zzzz": "jjjjj",
                    "position": {
                        "type": "Point",
                        "coordinates": [
                            6.835789963766803,
                            53.46926328556648
                        ],
                        "crs": {
                            "type": "name",
                            "properties": {
                                "name": "urn:ogc:def:crs:EPSG:4326"
                            }
                        }
                    },
                    "image_position": {
                        "type": "Point",
                        "coordinates": [
                            22310.318371607515,
                            243069.73173277662
                        ],
                        "crs": {
                            "type": "name",
                            "properties": {
                                "name": "urn:ogc:def:derivedCRSType:OGC:1.0:image"
                            }
                        }
                    },
                    "detection_time": "2021-02-28T05:41:52.644086",
                }
            }
        }
    }

I also need to add that even without the image CRS it still fails, looks like GeoTools cannot parse properties that are not simple dictionaries of values?
Edit 2:
Here is some code, it's very simple though:
    FeatureJSON featureJSON = new FeatureJSON();
    FeatureCollection featureCollection = featureJSON.readFeatureCollection(new FileInputStream(geojsonpath));

The returned feature collection is empty but the geojson has features.


texts: geoprocessing tool : solving a vehicle routing problem

texts:
i am trying to implement a geoprocessing tool in my silverlight application witch use the "solving a vehicle routing problem", i tried to do that in ArcMap it work fine but i have no idea how to do it in a model builder and publish it in ArcgisServer.

need your help.


texts: Plot Circles / Polygons with KM radius from CSV in Google Earth Engine

texts:
I'm looking to plot a series of points with a set radius in google earth engine from a CSV. The table.draw function I'm currently using plots a pixel based radius, how can I convert this to a Km based radius?

After uploading my csv as an asset, importing into the console and using the following code to plot:

Map.addLayer(table.draw({color: 'FF0000', pointRadius: 5, strokeWidth: 5}), {opacity: 0.8});


texts:
I downloaded LISS4 satellite imagery of my study area, consisting of multiple tiles. I digitize some area from one of the tiles but I found out that the other tiles don't seem to align with the first one when I started to digitize an area which falls on the 2nd tile .

I am new with satellite imagery, I just want to know if there are any processes that need to be done before working with satelite imagery.

Is there any way of correcting this problem?

Note: Both the tiles are of same projections .


texts: How do I describe a special relationship between connected edges?

texts:
Consider this simple situation where three edges connect at a node:

I’d like to write a succinct and clear description of the relationship between A and B in such a way that differentiates it from the relationship between A and C. Something like “when traversing the node in the clockwise direction, A is adjacent? to B, but A is not adjacent? to C.” But it’s not really adjacency.

Said in a different way: imagine you’re standing on the node and you’re facing towards A. You start spinning yourself clockwise. The next edge you’ll come to is B, not C.

Is there a way to describe this relationship between A and B in a more succinct, formal, or correct way than I’ve written above?

It must be directional (one relationship of this type exists in the clockwise direction from A, and another exists in the counter-clockwise direction). And it must scale up to cases where more than three edges are connected at the node. Maybe it has something to do with routing? (I am thinking about this in the context of road networks.)

Two approaches I’ve tried already but haven't gotten far with:

9IM-like topology references: I’ve looked at the DE-9IM, and
even though I'm not a mathematician, I think I can still tell from the diagrams
and terms that it doesn’t cover this type of relationship. Neither
do I find it yet in the topology descriptions at the ESRI help or Oracle help. (Maybe there's something there but I'm just not finding it yet!)
Faces: I’ve played around with the fact that the face on the “north” side
of A might also bounded also be by B, but not C. However, as you can
see in the diagram here, that’s not always true. Imagine my diagram is an extract from a road
network where A and C are arterial roads and B is a short dead-end
road.

I suspect there may not be a single term for what I’m trying to say; at a minimum I’d like to be able to describe such a relationship in a simpler way than I’ve done above. This is a platform-independent question. Right now, I’m just looking for the right words. Later I will try to implement the concept in python (pyqgis or arcpy) on a shapefile, so any answers with that end point in mind will be particularly interesting, but not necessary.


texts: Why does ArcGIS Server not return data when querying joined fields?

texts:
I've been trying to help troubleshoot an issue that is very strange, and I think we've boiled it down to an issue with ArcGIS Server.

We have a map service that shows up in the services directory as:

Layers:
    Some-LayerA (0)
    ...
    Some-LayerE (4)

Tables:
    databaseName.SchemaName.TableName (5)

And the layer we are querying against looks roughly like this:

Fields:
    databaseName.SchemaName.LayerName.ObjectID (type: esriFieldTypeOID, alias: OBJECTID)
    databaseName.SchemaName.LayerName.SomeField (type: esriFieldTypeString, alias: SomeField, length: 3)
    databaseName.SchemaName.TableName.ProjectKey (type: esriFieldTypeString, alias: ProjectKey, length: 18)
    databaseName.SchemaName.TableName.Spent To Date ( type: esriFieldTypeDouble , alias: Spent To Date )
    ... More joined fields (No duplicate names or aliases, all strings, dates, doubles)

When issuing a query with only a where clause against the joined table, I can do something like this and have it bring results up:

Where: databaseName.SchemaName.TableName.ProjectKey = 'ValidProjectKey'

--> Returns 2 results
-----> databaseName.SchemaName.TableName.ProjectKey: ValidProjectKey

-----> databaseName.SchemaName.TableName.ProjectKey: ValidProjectKey

However, if I add some extra parameters, then I get no results at all:

Where: databaseName.SchemaName.TableName.ProjectKey = 'ValidProjectKey'
Out Fields: *

---> Returns 0 results

Apparently, if specifying the out fields manually, it will give results:

Where: databaseName.SchemaName.TableName.ProjectKey = 'ValidProjectKey'
Out Fields: databaseName.SchemaName.TableName.ProjectKey

---> Returns 2 results
-----> databaseName.SchemaName.TableName.ProjectKey: ValidProjectKey

-----> databaseName.SchemaName.TableName.ProjectKey: ValidProjectKey

But, if I change that second query to go against one of the layer's native (non-joined) fields instead, then it seems to return an appropriate number of results, but with all of the join fields set to null:

Where: databaseName.SchemaName.LayerName.SomeField = 'SomeValue'
Out Fields: *

---> Returns 514 results
------> databaseName.SchemaName.LayerName.ObjectID: 553
------> databaseName.SchemaName.LayerName.SomeField: SomeValidValue
------> databaseName.SchemaName.TableName.ProjectKey: null
------> databaseName.SchemaName.TableName.ProjectManager: null

------> databaseName.SchemaName.LayerName.ObjectID: 1938
------> databaseName.SchemaName.LayerName.SomeField: SomeOtherValidValue
------> databaseName.SchemaName.TableName.ProjectKey: null
------> databaseName.SchemaName.TableName.ProjectManager: null

Any clues why going against the joined table would lead to such different behavior? I don't have much control over the actual query, because it is being issued through geocortex.


texts:
I'm currently working on a dataset with 3.8m records. I have a very simple expression to update a field that I want to execute using the Field Calculator, while having the Attribute Table open.

Do I need to wait for the Loading Features prompt to complete for every Field Calculator operation? This is very slow on my computer.


texts: LULC changes in Google Earth Engine using Sentinel-2 data

texts:
I'm looking to classify land changes on river islands in the Congo River using Google Earth Engine with Sentinel-2 data. I need to obtain two cloud-free images for classification and then clip my images so that it is only a certain area in my tile that I analyse. I then need to classify and train the data using supervised classification and finally do a confusion matrix. I'm not sure how to start as I'm new to GEE and JavaScript.
This is the code I'm trying to use.
var START_DATE = '2016-7-4';
var END_DATE = '2021-3-3';
var MGRS_TILE = '33MYU';
// Malebo pool should be '33MWR', but for some reason it only shows the top corner of the tile?

// get tile
var s2_unsorted = ee.ImageCollection('COPERNICUS/S2')
                    .filterMetadata('MGRS_TILE', 'equals', MGRS_TILE)
                    .filterDate(START_DATE, END_DATE)
                    .filterBounds(Islands);
var s2_sorted = s2_unsorted.sort('CLOUDY_PIXEL_PERCENTAGE');

// Convert the collection to a list and get the number of images.
var size = s2_sorted.toList(1).length();
print('Number of images: ', size);

// Get the number of images.
var count = s2_sorted.size();
print('Count: ', count);

// Get the date range of images in the collection.
var dates = ee.List(s2_sorted.get('date_range'));
var dateRange = ee.DateRange(dates.get(0), dates.get(1));
print('Date range: ', dateRange);

// Get statistics for a property of the images in the collection.
var sunStats = s2_sorted.aggregate_stats('SUN_ELEVATION');
print('Sun elevation statistics: ', sunStats);

// Sort by a cloud cover property, get the least cloudy image.
var image = ee.Image(s2_sorted.sort('CLOUD_COVER').first());
print('Least cloudy image: ', image);

// Limit the collection to the 10 most recent images.
var recent = s2_sorted.sort('system:time_start', false).limit(10);
print('Recent images: ', recent);

var least_cloudy_image = ee.Image(s2_sorted.first());
//var least_cloudy_image = ee.Image(s2_sorted.mean());

// Visualization parameters
var visParams = {bands: ['B8', 'B4', 'B3'], max: 3048, gamma: 1};
var visParams_ndvi = {min: -0.2, max: 0.8, palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
    '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'};

//Sandy Beach area
Map.setCenter(16.86, -1.25, 15);
// Malebo pool should be below

//Map.centerObject(least_cloudy_image);
Map.addLayer(
  least_cloudy_image,
  {bands:'B4,B3,B2', min:0, max:3048, gamma: 1},
  'least cloudy image'
);
print('least cloudy image', least_cloudy_image);

// Calculate NDVI
var image_ndvi = least_cloudy_image.normalizedDifference(['B8','B4'])
.filterBounds(Islands);

// Map results
//Map.centerObject(image,9)
Map.addLayer(least_cloudy_image,visParams,'Sentinel-2 False Color Infrared')
Map.addLayer(image_ndvi,visParams_ndvi,'Sentinel-2 NDVI')

var newfc = Forest1.merge(Agriculture1).merge(Burnt1);
print(newfc);

var bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7','B8'];

var training = image.select(bands).sampleRegions({
  collection: newfc,
  properties: ['landcover'],
  scale: 30
});
print(training);

var trainingNoNulls = training.filter(
  ee.Filter.notNull(training.first().propertyNames())
);

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

var classified = image.select(bands).classify(classifier);

Map.centerObject(newfc, 11);
Map.addLayer(classified,
{min: 0, max: 2, palette: ['green', 'blue', 'red']},
'classification');
Map.addLayer(newfc);


texts: Can't install support for ECW in QGIS 3.6 / 3.4 on Ubuntu 18.04

texts:
I've been trying to install support for ECW so I can load them in latest QGIS 3.6 (also tried in QGIS 3.4.4)

Having tried the following tutorials:

https://wiki.openstreetmap.org/wiki/ECW   (gdal 2.2.3)

https://gist.github.com/1papaya/568c4580b1909071696c1cb119101823   (gdal 2.3.1)

Both appear to work and report back ECW support against gdal

  gdalinfo --formats | grep ECW
  ECW -raster- (rw): ERDAS Compressed Wavelets (SDK 3.x)
  JP2ECW -raster,vector- (rw+v): ERDAS JPEG2000 (SDK 3.x)

However when opening QGIS I notice that key UI is missing - the little triangles next to all the processing toolbox processes are missing, and the whole of QGIS becomes unstable - for example it crashes to desktop when selecting plugins within the plugin manager.

I then have alot of trouble recovering QGIS to get the UI back - full removal and purge of all dependencies + QGIS reinstall does not restore it! have to carry out full ubuntu reinstall to get it all back, which is not ideal.

Is gdal completely independent per version install. as above tutorials focus on particular gdal download which you then appear to manipulate?

QGIS 3.6 appears to use gdal 2.2.3

QGIS 3.4.4 appears to use gdal 2.4.0

so far I've carried out above tutorials first then installed QGIS, should I perhaps install QGIS first then attempt to adjust the GDAL that QGIS installs as a dependency?


texts: Con function loop not calculating for all fifty rasters in list

texts:
I have been trying this for hours and I know it should be easy but I can't make it work yet.
I have a list of 50 rasters that I need to change the NoData to "0" values. I have the right arcpy function right but I can't get the loop to work yet (I'm a beginner)
So far what I have is"

env.workspace = "My raster folder"
rasterList = arcpy.ListRasters("*.img")
print rasterList
#works
OutDir= "My output raster folder" + raster [0:7] + ".img"
#Here is my loop that won't work
for raster in rasterList:
    OutRaster = Con(IsNull(raster), 0, raster)
outCon.save(OutDir)
#this prodices only one raster out of the 50 rasters in the List.

I need to do the same iteration for all the 50 rasters and making sure that on each iteration both of the rasters entered in the Con function is the same each run, This is something like this:

Con(IsNull(raster1), 0, raster1)
Con(IsNull(raster2), 0, raster2)
Con(IsNull(raster50), 0, raster50)

Can someone help me?


texts: How to define a projection for a file with geo points?

texts:
I have a HDF file that I want to convert to tif or ENVI file.

I know that the +Proj of the original file is:
"+proj=cea +lon_0=0 +lat_ts=30 +ellps=WGS84 +units=m"

so:

  gdal_translate -of "ENVI" -a_srs "+proj=cea +lon_0=0 +lat_ts=30 +ellps=WGS84 +units=m" HDF5 file.h5  output.img

but I found that the projection was not defined but instead I got this:

       ENVI
     samples = 1000
     lines   = 500
     bands   = 1
     header offset = 0
     file type = ENVI Standard
     data type = 4
     interleave = bsq
     byte order = 0

My question how to deal with geo points?


texts: How to divide map into square grid using ArcGIS Desktop?

texts:
I want to divide a map in square grid using ArcGIS Desktop.

Is it possible?


texts: Multiple Line in single stacked text box using ArcGIS Desktop?

texts:
I am having difficulties with creating a stacked label based on different soil boring depth.

Below is my attribute table:

I want to display something like this stacked in 1 text box:

Every time I use my script I get multiple text boxes and then have to copy and paste them in to 1 text box.

My VBScript is:

[GIS_SB_Ben] &vbCrLf & [GIS_SB_B_1] +" - "+ [GIS_SB_B_2] &vbCrLf & "Benzene = "+ [GIS_SB_B_3]

Any suggestions?