texts: How to use SAGA Raster processing in processing scripts (QGIS3)

texts:
I want to make a processing script which uses the saga:inversedistanceweightedinterpolation algorithm (Interpolating a point layer - magnetometry measurements - to a raster. The SAGA runs so far but the result could not be loaded by QGIS. I have the following script:

from PyQt5.QtCore import QCoreApplication
from qgis.core import (QgsProcessing,
                   QgsProcessingAlgorithm,
                   QgsProcessingParameterFeatureSource,
                   QgsProcessingParameterRasterDestination)
import processing

class ExampleProcessingAlgorithm(QgsProcessingAlgorithm):

   INPUT = 'INPUT'
   OUTPUT = 'OUTPUT'

   def tr(self, string):
       return QCoreApplication.translate('Processing', string)

   def createInstance(self):
       return ExampleProcessingAlgorithm()

   def name(self):
       return 'test_IDW'

   def displayName(self):
       return self.tr('Test IDW')

   def group(self):
       return self.tr('Example scripts')

   def groupId(self):
       return 'examplescripts'

   def shortHelpString(self):
       return self.tr("Example algorithm short description")

   def initAlgorithm(self, config=None):
       self.addParameter(
            QgsProcessingParameterFeatureSource(
            self.INPUT,
            self.tr('Input layer'),
            [QgsProcessing.TypeVectorPoint]))
        self.addParameter(
            QgsProcessingParameterRasterDestination(
                self.OUTPUT,
                self.tr("Output Raster"),
                None, False))

    def processAlgorithm(self, parameters, context, feedback):
        output =  processing.run("saga:inversedistanceweightedinterpolation", {
            'DW_BANDWIDTH' : 1,
            'DW_IDW_OFFSET' : False,
            'DW_IDW_POWER' : 1,
            'DW_WEIGHTING' : 1,
            'FIELD' : 'READING',
            'SEARCH_DIRECTION' : 0,
            'SEARCH_POINTS_ALL' : 0,
            'SEARCH_POINTS_MAX' : 20,
            'SEARCH_POINTS_MIN' : -1,
            'SEARCH_RADIUS' : 1,
            'SEARCH_RANGE' : 0,
            'SHAPES' : parameters['INPUT'],
            'TARGET_DEFINITION' : 0,
            'TARGET_OUT_GRID' : parameters['OUTPUT'], #'memory:',
            'TARGET_TEMPLATE' : None,
            'TARGET_USER_FITS' : 0,
            'TARGET_USER_SIZE' : 0.25,
            'TARGET_USER_XMIN TARGET_USER_XMAX TARGET_USER_YMIN TARGET_USER_YMAX' : None
        }, context=context, feedback=feedback)['TARGET_OUT_GRID']
        return {self.OUTPUT: output}

I get the following error after execution:
Loading resulting layers
The following layers were not correctly generated.C:/Users/Stefan/AppData/Local/Temp/processing_03bc47ec83e54650bd96292f861633b5/11bb8145ca214423b41e7c1c33c3965c/OUTPUT.tifYou can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.
the Log Message Panel shows: 2019-01-30T18:06:45     INFO    Task complete : Executing “Test IDW”


texts: I have about 16000 points with speed on the roads and some roads (lines). I need to join points to lines and their attributes so that there are speed values from points in the line attributes. I have no idea what to do.


texts: QGIS Time Manager animation does not run

texts:
I'm encountering some problems with Time Manager that after adding layer to Time Manger settings the animation doesn't run. Nether the play button or trough slider. I'm following steps of the video (https://www.youtube.com/watch?v=nHrFOPf1UGw&t=283s). What should be the field type for changing attribute?

Sample data:

x       y       aika                haitta10

550144  7116401 01-01-2001 00:00:00 0
550144  7116401 02-01-2001 00:00:00 2
550144  7116401 03-01-2001 00:00:00 4

QGIS 2.14.10, Time Manager 2.4, Windows 7


texts: Finding line with minimum distance to point and update attribute of points?

texts:
I am currently working with QGIS 2.18 and I have two layers: one layer is filled with lines which represent streets. The second layer describes the location of sensors which are represented as points like in the image below:

Each point has its own attribute values like followed:

As you can see, the attribute edgeid is still empty and I want to fill it with the id of the nearest street. Each line in the street-line-layer has an edgeid as followed:

Now I will have to determine the street with the shortest distance to the point and store its id into the edgeid-attribut of my point. The distance between point and line should be calculated as illustrated (d=distance; AB is the street vector and P is the point):

I have already tried the Distancematrix in QGIS but it seems that it can only be done between points. And I do not want to convert the lines into points.

Another thing I tried was the v.distance()-algorithm (GRASS) but I couldn't get it to work. I tried the following configurations for executing the algorithm:

For your better understanding: The tiger:graph is the line-layer and the EVIS_Stat_Sensorik_... is the point-layer. I want to calculate the nearest line (or edge) for each point and write its id into the attribute "edgeid" of the point. Is this even possible with this function? As I am quite new to this tool, I do not really know how to debug this.

Is there any other possible way to do this?


texts: Geotools produces different GEOGCS value for ESRI WKT

texts:
ArcMap creates following contents in .prj file for WGS 84.

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

While Geotools creates following when used ESRI citation, as per documentation: (http://docs.geotools.org/stable/userguide/library/referencing/crs.html)

GEOGCS["WGS 84", DATUM["D_WGS_1984", SPHEROID["D_WGS_1984", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS["Geodetic latitude", NORTH], AXIS["Geodetic longitude", EAST], AUTHORITY["EPSG","4326"]]

Can anyone please tell me that why value of GEOGCS is appearing as "WGS 84" with geotools and not "GCS_WGS_1984". Because of this difference, ArcMap is not recognizing the projection.

Or is there any other way that I can create .prj file which will reflect  GCS_WGS_1984?

Geotools version: 13.5

Sample code:

Formattable f = (Formattable) CRS.decode("EPSG:" + epsg.getCode(), true);
String wkt = f.toWKT(Citations.ESRI, 0);

CoordinateReferenceSystem crs = CRS.parseWKT(wkt);


texts: Making leaflet handle tiles that have nodata values transparent?

texts:
I have a tif that I cut out the two sections of and merged them together into one tif and then processed it.  Shown below is the original and the processed cuts (not to scale).  Between the two processed cuts are nodata values. My original mosaic is already tiled.  How would I go about displaying the tiled cuts that I have over the original tiled mosaic so I can see both of them at once in leaflet?

Update: For my example above if I want to overlay the red images over the mosaic above it all that is needed to add it to the map layers variable in leaflet and make sure that the red images have no data values for the pixels in between them.


texts: Missing columns on WFS (GeoServer/MongoDB)

texts:
We are currently checking up on GeoServer over MongoDB on our dev-environment. We are having some difficulties: we have a collection, containing different documents which answer to the same structure ("schema-like"). When querying a WFS service from GeoServer - some of the columns are missing.  We discovered that if one of the "columns" on one of the documents contains null value - the whole column does not appear on the query response.

Is there some configuration we are missing on the GeoServer to return null values on WFS response?

Even if one instance of an attribute value equals null, the whole column is missing. I don't see how DescribeFeatureType may help as many Q&aA already suggested over the web.


texts: How to get the EGM96 and EGM2008 elevations of a list of WGS84 lat/lng coordinates from a file using GDAL?

texts:
I'm trying to do the following:
Given a list of WGS84 latitude / longitude pairs in a text file (one pair per line), I'd like to compute the EGM96 and EGM2008 elevations corresponding to those points using GDAL. Unfortunately, it seems gdalwarp is usually used to convert from one SRID to another, not the use case I'm looking for.
I'm sure there must be a way - I'd prefer to use GDAL for this task (as I've been recommended to use it), but I'm happy to consider other (command-line) alternatives to solving the same taks.
EDIT: MATLAB gives me the exact API I'm looking for
https://www.mathworks.com/help/aerotbx/ug/geoidheight.html
But I don't have access to MATLAB and it seems GDAL does have the same capabilities, if only I could find it!


texts: Returning top three highest numeric values using Arcade

texts:
I was introduced to Arcade a year or so ago and haven't implemented it into any of my projects as I hadn't had a need to, but I think I might have a use for it, if it can work.
I have a feature class that summarizes all of our mosquito trap data for our season. And in this dataset we have roughly 30 types of mosquitoes we catch throughout the season.  I wanted to create a pop up that only shows the top 3 mosquito species (in terms of amount caught) for that particular trap.
Is this possible?
Before I research this any further I'd figure I'd ask the pros!


texts: Problem in distributing GDAL in stand alone application

texts:
I'm desperate with the distributing problem of gis packages each time. It is the third day that I'm desperately trying to make an .exe using Pyinstaller to a small application I developed using GDAL and PyQt5 with Python 3.7.
I use Pycharm IDE with venv based on the qgis_ltr 3.4 Python as a project interpreter.

First to make the dev environment I had to run Pycharm with a batch that sets environments variables. Now that I want to make the install, I know I must reference GDAL somehow, after trying many alternatives ( installing gdal in a venv )  Nothing seems to work.

I get the exe file but that gives an error message not founding the GDAL module.
Is there any simple way to correctly generate the .exe file?


texts: QGIS 3.10 send incorrect ArcGIS MapServer requests (pkk.rosreestr.ru)

texts:
I am trying to add a web service (ArcGIS MapServer) in QGIS 3.10 here is the url:

https://pkk5.rosreestr.ru/arcgis/rest/services/Cadastre/BordersGKN/MapServer/

I see a list of layers, but QGIS gives the following error:

  Network error: Error transferring

  https://pkk5.rosreestr.ru/arcgis/rest/services/Cadastre/BordersGKN/MapServer/export?bbox=8097219.581788,8629436.355556,8242578.042442,8741170.443531&size=977,751&format=PNG32&layers=show:33&transparent=true&f=image
  - server replied: Bad Request

Through the browser, I examined the requests that are sent when viewing the original site. I noticed that QGIS requests are different from requests received through the browser.

Requests sent through the browser (successful):

https://pkk5.rosreestr.ru/arcgis/rest/services/Cadastre/BordersGKN/MapServer/export?layers=show%3A&dpi=96&format=PNG32&bbox=8061966.247258484%2C8688138.383019326%2C8140237.764222398%2C8766409.899983242&bboxSR=102100&imageSR=102100&size=1024%2C1024&transparent=true&f=image

1. Is it true that QGIS 3.10 does not work correctly with ArcGIS web services (MapServer)?

2. How to fix this behavior?


texts: What causes 503 Errors when loading tiles from CartoDB?

texts:
I'm currently having great difficulty loading my visualizations onto my webapp from the CartoDB cdn. A majority of the tile requests are returning 503 errors with this response :

<?xml version="1.0" encoding="utf-8"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">

 <title>503 backend read error</title>

Error 503 backend read error

backend read error

Guru Meditation:

XID: 2031602388

Varnish cache server

Is this perhaps due to bad data migration?


texts: Polygons not merging in QGIS when dissolving?

texts:
I am new to QGIS and am trying to dissolve polygons in my shapefile using one of the fields in the attribute in QGIS 2.8.5. The process runs, and it does add a new map but nothing seems to have changed in the new map and the polygons chosen on the basis of the field are not merged.

If this is a known problem, how can I work around it?


texts: Mapbox GL JS center map on click near the boarder

texts:
The code below centers the map on a click, how could this code be edited it so it only centers the map if the click is near the boarders (say within 100px) of the map?

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

    <link href='https://api.tiles.mapbox.com/mapbox-gl-js/v0.26.0/mapbox-gl.css' rel='stylesheet' />


texts: Modifying attributes of a shapefile using Java

texts:
I have a polygon vector shapefile in which each feature has an integer attribute. I want to update this integer attribute (for all features) programmatically using Java on regular intervals of time.

How can I do this?

EDIT: I am using GeoServer to display my shapefile over the network. Since one of the attributes is continuously changing on a daily basis, I need to reflect that by modifying the shapefile. I want these changes to further reflect in GeoServer without restarting the GeoServer.


texts: Map Algebra expression in a model fails on open workspace

texts:
I have model which runs 2 processes that create rasters and a third process that (attempts to) add the two of them using Map Algebra. The expression is very simple, just "%Raster_A%" + "%Raster_B%".

This third process aborts with ERROR 000539: Error running expression: rcexec() <type 'exceptions.RuntimeError'>: ERROR 010167: Could not open workspace .

The help file for Error 000539 : Error message from Python does not appear relevant, this isn't a problem with python syntax or type use. There is no help page for #010167, but  the "not open workspace" is self explanatory.

I haven't found any info on rcexec() but name and context suggest it's a method for shelling out to an external python process and getting something back. It's reasonable to assume the external process doesn't know it's being executed from a model and that the workspace lock it's encountering is from the parent, and that this is okay. Or, maybe rcexec is not encountering a lock and what's really having trouble with is finding the workspace. I have no way of knowing (I don't think).

The model has an Output Workspace input parameter, which the two antecedent processes use just fine. I tried connecting the Output Workspace input parameter to the map algebra expression as Enviroment > Workspace > Current Workspace to no effect.

Any ideas how to track the source of this error down and fix it?

update: Aha! The environment is not clean. I edited the model and deleted the Current and Scratch workspace settings, and they reverted to the default C:\Users\matt\ArcGIS\default.gdb. At that point the error changed from can't open to invalid workspace:

ERROR 000539: Error running expression: rcexec() <type 'exceptions.RuntimeError'>:
ERROR 000875: Output raster: T:\env298.gdb\A_plus_B's workspace is an invalid output workspace.

The curious thing is that the workspace specified in the error was the model's workspace at one time. So rcexec() is not respecting the value set in the model's environment properties, nor the process's workspace environment property.

Now I just need to figure out how to really reset the environment to defaults.


texts:
I'm having an issue where the vector layer I am using is not working with with the basemaps (see figure). I know the issues is because the vector data is trying to draw a shape from the end of one side of the map to the other. The original data is from here (https://datasets.wri.org/dataset/tropical-coral-reefs-of-the-world-500-m-resolution-grid) and is the most complete data set I can find on coral reefs (I know natural earth has a set but it isn't as complete). All I want to do is display the shape/vector file (.shp) on the world basemap.

I'm very new the GIS so hopefully its just a simple thing that I have overlooked but I've been trying for hours and I can't even think of the right keywords to search.


texts: How to assign OGR geometry a spatial reference when writing to a new layer using the Python API?

texts:
I'm using the Python OGR API to add new features to a PostGIS layer. However, despite setting the geometry spatial reference, when the geometry is written out the spatial reference is missing. This is a problem as the destination table has check constraint forcing geom_srid to be a particular type.

For example:

ogr.UseExceptions()

# Create a feature in an existing layer (where layer is a PostGIS table)
feature = ogr.Feature(layer.GetLayerDefn())

# Create spatial reference
out_srs = ogr.osr.SpatialReference()
out_srs.ImportFromEPSG(27700)

# Assign to geometry
geom.AssignSpatialReference(out_srs)

# Set feature geometry
feature.SetGeometry(geom)

# Create the feature
layer.CreateFeature(feature)
feature.Destroy()

Returns:

Error: new row for relation "table" violates check constraint
"enforce_srid_geom"
Command: INSERT INTO "table" ("geom", "attr1") VALUES (GeomFromEWKT('SRID=-1;POINT (5000, 1000)'::TEXT), 'row1')

At this point it is clear that OGR is converting the geometry to an EWKT string without pulling in the spatial reference attribute.

Specific questions:

1) Does anyone know how I can force the geometry output from SetGeometry() to include an SRID when it is converted to EWKT?

2) Has anyone found this with the CPP bindings?

3) Is this a bug in OGR?
From OGR source the output string is being created in ogrpgtablelayer.cpp where nSRSID is being set to -1 if no valid SRS input is found.

Many thanks,

Tom


texts: Restricting portions of road network from analysis in ArcGIS Network Analyst?

texts:
I'm trying to sort out network analysis in ArcGIS 10.

This example of finding shortest distance between two points looks something like that :

Where squares are two locations, and the green line is an output from Closest Facility tool.

What I need to do, is to restrict motorways and slip roads from the analysis. So the new route should avouid the red segments:

What is the best strategy to proceed in such case?

Should I edit my network and delete all the segments I do not want to have in analysis?

Or perhaps use linear Barriers?


texts: MMGIS not able to connect to maps.googleapis.com

texts: I've been trying to geocode a few addresses from a .csv file to a shapefile.  Been using MMQGIS on Qgis, but, even though I seem to have everything set up properly, I continually get the error message that there was a failure to connect to maps.googleapis.com and therefore it fails to geocode any of the addresses from my excel spreadsheet.  Any suggestions on how I can fix this?


texts: Rotating symbols with data frame in ArcGIS Pro?

texts:
I have a Map Series in ArcGIS Pro that was imported from ArcMap (with Data Driven Pages) to produce several different orientations. In Pro, symbols are not rotating with the data frame-- they are always oriented North up. Symbols with a rotation value in the attribute table are being rotated by these values, but nothing is rotated with the data frame.

Is there something I need to do (along the lines of checking "rotate symbols with data frame in Utilities for ArcMap) to get this to work?

I am totally new to ArcGIS Pro.


texts: OpenLayers load vector data from web API

texts:
We have a web application that displays a map with lots of point and polygon features. The features are loaded as WKT from our web API along with other data.

Currently we just fetch all the features when the web first loads, but naturally this leads to a lot of unnecessary data being transferred. The polygon geometries in particular take a lot of bandwidth, and most of them will never be displayed to the user since the polygon layer is only visible when zoomed in.

Is there a way to have OL automatically make API queries for only the necessary extent?

Ideally we'd like OL to make an API query for the current extent and visible layers and merge the resulting data to the already loaded features.

We use OL 5.


texts: Can FME Desktop (node-locked) be installed on an Amazon EC2 or VPC instance?

texts:
I did do some searching in the FME help and the legal notices, but I'm not getting any hits (probably due to vocabulary). Are there any legal restrictions for FME Desktop that would prevent me from installing onto an EC2 or VPC Amazon instance?  I'm not sure if (or how much) that would compete against the FME Cloud offering.  I'm mostly interested for test and evaluation, but I am also curious if anyone is using an EC2/VPC as a permanent home for FME Desktop.

This another question that is inspired by this closely related question:
Is there an optimum RAM configuration to use for FME Desktop?


texts: Con statement in raster calculator ArcGIS 10.3

texts:
I have four rasters (A, B, C, D) based on which my classification works. The con statement I have been using is posted below.

My problem is when A_01 is 2. Area when A=2, if value of B is less than 3, it categorizes the ares under 2, but when B is between 3 and 12, and value for C is exactly 0.24, it is categorizing the area under 3.

For eg: any area where A=2, B=3.56 and C=0.24, I assume it should fall under category 2, since value for B is less than 12 and C is exactly 0.24, but it is being categorized as 3. Can anybody help me figure the error in this statement?

Con(("%A_01.tif%" == 1) & ("%B_01.tif%" > 12), 3,
Con(("%A_01.tif%" == 1) & ("%B_01.tif%" <= 12), 4,
Con(("%A_01.tif%" == 1) & ("%D_01.tif%" == 1), 4,
Con(("%A_01.tif%" == 2) & ("%B_01.tif%" <= 12) & ("%C_01.tif%" >= 0.24), 2,
Con(("%A_01.tif%" == 2) & ("%B_01.tif%" > 12), 2,
Con(("%A_01.tif%" == 2) & ("%B_01.tif%" >= 3) & ("%B.tif%" <= 12) & ("%C_01.tif%" < 0.24), 3,
Con(("%A_01.tif%" == 2) & ("%B_01.tif%" < 3) & ("%C_01.tif%" < 0.24), 4,
Con(("%A_01.tif%" == 2) & ("%D_01.tif%" == 1), 4,
Con(("%A_01.tif%" == 3) & ("%D_01.tif%" == 0), 2,
Con(("%A_01.tif%" == 3) & ("%D_01.tif%" == 1), 4,
Con(("%A_01.tif%" == 4) & ("%D_01.tif%" == 0), 1,
Con(("%A_01.tif%" == 4) & ("%D_01.tif%" == 1), 4))))))))))))


texts: ol is not defined Error while trying to use Openlayers 3 in Smarty template

texts:
I'm trying to integrate an Openlayers 3 map into an Smarty-template.
The parts of the code work fine when used in a single html-(php)-page. When I integrate this code in the Smaty template, the sites displays without a map and in Firebug I get an 'ReferenceError: ol is not defined'.
All the libraries are loaded correctly and I verified, that the order of the loaded files is the same as in my working single file.

Here parts of the Smarty-Template (I left the Smarty part away, but if needed, I could provide also the entire code).
The css - files as well as the .js-files are loaded via href / src, so they are not enclosed by {literal} {/literal}. This shouldn't be the problem, since in the original Smarty-code also other css- or .js-files are implemented that way.

I use a similar architecture with the google maps api, which work with no probs, but I wan't to change to OpenLayers. I don't see, why OL isn't accessible?
The Errors I get relays to the ol-objects defined in the ol3-layerswitcher.js as well as in the mapOL.js.

<!DOCTYPE html>

// other head-tags and Smarty-code

    {* ------------------- Add js- and css-files for map-display ----------------------- *}

    <link rel="stylesheet" href="./_css/ol.css" type="text/css">

    <link href="./_css/fa/css/font-awesome.min.css" rel="stylesheet">

    <link href="./_css/ol3-layerswitcher.css" rel="stylesheet">

    <link href="./css/custom.css" rel="stylesheet">
    <link href="./css/map.css" rel="stylesheet">

    {* ------------------- End Add js- and css-files for map-display ----------------------- *}

{* ------------------- Body starts here ----------------------- *}

// Navigation code from Smarty

 {* -------------- Add map div -------------------------- *}

 {* -------------- End Add map div ---------------------- *}

// other code from Smarty


texts: How to stop transparent MapServer map images becoming opaque when a tif image is added?

texts:
How do you add a tif image to a MapServer .map file and keep the MapServer image transparent beyond the extent of the tif image?
I have a MapServer .map file with the following OUTPUTFORMAT declarations:
  OUTPUTFORMAT
    NAME "png24"
    MIMETYPE "image/png"
    DRIVER "AGG/PNG"
    EXTENSION "png"
    IMAGEMODE RGBA
    TRANSPARENT TRUE
    FORMATOPTION "DPI=100"
  END # OUTPUTFORMAT

  OUTPUTFORMAT
    NAME "print"
    MIMETYPE "image/png"
    DRIVER "AGG/PNG"
    EXTENSION "png"
    IMAGEMODE RGBA
    TRANSPARENT TRUE
  END # OUTPUTFORMAT

When I add a tif image as a raster layer, the mapfile becomes opaque and I can't see any background colour or image through it.  This opaque change happens outside the extent of the raster image.  My raster layer declaration is:
  LAYER
    COMPOSITE
    END # COMPOSITE
    DATA "D:\\raster.tif"
    METADATA
      "epsg"    "28356"
      "coordsys_name"   "GDA94 /  MGA zone 56"
    END # METADATA
    NAME "my raster image"
    PROJECTION
      "proj=utm"
      "zone=56"
      "south"
      "ellps=GRS80"
      "towgs84=0,0,0,0,0,0,0"
      "units=m"
      "no_defs<>"
    END # PROJECTION
    STATUS OFF
    TILEITEM "location"
    TYPE RASTER
    UNITS METERS
  END # LAYER

I have tried adding an alpha band to my raster, and, I have tried making a PROCESSING "BANDS = 1,2,3,4" declaration but these don't appear to help.


texts: Dynamic labeling looks awful in flex map

texts:
I published a service to ArcGIS Server 10.1 from an .mxd with no halos around the labels and created a flex application using the flex viewer. When I browse to the map here and zoom in on an individual trail segment, the labels become very grainy and appear to have a white halo around them, making them very difficult to read. It becomes exceptionally noticeable when imagery is used as a basemap.

Is there any thing I can do to stop the labels from getting all messy looking or is this just how they look when using flex?


texts: Moving attribute data to the objects geometry

texts: I have contour lines with the elevations of them in an attribute called 'CONTOUR' however these do not actually elevate the lines when I bring them into AutoCAD.  To do this I want to copy and paste the 'CONTOUR' attribute information and paste it into the 'ELEVATION' attribute as shown in the table.  When I Join the excel sheet to the shapefile it doesn't bring in the elevations though and I am not sure why.


texts: Overlay raster image over a polygon shapefile

texts: I have a raster image which is the coordinate system WGS 1984 UTM 44 of the inundation map and I have a polygon shapefile of same coordinate system WGS 1984 UTM 44 containing the shapefiles of districts. When I am superimposing them, they are not able to superimpose on each other even though the area covered by raster layer is enclosed within the shapefiles of districts.


texts: 'OnStartup()' method in ArcGIS extension no longer running when ArcGIS starts up

texts:
I have a couple 'large' extensions for ArcGIS that have been in development for a couple of years now. Recently I've found that the "OnStartup()" method in the Extension class is no longer getting hit when ArcGIS is launched. I have a feeling that something has gone with my COM registration.

Note that this is an extension, not an add-in.

The major issue I have run into is my default workspace for storing temporary files is being setup from this "OnStartup()" method (checks to see if the workspace already exists and has files in it, if not it will create the workspace at a default path). So, when I attempt to run a tool from the extension it crashes as it cannot find the path it expects to store temporary files. Interestingly, if I open the Customize->Extension menu in ArcGIS and un-check/check the extension the "OnStartup()" method runs. So, it's almost as if ArcGIS is treating the extension as if it were "off", even though it is checked in the extension menu.

I've gone way back and compared the extension class in source control and cannot find any changes that could have caused this error. This problem may have occurred when I installed the latest service pack, but I can't be sure.

Has anyone else encountered this? Any help would be greatly appreciated!


texts: QGIS 2.12 How can I assign the attributes of the polygons of the layer A containing the geometry of layer B?

texts:
In other words, I have two shapefiles:

Layer A: not-intersected Areas (1 attribute: ID);
Layer B: Buildings (Several attributes).

I want to identified by the related ID attribute all the buildings of layer B included in the area of Layer A.

If the area represents the boundary of several towns and the ID attribute represents the name of each town, I want to assign to each building an attribute representing which town they belong to.


texts: How to get WFS GetCapabilities information as JSON

texts:
I can get a WMS GetCapabilities response via AJAX call and I can convert response to JSON using ol.format.WMSCapabilities tool.

But I need to make a WFS GetCapabilities query via a AJAX request and get  a JSON result.

How can I do it?

My GetCapabilities operation metadata is like this:


texts: Merging and generalising polygon data with ogr2ogr and SQL

texts:
I have some building data that is already generalised to blocks of buildings rather than individual buildings, and I want to create another set of data from it that is generalised again to urban areas.

Here's a sample of what I have:

You may be able to make out that there are spaces between the blocks, which account for roads/parks/etc. So the polygons I have are disconnected, so a simple merge/dissolve won't do quite what I want.

What I want to end up with is something like this (excuse the poor Photoshop skills):

Basically, I want to end up with polygons like the ones I have crudely drawn in black. So the buildings data is generalised and merged into urban areas where the buildings are pretty dense. I don't mind if buildings in sparse areas are either drawn as-is or ignored entirely, either is fine (but the latter is preferable).

I am using ogr2ogr to preprocess all of my data along with the -sql flag to filter my data and create new types (for example, I've processed point data representing roundabouts into polygon data, removed unwanted metadata etc), so if I can achieve my goals using SQL rather than QGIS that would be preferable.

Any ideas?


texts: How to extract bounding boxes of cities from OSM?

texts:
For simple geocoding purposes I want to extract all cities from a single country of openstreetmap.

They only queries I want to do are:

ZipCode -> (List of) BoundingBox of City/Village
City/Villagename -> (List of) Boundingbox of City/Village

I'm not interested in streetnames, rivers etc.
Is this possible using osmfilter. Can anyone point me towards the tags of the data i'm looking for?


texts: Technique for splitting raster image into tiles

texts:
I have written a standalone application which renders a globe similar to GoogleEarth/WorldWind which is currently able to display raster imagery from WMS servers. The native geodesy of this application for rendering purposes is WGS84.

The tiling system employed is similar to this: WorldWind tiling system

I now want to add support for importing geo-referenced rasters (supported by GDAL) directly.

I can successfully import and render world files which have a direct 1:1 mapping to my tiling scheme. I'm using the GDALBand ReadRaster method to read the raster into memory.

The problem I am facing is that I am currently stuck at trying to figure out how to import rasters that do not map nicely to my tiling scheme (e.g a tif that has it's coordinates not match up nicely to any tile corner/boundary).

I know WMS can deal with this as I can get it to render rasters with WMS that do not nicely match up to my tile system. It looks like it just pads everything out with a default colour.

I'm interested to learn of possible techniques to deal with this problem.

Will it require padding out the gaps?

Do I need find the level which will result in tiles smaller than my raster and then split up the raster that way? This still results in the same question of what do I do with the corners and edges of the raster which are only partially covered by a tile.


texts: Drawing pie charts with vectors with OpenLayers?

texts: Has any one got any examples of drawing pie charts with vectors as markers with OpenLayers?


texts: Extracting information from buffer in ArcMap?

texts:
I have a buffer in ArcMap 10 extending from a polyline, with numerous 5 km intervals (multiple buffer tool). I also have many layers already in raster format in ArcMap, including slope angle, geology, precipitation data etc.

Basically I want to be able to get a value for each of these layers for each buffer polygon. The area is a transect through Nepal, and thus there are approximately 30 separate buffers extending 155 km towards the Himalayas.

Ultimately I am aiming to produce a graph with distance up Nepal (x-axis) and then elevation (y-axis) and another variable such as precipitation (opposite y-axis).


texts: Creating a drop down list in ModelBuilder

texts: I have a model that exports GIS data to CAD.  The model has two nested models, the first collects the data to be exported and the other exports it.  The export to cad tool has a seed file parameter and I would like to make this a drop down box of two seed(.dwg) files that I have.  I am able to make it a parameter but not a drop down list.  I thought I could go into properties-->parameters-->filter but nothing happens when I click on the "filter" box which is set to none.  What am I missing?


texts: How Rational Polynomial Coefficients(RPCs) are calculated? Need references

texts:
I wonder how RPCs are created for satellite images. I read that RPCs are used to correct for systematic error in the image which arise due to Surface relief, camera orientation wrt nadir, camera tilt and geometric errors. And that it ensures uniform scale for all pixels in the image i.e., orthorectification. Now I wonder how these rational polynomial coefficients are derived for a region on earth.
Kindly give me relevant references


texts: what does it mean to store cross products in SAR product data file?

texts:
I'm reading the UAVSAR data format
For definition of two types of file we have:

MLC files (.mlc): calibrated multi-looked cross products, floating
point format, three files8 bytes per pixel, three files 4 bytes per
pixel, little endian.
Ground projected files (.grd): calibrated complex cross products
projected to the ground in simple geographic coordinates (latitude,
longitude). There is a fixed number of looks for each pixel. Floating
point, little endian, 8 or 4 bytes per pixel.

I know that the schematic form of SLC file is like this:

where I and Q are successively the real and imaginary part of the complex signal for each pixel after sampling and analog to digital conversion but if I take cross product as explained here, how can I imagine MLC or GRD format schematically?
I mean which components in the SLC file are multiplied by each other to form the cross product of each pixel?
What does cross product of the pixel mean in a physical sense?


texts: Trying to bind pop-up in leaflet via GeoServer GetFeatureInfo. Is there anything wrong in my code?

texts:
Simple Leaflet Map

<script
    src="js/leaflet.js"> </script>

<script src="http://ajax.googleapis.com/ajax/libs/jquery/1.9.1/jquery.min.js"></script>

<script>

     var osm = L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png', {
    maxZoom: 18,
    attribution: "Data by OpenStreetMap"
    });

var dep = L.tileLayer.wms("http://localhost:8080/geoserver/TestpostGIS/wms", {
    layers: 'TestpostGIS:departments',
    format: 'image/png',
    transparent: true,
isBaselayer:false,

});

var ground = L.tileLayer.wms("http://localhost:8080/geoserver/ucity/wms", {
    layers: ' ucity:ground',
    format: 'image/png',
    transparent: true,

});
var map = new L.Map('map', {
    center: new L.LatLng(15, 0),
    zoom: 2,
    layers: [osm, dep,ground],
    zoomControl: true
    });

var departments={
    "Departments":dep
}

var grounds={
    "Grounds":ground
}

var layerControl = L.control.layers(departments, grounds);
    map.addControl(layerControl);

ms_url="http://localhost:8080/geoserver/TestpostGIS/wms";
map.addEventListener('click', Identify);
function Identify(e)
{
    // set parameters needed for GetFeatureInfo WMS request
    var BBOX = map.getBounds().toBBoxString();
    var WIDTH = map.getSize().x;
    var HEIGHT = map.getSize().y;
    var X = map.layerPointToContainerPoint(e.layerPoint).x;
     var Y = map.layerPointToContainerPoint(e.layerPoint).y;
     // compose the URL for the request
    var URL = ms_url + 'SERVICE=WMS&VERSION=1.1.0&REQUEST=GetFeatureInfo&LAYERS=TestpostGIS:departments&QUERY_LAYERS=TestpostGIS:departments&BBOX='+BBOX+'&FEATURE_COUNT=1&HEIGHT='+HEIGHT+'&WIDTH='+WIDTH+'&INFO_FORMAT=text%2Fhtml&SRS=EPSG%3A4326&X='+X+'&Y='+Y;

    //send the asynchronous HTTP request using jQuery $.ajax
    $.ajax({
        url: URL,
        dataType: "html",
        type: "GET",
        success: function(data)
        {
            var popup = new L.Popup
            ({
                maxWidth: 300
            });

            popup.setContent(data);
            popup.setLatLng(e.latlng);
            map.openPopup(popup);
        }
    });
}

</script>


texts: How to save the files output of pgsql2shp command(shp,shx,dbf,cpg,prf) to zip file. I want all the files in one zip file

texts:
I have a script from which I have to execute the pgsql2shp command and prepare a zip file which includes all the files output of pgsql2shp command:

pgsql2shp.exe -f "C:\locations" -h localhost -u postgres testDB "select * from locations"

Below is the output:

Initializing...
Done (postgis major version: 2).
Output shape: PolyLine
Dumping: X [36 rows].

Files will be created in location given in cmd but unable zip them with the single command.
How can I achieve this? Does anyone have any idea? I tried redirect (>>) option but with no luck.


texts:
I'm trying to resolve an error when accessing a MapServer layer via MapCache.

I'll add excerpt of the relevant files below, but the basic situation is:

When accessing a particular WMS layer via MapCache, there is one small area of the map with missing tiles. The rest of the layer works fine.

Apache log shows the following:

[Tue Jan 14 13:23:22 2014] [error] [client ::1] Premature end of script headers: mapserv

[Tue Jan 14 13:23:22 2014] [error] [client 192.168.56.1] wms request for tileset Aerial returned an unsupported format:\n<!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.0//EN">\n\n<title>500 Internal Server Error</title>\n\n
Internal Server Error
\n
The server encountered an internal error or\nmisconfiguration and was unable to complete\nyour request.
\n
Please contact the server administrator,\n [no address given] and inform them of the time the error occurred,\nand anything you might have done that may have\ncaused the error.
\n
More information about this error may be available\nin the server error log.
\n
\n<address>Apache Server at localhost Port 80</address>\n\n

When building the cache for this layer, no errors are reported.

Even stranger, accessing the WMS layer directly through MapServer works perfectly.

The above error is the only meaningful error I've been able to find.

To debug this I've set DEBUG 5 on both MAP and LAYER properties and tried accessing the 'broken' tile both via MapCache and MapServer. The results are... weird.

I'll post the files below, but note that:

The log for the MapCache request ends prematurely: msDrawRasterLayerGDAL(): red,gree
The source and other attributes referenced in the log are vastly different between the MapServer and MapCache request. Everything about the request is the same (projection, bounding box, etc), except the 'bad' request is passed through MapCache.

Files / URLs

MapCache request (broken)

URL: http://MYHOST/maps?LAYERS=Aerial&BGCOLOR=0x444444&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&FORMAT=image%2Fjpeg&SRS=EPSG%3A27700&BBOX=363840,145762,365600,147522&WIDTH=512&HEIGHT=512
Log: apache.mapcache.log

MapServer request (works)

URL: http://HOST/cgi-bin/mapserv?MAP=/path/to/Aerial.map&LAYERS=Aerial&BGCOLOR=0x444444&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&FORMAT=image%2Fjpeg&SRS=EPSG%3A27700&BBOX=363840,145762,365600,147522&WIDTH=512&HEIGHT=512
Log: apache.mapserver.log

MapCache config excerpt - (Paths / URIs altered)

mapcache.xml

Map file

Aerial.map


texts: Integrating numeric scales into a text or HTML box in QGIS print composer.

texts:
I am trying to embed the scale of a map in the print composer expression builder or HTML frame. Effectively the value of the scale in the map item properties would come through into the text and be used as part of a standard text or HTML. This seems to be related to question:

  Is it possible in QGIS print composer to include a numeric scale in an html box?

In QGIS 2.6 (and earlier versions) there is an option to use the scale as part of an expression, however, it isn't possible to use this function in the Print Composer - I think it is probably to do with the need to stipulate which map item you are interested in the scale of (this doesn't seem to be possible using the scale function in the expression builder).

Can anyone suggest an alternative way of doing this or is what I suggest not possible and one to be captured as an enhancement request?

I am using QGIS 2.6 in Windows 7.


texts: Creating print layout with 3 different maps of same area using QGIS

texts:
I have a problem with the Qgis print composer: I want to create a print layout with 3 different statuses of the same geographical area, stylizing with 3 different fields. So far everything works but when I do the export to PDF or image format all 3 different maps are correct at the last style used. I can not figure out how to solve the problem.

(same problem here: How to avoid map refreshes in print composer? but this solution seems to me it does not work)

I tested on the following systems:

Windows XP + Qgis 1.8 (postgres/postgis data)
Windows 7 + Qgis 1.8 (shape data)
Ubuntu 13.04 + 1.9 Qgis nightly b. (postgres/postgis data)


texts:
When I add an ArcGIS feature service in QGIS, the Digitizing Tools become disabled.

Is there any way to edit data from ArcGIS Rest Services (Feature Services from ArcGIS Server 10.5.1 like WFS services) with QGIS 2.18.20 ?


texts: Leaflet Plugin 'leaflet-icon-layers' doesn't show layer icons in correct size

texts:
I'm using the leaflet-icon-layers plugin to switch between a basic map and a satellite map, but their icon size isn't correct, although the set images are 80x80px. on my map they looks like this. Also, I use Angular.
The control:
const iconLayersControl = new (L as any).Control.IconLayers(
  [
    {
      title: 'Base Map',
      layer: tiles,
      icon: basic,
    },
    {
      title: 'Satellite',
      layer: esriWorldImagery,
      icon: satellite,
    },
  ],
  {
    position: 'bottomleft',
    maxLayersInRow: 2,
  }
);

iconLayersControl.addTo(this.map);


texts: Lat/Lon/Alt to spherical or cartesian coordinates

texts: On my previous question, it was advised I convert lat/lon/altitude to spherical or Cartesian coordinates. I'm not working near the poles and it would be safe to assume a spherical earth. What would be the best way to go about this with the minimal amount of operations? Would it be best to use x/y/z or phi/theta/rho? I'm working on a small microcontroller, with no hardware FPU (software FPU only!) so every cycle counts.


texts: Generating polygonal blocks using growing lines

texts:
We first generate a set of locations (i.e., points(x,y)), randomly. Then we assign for each point a number showing the angle of a line being centred at that point. The criteria for the length of lines is that, the length of lines (all together) grows slowly until a line intersects another line (as shown below). Then that end-point of growing line will be frozen (no further growth) for that line. The final result will be a set of closed areas.

We are looking for algorithms and ideas regarding the implementation of a solution for the mentioned problem.

Update:
Here we add the result of our implementation in Python of the R code developed below by whuber, the accepted answer. His mathematical solution is very nice, such that even its Python edition took only 186s (~3m) for 500 points, Wow! This is awesome. Note that in this solution there is no further requirement for refining the resulting lines.

Final Note:
In order to get lines instead of rays you can replace (in R):

rays <- array(as.vector(cbind(pts,pts+limits*vec)),dim=c(n.points, 2, 2))

with something like (in Python):

pts += limits*vec
n = len(pts)//2
lines = hstack(pts[:n],pts[n:])

Really Final Note
To mention that to make blocks from touching lines (the product of the post) if you are interested in using of DCEL there is a simple (with some bugs) implementation at {here} for Python lovers. Also you may follow {this post} for other methods in generating blocks (polygons) from touching (intersecting) lines.


texts: Measurement errors/differences in qgis/postgis

texts:
I get different results when measuring a object in QGIS and/or Postgis!

To have something exact I did this:

In QGIS:

create a point shp-layer (epsg 25832)
add a point
created a buffer around it with 100 m
made a polygon from layer extents

So I should have a square with 200 m x 200 m right and 40000 qm?
And that is what is shown in the attribute table which is created automatically by the "Polygon from layer extent"-Function!

But:
When I create a new column and measure the area with $area in QGIS the result is: 39892.72715?

When I import the layer to a postgis (keeping the epsg) table and measure the area with ST_area the result is: 40000.0000000116?

The difference between QGIS and postgis becomes much larger when polygons are involved!

Do I miss something?

Well I/You can reproduce the the result by doing this:

Create a .csv file (called it Point.csv) with this content:

nr,x,y
1,437850,5795350

Add this .csv to a new QGIS project (Set to EPSG:25832) as new Layer with "Add Delimited Text Layer ..." again with EPSG:25832. Now there is point layer called Point.

Now I use Vector > Geoprocessing Tools > Buffer(s) ... with this settings:

Input Vector layer: Point
Segments to approximate: 25
Buffer distance: 50
Output shapefile: BufferPoint.shp

Next I create the square: Vector > Research Tools > Polygon from Layer Extent... with this settings:

Input: BufferPoint
Output: SquarePoint

When I now look into the attribute table (the columns and values are created automaticly) everything looks fine:

MINX      MINY    MAXX    MAXY      CNTX     CNTY         AREA    PERIM
437800.0 5795300.0 437900.0 5795400.0 437850.0  5795350.0   10000.0   400.0

HEIGHT WIDTH
100.0   100.0

But when I create new columns within this attribute table with the field calculator ($area and $perimeter) (both "real", length 15, precision 5)
I get this values:

area1         border1
9973.36847    400.14108

And this is what makes me confused?

PS: the Values I get in PostGis (ST_Area, ST_Perimeter) also differ from the expected and right values!


texts: Compare only valid polygon geometries in PostGIS

texts:
I'm trying to run a fuzzy geometry comparison on a large polygon dataset (>2 million rows), but am getting an error like so:

ERROR: ptarray_contains_point called on unclosed ring
SQL state: XX000

The problem lies with a few annoying invalid geometries being present in one of the tables, and this error is raised on the 'ST_DWithin' statement (see code excerpt below):

SELECT
    DISTINCT ON (t1.gid) t1.gid, t1.id, t2.gid, t2.id,
    (CASE WHEN ST_INTERSECTS(t1.geom, ST_MakeValid(t2.geom)) IS false THEN '9999.99'
          ELSE ST_HausdorffDistance(t1.geom, t2.geom)
    END) as hausdorff
FROM
    (SELECT id, gid, geom FROM new_lyr WHERE gid = 12345) as t1,
    (SELECT id, gid, geom FROM old_lyr) as t2
WHERE
    ST_DWithin(t1.geom, t2.geom, 2)
ORDER BY
    t1.gid, hausdorff;

My question is: How can I somehow bypass the error, while still making use of the in-built PostGIS spatial index (i.e. the ST_DWithin statement). I have tried to put ST_MakeValid(t2.geom) into the WHERE clause, but this takes almost forever to run since it tries to make all the geometries in table t2 valid! Is there any other workaround way to resolve this, either in SQL, or would it be better just to clean up the geometry data in those tables and re-run the query?


texts:
I am developing a district editor for cities. The idea is that we divide the city into districts. I am using GoogleMaps, but any other js library is acceptable.

Problem 1: districts are connected polygons, so changing one's border changes it's neighbour. I can not imagine a simple interface to join the polygons, and the ones, that I can think of are hard to implement (like sticky borders idea where borders 'stick' together when are near and 'unstick' when stretched far away.)

Question: Is there some library or an example of an Editor that can join Polygons, and changing one's border changes the neighbour? (Even non web examples would be of use to give me an idea of the user interface.)

Problem 2: I have a dataset based on data from wikimapia, that contains district polygons. Problem is that those polygons do not join by points and somehow overlap or have space between them.

Question: How can I process those polygons to gain ones, that have joined borders. I can use Python scripts or any other if needed instructions are provided. This is a one time operation, no need to implement a graceful solution.


texts: Updating overview map in Layout View of ArcMap

texts:
I am using ArcMap 10.5.

I am currently working on changing the data driven pages for a map that was last updated years ago. The project requires me making new data driven pages covering a wider area. Part of the legend that I imported from the previous layout (which otherwise needs nothing else changed) is an overview of the map showing the prior grid pattern which was used to make the data driven pages years ago. I have created a new grid pattern to encompass the whole area, but the overview in the legend still only has the old layout, but shows the current data driven page I am looking at (red rectangle) despite it not being part of the grid.

How can I update this overview to represent the changes I have made in the grid pattern for the data driven pages?


texts:
I would like to use different generalization tools for polygons in QGIS 3.0.
Especially for aggregation, simplifying and smoothing.
Back using QGIS 2.8 there were the following possibilities:

Vector -> Geometrie Tools --> Simplify
Using the Toolbox which included the GRASS GIS generalization methods
The Generalizer extension --> https://plugins.qgis.org/plugins/generalizer/

In Version 3.0 I only found the first option. The second one is merged into v.generalize, am I right?

The third option seems to be not compatible with the 3.0 version.

Have I missed something and is there a good alternative I could use instead of the options above?


texts: QGIS Points in a Polygon extract name information from Individual polygon

texts: I am new to QGIS and need to know how to extract names and address' from the points in certain polygons.  My main layer has all the points on the maps and then I drew the polygons on the map.  I then counted all the points in the polygons and that has all worked fine but now i need to know the names and address' in the polygons.  It is probably a simple fix but I have tried for a day and cannot figure it out.  It has to be something in the attributes table and joining the information but I am lost.


texts: Creating an NDVI time series in GEE using Landsat 5 & 7

texts:
I have used code from an example by Tim Assal, who used Landsat 8 imagery, but adapted this for my area (the Eastern Cape, South Africa).

Everything works fine but I need to go back as far as possible for my time series (late 1970's) so need to use Landsat 5 and/or 7 in my study area.

However, when I change the code to use Landsat 7 data (see below), it doesn't show the NDVI.

//Import Macubeni boundary
//Add Macubeni boundary to layer

Map.addLayer(boundary,{color: 'green'},'boundary');

// load all Landsat 7 SR image within boundary
var l7collection = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
          .filterBounds(boundary)
          .filterDate('1999-01-11', '2013-03-30');

print(l7collection);
Map.addLayer(l7collection);

// Mask pixels with clouds and cloud shadows using the 'pixel_qa' band
// that is based on the cfmask for surface reflectance.
// https://landsat.usgs.gov/landsat-surface-reflectance-quality-assessment

// create function to mask clouds, cloud shadows, snow
var maskClouds = function(image){
// make a new single band image from the pixel qa band
var pixel_qa = image.select('pixel_qa');
// retain clear (0) and water (1) pixels
return image.updateMask(pixel_qa.eq(322));
};

// use "map" to apply the function to each image in the collection
var l7masked = l7collection.map(maskClouds);
print(l7masked);

// // create function to add NDVI using NIR (B4) and the red band (B3)
var getNDVI = function(img){
return img.addBands(img.normalizedDifference(['B4','B3']).rename('NDVI'));
};

// // map over image collection
var l7ndvi = l7masked.map(getNDVI);

// // for each pixel, select the "best" set of bands from available images
// // based on the maximum NDVI/greenness
var composite = l7ndvi.qualityMosaic('NDVI').clip(boundary);

// //set visualization parameters for NDVI
var ndviParams = {min: 0, max: 1, palette: ['red', 'yellow', 'green']};
Map.addLayer(composite.select('NDVI'), ndviParams, 'NDVI image');


texts: Create MultiPolygon from QgsGeometry.fromMultiPolygon

texts:
I have a list of polygons' points like

[[x1,y1],[x2,y2],[x3,y3],[x4,y4]]   [[x5,y5],[x6,y6],[x7,y7].[x8,y8]]

now I want to create multipolygon from it by using

feature = QgsFeature()
feature.setGeometry(QgsGeometry.fromMultiPolygon([polygon_points]))

but it returns an error:

TypeError: QgsGeometry.fromMultiPolygon(): argument 1 has unexpected type 'list'

So how to make it right type for QgsGeometry.fromMultiPolygon()?

or maybe I am wrong with something?

======================Newly Update======================================

following is my point list

Points_list

My code read this list by

read_polygon=open('shapes_point.txt','rb')
polygons=[]
polygon_lines=read_polygon.readlines()  #read all polygons
for line in polygon_lines:              #read per polygon in polygons
  for points in line:                   #read per point in polygon
        x=points[0]
        y=points[1]
        Qgs_point=QgsPoint(float(x),float(y))
        polygons.append(Qgs_point)
feature = QgsFeature()
feature.setGeometry(QgsGeometry.fromMultiPolygon((QgsGeometry.fromMultiPolygon([polygons for polygon in polygons]))))

and returns me the error above

=====================================================

I will try to make it more details

I want to read a shapefile in a dialog box, after some calculate progress, I will put values in each polygon features.  The following is the main problem I have now!

I have the shapefile, and I want to fill it with color by values, but the method I only know is make a new polygon shapefile and append value to it

I don't know how to append the value in the selected polygon shapefile

So I use

feature = QgsFeature()
feature.setGeometry(QgsGeometry.Polygon([polygon_points]))

but it's always crash in multipolygons

so I want to know how to make multipolygons function need?


texts: Open Topography Coordinate System xmin ymin xmax ymax

texts:
I am trying to work with some GeoTiffs downloaded from OpenTopography.

I need to know the x dimension and y dimension of the region extracted in units, for example meters, but I do not know how to get those from the values OpenTopography provides:

Xmin = -63.080749511718736
 Ymin = -38.666854934630194
 Xmax = -61.00982666015623
 Ymax = -37.28126487387608

Maybe there is a way to convert those to unit dimension using Qgis?


texts: SSL handshake error when downloading images

texts:
I'm learning the Planet API and going through the tutorials posted online.  When I attempt to perform an image search, download, etc. I get an SSL handshake error.

Some example commands I am calling are from here: https://www.planet.com/docs/api-quickstart-examples/cli/

I get the following error:

requests.exceptions.SSLError: HTTPSConnectionPool(host='api.planet.com',
port=443): Max retries exceeded with url: /data/v1/quick-search?
_page_size=100 (Caused by SSLError(SSLError("bad handshake: Error([('SSL
routines', 'tls_process_server_certificate', 'certificate verify
failed')],)",),))

Downgrading my python libraries and turning SSL cert checking haven't worked.  I am using Python 3 on a Windows machine.


texts:
I'm trying to compare two datasets; a satellite image and some known points. In the satellite image I've got 20 to 50 detections that should match up to a set of lat/lon points that I have. There is generally a unknown uniform offset between the two datasets. This problem isn't too difficult when the points are far apart from each other but when they get close I start running into mismatching problems.

So in this dummy pic I've got some satellite image detections (red circles) matched to some ground truth samples (green squares). I generally end up with some false red circles (bad detections, artifacts etc) and some unmatchable green squares (resolution too bad to detect them). There is often a uniform offset that's not predictable (top left corner) that causes some mismatches. This can't be corrected for since it's a function of a bunch of variables and not just a geo-coding error. For example the samples in the dotted ellipse would be matched by a simple threshold algorithm.

What's the correct algorithm to be using here? Some kind of global nearest neighbor matching with the input features being straight up x,y coords? Is there something like this in sci-kit learn or some other python lib?


texts: Adding coordinates that start with 0 in ArcGIS Desktop?

texts:
I have a list of coordinates of archaeological sites I'm working on in the south of Israel and Sinai. The points were taken in the old Israeli system (Palestine 1923 Israel CS grid) in which some of the points start with 0's. When I try to add the table from excel or create a new one in ArcMap it automatically deletes the 0's so the points aren't accurate. The only option I've found to keep the 0's at the beginning is if the field is marked as a text field, in which case it won't allow me to add the points as coordinates.

How can I add coordinates that start with 0?


texts: Calculating indices such as NDVI or NDSI on Landsat 8 imagery whose pixel type is 16 bit signed integer?

texts:
I am trying to calculate an NDSI (Normalized Difference Snow Index) using Landsat 8 OLI level 1 imagery atmospherically corrected in ERDAS ATCOR using raster calculator in ArcMap. NDSI is used to to detect the presence of snow cover and is calculated by using the following equation:

(band 3 - band 6)/(band 3 + band 6)

When I calculate index on the raw imagery I get a suitable index that discriminates snow cover from other land cover and ranges from -1 to 1. However, if I try to do same index on the atmospherically corrected imagery I get a range from 1189 to -1135 and a very weird output image.

I noticed that the imagery in the raw imagery had a pixel type and depth of 16 bit unsigned integer compared to the corrected imagery which had a pixel type and depth of 16 bit signed integer. So I'm assuming this is the case. Somehow in the process of atmospheric correction in ERDAS the imagery was converted from 16 bit unsigned to 16 bit signed.


texts: "Raster map not found" error using GRASS r.walk

texts:
I'm using GRASS' r.walk for the first time, and I'm trying calculate the best path between two different points. I have a DEM (from SRTM) and a slope map derived from that same DEM. Both are in the same directory, and I'm running it via command line.
The command is pretty simple, but I'm getting an error I'm not sure how to troubleshoot.
R.walk elevation=cut_n30e000.tif friction=Slope_Italy_Greece.tif output=pathing start_coordinates=36.440683,28.216723 stop_coordinates=36.195597,27.964125

I get the following from running the command:
Walking costs are a=0.72 b=6 c=1.9998 d=-1.9998

Lambda is 1

Slope_factor is -0.2125

ERROR: Raster map <cut_n30e000.tif> not found

What do I do about that error? How is it "not found"? It's in the same directory.


texts: Changing Polygon shape on Zoom in and Zoom out in google maps using javascript

texts: In my application iam using google maps v3, on page load with default zoom level of 15 i load few polygons by passing points(latlng coordinates) from database on the google map now my question is when i zoom out the map the polygons shape has to be changed and size of the polygons has to be reduced how can i do this using javascript.


texts: Unable to register gdb with the server due to sql server client incompatibility

texts:
There are two virtual machines:

One where I've created an enterprise geodatabase on a Microsft SQL 2014 server- with a Windows Server 2012 R2 64-bit Operating system
And a second machine where ArcGIS Server 10.3 and ArcGIS Desktop 10.3 are installed + Microsoft SQl Server 2012 Native Client 64-bit SP1.
This machine has a Windows Server 2012 R2 64-bit Operating system, as well.

I'm trying to publish a map service on the server with data from the enterprise gdb. When trying to register my geodatabase with the server, I get the following error:

  " The connection property set was missing a required property or the
  property values was unrecognized. Connection was attempted with an
  older version of SQL Server client communications software that is not
  compatible with the SQL Server database server"

So I've search for Microsoft SQl Server 2014 Native Client but couldn't find it. There is no SQL Server 2014 Native Client.

It's the same problem as this one:  https://geonet.esri.com/thread/84553, but for me installing the 64 bit client didn't fix the problem.


texts: OpenLayers FrameCloud does not find JavaScript functions

texts:
I am trying to use a function inside the FrameCloud function:

if (selectedFeatures.length > 0) {
                    popup = new OpenLayers.Popup.FramedCloud(
                        "featurePopup",
                        map.getLonLatFromPixel(e.xy),
                        null,
                        buildSelectTagWithFeatures(),
                        null,
                        true
                    );

                    map.addPopup(popup);
                }

Which then calls buildSelectTagWithFeatures:

function buildSelectTagWithFeatures()
            {

                var record;
                var info = "
";
                info += "<select id='featureSelect'>";

                // Loop through all objects found
                for (var i = 0 ; i < selectedFeatures.length; i++) {

                    // Setup next object
                    record = selectedFeatures[i];

                    // Setup select option
                    info += "<option value='" + i + "'>" + record[0].Value + "</option>";

                }

                info += "</select>";
                info += "<button onclick='buildSelectedFeature();'>Show</button>";
                info += "
";
                info += "

";

                return info;
            }

At this point, the user selects a feature from the drop-down and clicks on "Show", at this point the onclick function is suppose to fire.  The funtion buildSelectedFeature() is there but it cannot find it.  In Firefox, it is giving me an error "ReferenceError: buildSelectedFeature is not defined".

Why does it not find the function and how do I get it to find the function?


texts: Trying to make private land white overlay hillshade?

texts:
I'd like to have anyone's feedback on this draft map I created. Tried to make the private ownership white but since it's RGB is not near to 255.

I've seen someone created this map similiar to it and leaves it as light grey color.

http://www.cartotalk.com/index.php?showtopic=7125

I am using ArcMap 9.3.1. I have Photoshop 7 on my computer but no training background on it. Have use couple of times in the past not related to GIS.


texts:
internet -- router -- reverse proxy server -- arcgis server

Had this config. working great, after following this proxy server setup tutorial:
http://support.esri.com/en/knowledgebase/techarticles/detail/39568

I was able to see all services by visiting my external facing FQDN:
hxxp://mydomain.com/arcgis/rest/services

I could click on the "view in Javascript" link on the services page and was able to preview the service; all was well.

I then installed SSL cert. on proxy server, communication from proxy to gis server is over port 80, not encrypted. I updated the rest.conf file on the gis server and an now can browse: hxxps://mydomain.com/arcgis/rest/services, all the services are listed.

however, the view in javascript does not preview the service. the viewer loads, but the service's content does not load.

Here's the firebug error:
GET MapServer?f=json&dpi=96&transparent=true&format=png8&callback=dojo.io.script.jsonp_dojoIoScript1._jsonpCallback 404 not found

the URL for the above GET request is wrong,

The issue is http://
If it were https:// it would work

In the meantime, I just added a premanent redirect to the proxy server httpd.conf, but this is a work around.

anyone have any advice?


texts: Label only features that intersect a specific line

texts:
hello I have two lines layers in QGIS line1 and line2 with field my_field(line1) and label(line2).

I have a global project variable label_filter=one

now I want to create a base label rule,first I want to filter line1 with project variable label_filter  using something like this my_field like label_filter to show only that labels and after this filter to create addition rule to show only the labels from line2 where have intersect with line1.

i try this but without success :

 geomintersects("line1",'my_field') and dbquery('line1','my_field','my_field like  @label_filter ')

first image :

that I want :


texts: Resampling raster based on more than 16 nearest cells?

texts:
I have a habitat suitability map, with cells of 30x30, each with a score of relative suitability for fox. The problem is that, with 30x30 cells, the grain is too fine for getting a broad sense of what general areas are likely to support the species. Whether or not a region will be suitable will also depend on the suitability of the nearby cells, not just the immediate cell.

In trying to best represent this visually I think I need to do some sort of resampling to get the smoother effect I'm looking for. The issue though, is that the built in methods don't take into account enough of the neighbouring cells to accurate reflect the score. Cubic appears to take the average of the nearest 16 cells into account, but I'd like to incorporate more. I don't actually want to change cell size, just the values within. To that end, I'm not even sure that "resampling" is what I need.

I'm using ArcGIS 10.3


texts:
I am receiving a table with customer names, group of customers and some data.
I need to create dynamically Layers of these groups, then add these customer to these respective groups, then display a default base layer with the possibility to tick also a group layer.
I am good with most of it except "create dynamically Layers from a table of string".
The declaration or creation of Layers is wrong the syntax probably or a missing part.

...

    function hide_marker( xxLayer, xxChecked ) {
        console.info("hide_marker: "+xxLayer);
         if( !xxChecked ){
            mymap.removeLayer(xxLayer);
         }else
         {
            mymap.addLayer(xxLayer);
         }
    }

    var locationsCust = [
                     ["Group1", "Cust1", "Address Cust1", 50.00000000, 10.00000000,"00FF00"]
                    ,["Group1", "Cust2", "Address Cust2", 50.00000000, 11.00000000,"00FF00"]
                    ,["Group2", "Cust3", "Address Cust3", 50.00000000, 12.00000000,"00FF00"]
                    ];

    var legendrows = '';
    var rowidNameGroup;
    var layerAllCust= L.layerGroup();
    var mbUrl  = 'https://api.mapbox.com/styles/v1/{id}/tiles/{z}/{x}/{y}?access_token=pk.eyJ1IjoibWFwYm94IiwiYSI6ImNpejY4NXVycTA2emYycXBndHRqcmZ3N3gifQ.rJcFIG214AriISLbB6B5aw';
        mbAttr = 'Map data &copy; OpenStreetMap contributors, ' +
                    'CC-BY-SA, ' +
                    'Imagery © Mapbox';

    var mapView        = L.tileLayer(mbUrl, {id: 'mapbox/streets-v11', tileSize: 512, zoomOffset: -1, attribution: mbAttr});

    var layerName = '';

    var mymap = L.map('mapid', {
        center: [50, 10],
        zoom: 5,
        layers: [mapView, layerAllCust]
    });

    L.tileLayer(mbUrl, {
        maxZoom: 10,
        attribution: mbAttr,
        id: 'mapbox/streets-v11',
        tileSize: 512,
        zoomOffset: -1
    }).addTo(mymap);

    var overlays = {
        "Cust": layerAllCust
    };

    var baseLayers = {
        "Map": mapView
    };

    **var overlayMaps = {};**

    var legend = L.control({position: 'bottomright'});

    for (var i = 0; i < locationsCust.length; i++)
    {
        var locData = locationsCust[i];
        var Icon = L.IconMaterial.icon({
            markerColor: '#'+locData[5],  // Marker fill color
            outlineColor: 'black',        // Marker outline color
            outlineWidth: 1,              // Marker outline width
        })

        L.marker( [locData[3], locData[4]] , {icon: Icon})
            .bindPopup(locData[1])
            .addTo( layerAllCust);

        if (rowidNameGroup != locData[0])
        {
            rowidNameGroup = locData[0];
            **overlayMaps[layerName] = "POP"+i;**
            //add the layer to the overlayMaps with a dynamic layer name
            legendrows = legendrows + ""
                                    +"<input type=checkbox class='xxVBID' onClick='hide_marker(overlayMaps[layerName], this.checked);' />"
                                    +""+rowidNameGroup+"";
            console.info("Layer :"+overlayMaps[layerName]);
            ****
        }
    }
    ...


texts: Install other python distributions [e.g. Python(x,y)] on top of ArcGIS installation ?

texts:
I'm using ArcGIS 10 and would like to use various additional libraries that are not included in the delivered version of ArcGIS (esp. numpy and matplotlib).
First I tried to install them manually as described here: http://matplotlib.sourceforge.net/users/installing.html

But because this seem not to work

I tried various different mailing lists and forums but so far it seems that this problem cannot be solved; see also my post here: matplotlib does not work in ArcGIS 10 and Eclipse (Python)

I'm wondering if it may work if I just install the distribution PYTHON(X,Y) which has all packages I need for development. (see http://www.pythonxy.com/)

Does anyone have experience and knows if I can just install another distribution simply on top of the existing version installed by ArcGIS 10 without loosing my ArcGIS Python functionalities. Maybe someone did even try this and installed the PYTHON(X,Y) distribution to use with ArcGIS already? What version would be best? The newest one runs with Python 2.6.6.


texts: I am pretty new to QGIS and was hoping some one could point me on how to put 3d spikes on to my 2d map? I have lat, long data and a value (lets say population] for each of the record. I am trying to map the population as a 3d spike.


texts: Change the label of the ui.Chart.image.seriesByRegion in Google Earth Engine

texts:
I would like to change the label of the chart in Google Earth Engine. I had the code:

var image= ee.FeatureCollection('ft:1OXpMKiwvIBC__iYcgjyW23A2oMr89IQlC3cLqk3q')
     .filter(ee.Filter.or(ee.Filter.eq('name', 'Lake Tekapo')));
print (image);
var collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
            .filterBounds(image)
            .filterDate('2017-01-01','2017-05-01');
// Map a function over the Landsat 8 TOA collection to add an NDVI band.
function test(image) {
  var equation = image.expression(
    '(7.27 * B4/B3) - 1.7',
    {
        B3: image.select('B4'),
        B4: image.select('B3'),
    }).rename('TEST').float();
  return image.addBands(equation).set('system:time_start',
image.get('system:time_start'));
}
var test =collection.map(test).select(['TEST']);
print (test,'median');
//CHART TIME SERIES
 // Define customization options.
var options = {
  title: 'TEST over time in regions',
  hAxis: {title: 'Time'},
  vAxis: {title: 'TEST values'},
  lineWidth: 1,
  series: {
    0: {color: '00FF00'},
    1: {color: '0000FF'},
}};
// Create a time series chart.
var tempTimeSeries = ui.Chart.image.seriesByRegion(
    test, image, ee.Reducer.mean(), 'TEST',300)
  .setChartType('LineChart')
      .setOptions(options);
// Display.
print(tempTimeSeries);

I would like that the chart of the lable is the property name of the feature (name =  Lake Tekapo) instead of the Id of the feature (id = 2819) (see the figure attached).


texts: Why does my PyQGIS standalone script not work

texts:
I wrote a pyqgis script that should iterate over a folder with xyz files, check if the bouding box overlaps with a polygon in my QGIS project and if yes, it should copy the file into a new folder.

I don't get an error message but nothing happens either. My output Folder remains empty. I just don't see the problem. Can someone please help me find the error, I've already spent too much time on this.

def moveFiles(self):
    mask = self.getBoundingBox() #BoundingBox of polygon layer
    inputDir = self.dlg.input_folder.text() #directory of xyz files - "Users/denni/Desktop/xyz_data"
    outDir = self.dlg.output_folder.text()  #output directory - "C:/Users/denni/Desktop/output/"
    crs = "EPSG:25832"

    for file in glob.glob("*.xyz"):
        uri = "file:///" + inputDir +"/"+ file + "?type=csv&delimiter=%s&crs=%s&xField=%s&yField=%s" % (" ", crs, "field_1", "field_2")
        name = file.replace('.xyz', '')
        lyr = QgsVectorLayer(uri, name, "delimitedtext")
        lyr.selectAll()
        lyr.boundingBoxOfSelected()
        outputPath = outDir + file
        self.iface.addVectorLayer(file)

        if mask.xMinimum() <= lyr.xMinimum() <= mask.xMaximum() and mask.yMinimum() <= lyr.yMinimum() <= mask.yMaximum():
            QgsVectorFileWriter.writeAsVectorFormat(lyr, outputPath, 'utf-8', lyr.crs(), "CSV", layerOptions='GEOMETRY=AS_XYZ')
        if mask.xMinimum() <= lyr.xMaximum() <= mask.xMaximum() and mask.yMinimum() <= lyr.yMinimum() <= mask.yMaximum():
            QgsVectorFileWriter.writeAsVectorFormat(lyr, outputPath, 'utf-8', lyr.crs(), "CSV", layerOptions='GEOMETRY=AS_XYZ')
        if mask.xMinimum() <= lyr.xMaximum() <= mask.xMaximum() and mask.yMinimum() <= lyr.yMaximum() <= mask.yMaximum():
            QgsVectorFileWriter.writeAsVectorFormat(lyr, outputPath, 'utf-8', lyr.crs(), "CSV", layerOptions='GEOMETRY=AS_XYZ')
        if mask.xMinimum() <= lyr.xMinimum() <= mask.xMaximum() and mask.yMinimum() <= lyr.yMaximum() <= mask.yMaximum():
            QgsVectorFileWriter.writeAsVectorFormat(lyr, outputPath, 'utf-8', lyr.crs(), "CSV", layerOptions='GEOMETRY=AS_XYZ')

I'm using QGIS 3.4


texts: Making polygon from intersecting lines in QGIS using Polygonize or similar?

texts:
In QGIS, I want to make polygons from a set of intersecting lines. Looking at the screenshot here, I want chop the purple line (it's a closed line, not a polygon) where it crosses the green and blue lines. The result hopefully would be a polygon where I have hatched in red. I have quite a few to do.

I have tried the QGIS tool previously called Polygonizer (see An Introduction to the Polygonizer Plug-in), but now called Polygonize. This is a tool that, given a layer containing a set of crossing lines arranged essentially like a # symbol, should create a polygon from the enclosed space in the middle.

(In QGIS 2.01, Polygonize is in the top menu bar: go to Processing -> Commander, then type "Polygonize" on the command line, select "Processing algorithm: Polygonize" from the list that appears, and press Enter.)

At the momemt Polygonize doesn't work properly. Is there another way to do what I want to do?

For example, wouldn't it be great if the QGIS 'Split Feature' command was able to work like MapInfo's 'Polyline Split', where you select the polygons to split, then select an already-drawn polyline to split by, and voila. But Split Feature seems not to work that way.

UPDATE: Detail of problems with running Polygonize.

I have 64-bit QGIS 2.01 installed to Windows 7 from the installer package, QGIS-OSGeo4W-2.0.1-3-Setup-x86_64.exe.

The routine processed an output from my shapefile, but did not create a polygon from the region enclosed by intersections between the blue, green and purple lines. (Note: in the actual input file, all these lines are on the same layer; the screenshot just colours them differently for illustrative purposes). What I got instead was a simple conversion of the purple outline to a polygon, as though the program completely ignored the blue and green lines.

When I node-edited the closed purple line, to open it so that I had a set of open crossing lines like in the # symbol, and then ran Polygonize again, nothing happened. No crash, but no visible polygon in the output file.

When I used Nick Hopton's own Laxton.kml data, I got a crash; see screenshot:


texts: Drawing perfect circle with set radius in meters on map using QGIS?

texts:
Following this istructions on how to draw a circle I can make buffer layer on each point I need.

I am not able to draw perfects circle because the buffer will be elliptic?

The CSR is EPSG:4326 - WGS 84

And I am also not able to set unit misure in meters as I need (300 meters)

Where is my error? I use QGIS 1.8 on windowx 7


texts: Running v.clean tool in QGIS Python Console

texts:
How to use v.clean GRASS processing tool for break in QGIS from Python QGIS Console.I have tried using below command:

processing.runandload('grass:v.clean', layer,\
0,0.100000,'8132391.27018,8134827.49723,3988708.55368,3989168.7299',\
-1,0.0010,'memory','memory')

But it is giving me an error of "configuring your grass, it is not configured".

How can I run this command from the Python console?


texts: Image size limitation in QGIS print layout using HTML code

texts:
I'm trying to make an atlas with QGIS, with one image per page.
My goal is more complex so I have to use HTML widget in order to display the images.
The below code works, but only if the picture is lighter than 2 Mo. If bigger, I got a square with "?" in it.
Code used to display image :
[% '<img width="300px" src="' + "relative/path/to/image.jpg" +  '" />' %]

I've already tried to change width / height, without success.
If I change the resolution of a non-displaying picture, it is shown immediately in QGIS Atlas.
Export entire atlas doesn't change the result...
Has anyone already faced this issue?
I wonder if it is a hardware limitation, a QGIS limitation or a QGIS bug...
Configuration : Win10 + QGIS 3.16.2


texts:
I have a problem integrating WMS and WFS. Basically what i want is when user clicks on WMS then a clicked polygon is added to vector (and selected). The problem is that the WMS is in version 1.3.0 and WFS 1.1.0. OpenLayers automatically flips axis for WMS but when it calculates bbox for wfs request, it does not. I have the following code:

    var map = new OpenLayers.Map({numZoomLevels: 20});

             var wms = new OpenLayers.Layer.WMS(
                 "OpenLayers WMS",
                 "http://sdi.geoportal.gov.pl/wms_dzkat/wmservice.aspx?",
                 {layers: ['Dzialki','Numery_dzialek'], version:'1.3.0'}
             );
             select = new OpenLayers.Layer.Vector("Selection", {styleMap:
                 new OpenLayers.Style(OpenLayers.Feature.Vector.style["select"])
             });

    var protocol = new OpenLayers.Protocol.WFS({
                 version: "1.1.0",
                 //url:  "http://sdi.geoportal.gov.pl/WFS_DzKat/service.svc/post",
                 url:'/GSIP/proxy/dopost/',
                 featureType: "Dzialki",
                 featureNS:'http://www.intergraph.com/geomedia/gml',
                 featurePrefix:'gmgml',
                 geometryName:"GEOMETRIA_DEG"
             });

             control = new OpenLayers.Control.GetFeature({
                 protocol: protocol,
                 box: true,
                 //hover: true,
                 multipleKey: "shiftKey",
                 toggleKey: "ctrlKey"
             });

     control.events.register("featureselected", this, function(e) {
                 select.addFeatures([e.feature]);
             });
             control.events.register("featureunselected", this, function(e) {
                 select.removeFeatures([e.feature]);
             });

     map.addControl(control);
             control.activate();

And generated request after click:

<wfs:GetFeature xmlns:wfs="http://www.opengis.net/wfs" service="WFS" version="1.1.0" maxFeatures="10" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
<wfs:Query typeName="gmgml:Dzialki" srsName="EPSG:4326" xmlns:gmgml="http://www.intergraph.com/geomedia/gml">
<ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
<ogc:BBOX>
<ogc:PropertyName>GEOMETRIA_DEG</ogc:PropertyName>
<gml:Envelope xmlns:gml="http://www.opengis.net/gml" srsName="EPSG:4326">
<gml:lowerCorner>19.143386078712 52.371721634254</gml:lowerCorner>
<gml:upperCorner>19.143439722893 52.371775278435</gml:upperCorner>
</gml:Envelope>
</ogc:BBOX>
</ogc:Filter>
</wfs:Query>
</wfs:GetFeature>

The above request does not work (0 features returned) but when i flip coords so it looks like:

<gml:lowerCorner>52.371721634254 19.143386078712</gml:lowerCorner>
    <gml:upperCorner>52.371775278435 19.143439722893</gml:upperCorner>

it works. How to force OpenLayers to build getFeature request with flipped coords? In WMS Layer there is attribute reverseAxisOrder but can't find it in Vector Layer


texts:
I'm still learning the world of GIS reading this Introduction to PostGIS. So far I understand things but since I have an idea that I plan to later work on, I found something that I didn't quiet grasp. No where in the introduction I see the units used by points. Are they supposed to be longitudes and latitudes or some XY arbitrary points that are known to GIS elites?

My confusion results from my future plan: to make an app that can send location and the server API have to tell the app name of the location and some other fancy things about location. Android for example will require longitudes and latitudes. Is there any tutorial that talks about GPS in relation to PostGIS data?


texts: Mismatch in area calculation using R and ArcMap?

texts:
I am new to spatial statistics using R. In this scenario, I am trying to calculate the polygons area with in a shapefile by usingrgeos::gArea and compared with ArcMap calculate geometry value. both the values are of a lot.

NAME        RG_NAME     _ABBREV     AREA      AREA_cal_ByR
Addison      Texas        TX        3.718831  13719277
Carrollton   Texas        TX        17.061196 62942062
Carrollton   Texas        TX        12.675647 46807613

ArcMap Calculated area is of column name "AREA", which is in Square Miles and here is projected coordinated system of the shapefile CRS arguments:
 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0
+y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs

can someone guide me, where i am doing it wrong.

Here is the Rcode that used for area calculation -

inter_zips_within$AREA_clip <- round(rgeos::gArea(inter_zips_within, byid = T), 2)
inter_zips_within@proj4string


texts: OpenLayers. External ArcGISCache layer has different zoom levels/scales than current map in use

texts:
Current map has 40 zoom levels. Map uses mostly internal layers but for this one special layer an external service is used and that is where ArcGISCache comes in. This ArcGISCache layer has its layer info fetched with OpenLayers.Protocol.Script.createRequest(). When an attemp to draw the layer is made a method called ArcGISCache.getURL asks map for current zoom level and it reports 23. However ArcGISCache.lods property has only 13 elements and error is thrown Cannot read property 'startTileCol' of undefined (used caps in the code to mark the place)

getURL: function(a) {
   ...
   a = this.map.getZoom();
       if (this.lods) {
           if (b = this.lods[this.map.getZoom()], c < B.STARTTILECOL || c >
               b.endTileCol || d < b.startTileRow || d > b.endTileRow)
                return null
        } else {
           ...
        }
...
}

How do I integrate this new external ArcGISCache layer into the exising map?


texts: Increase the resolution of Google Earth Engine images

texts:
I am trying to fetch high quality images from Google Earth Engine. I want to understand the relationship between zoom_level option in:
Map.centerObject(image, zoom_level);

and scale option in below code:

Export.image.toDrive({
      image: img,
      description:'desc'
      scale: 1,
      maxPixels:1e13,
      folder:"japan"
    });

I want to get a high quality image of a specific location with as high zoom as possible. I believe scale 1 (1m per pixel) is the highest resolution possible, but the image is not at all clear for human eye. How can I increase the resolution of the image that I am fetching from Google Earth Engine?


texts: ArcMap keeps crashing whilst using simple C# .NET addin

texts:
I honestly have no clue what is causing this.

I want to have a message appear that tells me which layers are selected or not.

Here is my code:

using System;
using System.Collections.Generic;
using System.Text;
using System.IO;
using ESRI.ArcGIS.ArcMapUI;
using ESRI.ArcGIS.ArcMap;
using ESRI.ArcGIS.Carto;
using ESRI.ArcGIS.Framework;
using System.Windows.Forms;

namespace SelectedLayersTest2
{
    public class Button1 : ESRI.ArcGIS.Desktop.AddIns.Button
    {
        public Button1()
        {
        }

        protected override void OnClick()
        {

            IMxDocument mxdoc = ArcMap.Application.Document as IMxDocument;
            IMap map = mxdoc.FocusMap;

            ILayer layer;
            string message = "";

            for (int i = 0; i < map.LayerCount; i++)
            {
                layer = map.Layer[i];
                if (mxdoc.SelectedLayer == layer)
                {
                    message += layer.Name + "\n";
                }
                else
                {
                    message += layer.Name + " not selected\n";
                }
            }

            Message msgForm = new Message();
            msgForm.lbl.Text = message;
            msgForm.ShowDialog();

            ArcMap.Application.CurrentTool = null;
        }
        protected override void OnUpdate()
        {
            Enabled = ArcMap.Application != null;
        }
    }

}

This seems like such a harmless piece of code. What is causing this?


texts:
if I am merging my raster files or if I am creating a virtual raster, the quality is not as good, as the input rasters on their own. Before:

As a virtual raster:

Does somebody know any tips or tricks to maintain the quality ?


texts:
I am trying to geocode address onto a map in ArcGIS 10.3. The table is in dbf format and contains address, city, state, and zip.

I have attempted going to ArcCatalog and adding an address locator. I choose US address (Dual Ranges) and then select my table as the reference, but I receive a red "x" that says "Not enough reference data tables for this style"

Is there another way to display these addresses on my map? Or is there another way to use my table and use the address locator.


texts:
I have a small software that connects to a database (SQLServer) reads a few values and transfers them to a SDE Oracle instance (11g).

The problem is: I have over half million records that I need to transfer. My code reaches for the records without an issue, but when I try to insert...(my guess is that you guys already got throught this) I get a lovely COM error, screaming at me about Single Threaded and Multi Threaded apartments.

Here is a snippet:

        IFeatureWorkspace workspace = GetWorkspace() as IFeatureWorkspace;
        ITable newTable = workspace.CreateTable("blevers", fields, uid, null, "");

        IRowBuffer buffer = newTable.CreateRowBuffer();
        ICursor cursor = newTable.Insert(false);

        IMultiuserWorkspaceEdit workspaceEdit = workspace as IMultiuserWorkspaceEdit;
        IWorkspaceEdit2 workspaceEdit2 = workspaceEdit as IWorkspaceEdit2;
        workspaceEdit.StartMultiuserEditing(esriMultiuserEditSessionMode.esriMESMNonVersioned);
        workspaceEdit2.StartEditOperation();

        for (int i = 0; i <= stringList.Count - 1; i++)
        {
            buffer.set_Value(0, stringList[i]);
            cursor.InsertRow(buffer);
            Console.WriteLine("Escrevendo linha " + i.ToString() + " de " + stringList.Count);
        }
        cursor.Flush();

        workspaceEdit2.StopEditOperation();
        workspaceEdit2.StopEditing(true);

I usually get the error inside the for loop. How can I insert big big cursors without problems?

EDIT: i've also tried using a SchemaLock, but a COM error explodes in my face when converting the ITable to ISchemaLock. I'm out of ideas!


texts: GeoPandas sjoin returning fewer records than expected

texts:
I used GeoPandas sjoin() to merge two Geodataframes. The first one contains 15kk records of buildings and the second one contains all parts of ortofotomap of my country. I have got layer which contains 11.2kk records but I should get the same like the first layer. I was thinking about that and the problem could be when some building stands on 2 parts of map.
What if building stands on 2 parts of map?
Will the code return 2 records, 1, the newest or nothing?
I would like to get 4 records if buildings stands on 4 parts of map because it is possible. For example building can stand in a corner of every from 4 parts of map.
Maybe I should add some argument?
Code:
import geopandas as gpd

bud = gpd.read_file('/home/gis/budynkipolska/output.shp')
kaf = gpd.read_file('/home/gis/completed_kafs.shp')

kaf_bud_join = gpd.sjoin(bud, kaf)


texts: How to calculate cumulative angles on specific routes?

texts:
I am using Network Analyst extension in ArcGIS 10.4 (Windows 10 Pro 64bit) to calculate accumulative angles from Node 1 to Node 423 on the shortest path as shown in the left side of Figure 1.
Network is built from Voronoi poligons which is converted into center line colored green, so it is actually composed of short segments connected each other.
Length of the segments distributes from 0.00319 to 3.51153 which is used as distance impedance.

Figure 1. Network

Followings are the setting applied in Network Analyst:

Network attributes "Angle" is set as "Unknown" unit. Source values of "Angle" is set as 0 for all sources. Turn penalty is set as "Script" in default values.

VB script is set as below -

|||||||||||||||||||||||||

Pre-logic Script Code:

turnAngle = 0

a = Turn.Angle

If a > 330 And a < 30 Then

turnAngle = 0.001 'Straight

ElseIf a > 30 And a < 150 Then

turnAngle = 1 'Right-turn

ElseIf a >= 210 And a <= 330 Then

turnAngle = 1000 'Left-turn

Else
turnAngle=1000000 'U-turn

End If

|||||||||||||||||||||||||

Value = turnAngle

|||||||||||||||||||||||||

In "OD Cost Matirx", "Length (Meters)" is set as impedance and "Angle" is marked in accumulation tab.
The result accumulative angles (Name = Graphic Pick 1-423) shows 160003000, which means 160 times of U-turn and 3 times of Left-turn (there are 163 vertices between node 1 to 423).

However, as shown in the Figure 1, there is no U-turn observed in the network though Left-turn may be possible on the way from 1 to 423 as center line skews at some sections.

Do you have any idea which can fix this problem?

Reference

I also built another simple network to confirm the calculation, the network composed of 5 nodes which are connected by straight (0.001) -> left-turn (1) -> right-turn (1000) -> right-turn (1000000).  The result shows correct total angle as 1002, 1 time of right-turn and 2 times of left-turn, based on the same setting above.

I suppose short length of original network may cause the problem.


texts: Accessing a specific property/column in FeatureCollection in Google Earth Engine

texts:
I am trying to write a function that when a region is clicked on the map, a dot will show and will highlight the region, and show a property in its FeatureCollection (UHIAdded). Here is my code:
lat.setValue('Lat: ' + coords.lat.toFixed(2));
lon.setValue('Lon: ' + coords.lon.toFixed(2));
Dayall.setValue('Annual daytime UHI: ' + ee.Number(UHIAdded.filterBounds(ee.Geometry.Point(coords.lon, coords.lat)).get('UHIday')));
print(UHIAdded.filterBounds(ee.Geometry.Point(coords.lon, coords.lat)).getArray('Area'));

The result showed up as this:
Lat: 37.84
Lon: -120.81

Annual daytime UHI: ee.Number({ "type": "Invocation", "arguments": {
"object": { "type": "Invocation", "arguments": { "collection": {
"type": "Invocation", "arguments": { "collection": { "type":
"Invocation", "arguments": { "collection": { "type": "Invocation",
"arguments": { "tableId":......)

and the printed value is

Null

Why won't FeatureCollection.get work and how do I access the specific property/column


texts:
I am creating an atlas with multiple maps in QGIS.  Each map is at a scale of roughly 1:12,000 and also includes an overview map at a scale of 1:200,000 in the corner which shows the entire region of interest.

QGIS Composer allows you to add an "Overview" (AKA extent indicator in ArcGIS), which shows the extents of larger scale maps within the overview map.

I am hoping to show the extents of all of the maps in the each overview, and not just the extent current map.

It doesn't seem like this is possible in Composer so I tried doing it manually by creating a shapefile of all of the composer extents for the main map from the project.  Then I realized that there is no tool available to quickly generate rectangular features from extent values.

What method can I use to efficiently generate the overview indicators that I need, or polygons which could serve the same purpose?


texts:
I'm using QGIS 2.6.1 and R version 3.2.2.

I'm struggling to get the syntax on a R Script for creating plots via QGIS processing to work. My understanding is that R syntax in QGIS Processing is slightly different from other R interfaces ... and I specifically can't get facet_wrap of ggplot() to generate panels of histograms for my data.

In essence, I have a vector layer with polygon land cover features that I appended mean temperature values for each feature to. I used the "join attributes by location" function to calculate the average of all mean temperature tiles from the BioClim mean temperature raster dataset. A glimpse of the data is here:

I now want to create a panel of histograms using ggplot() that shows the distribution of mean temperature ("MEANDN") for each of the difference Forest Types ("Forest Typ").

The code for generating a histogram across all forest types works:

##Vector Processing=group
##showplots
##Layer=vector
##X=Field Layer

library("ggplot2")

ggplot() + geom_histogram(aes(Layer[[X]]/10))

However, if I add in a facet_wrap() call to try to generate distinct plots for each forest type to show:

# Example attempt... I added in a "##Y=Field Layer" call,
# just didn't show so here.
ggplot() + geom_histogram(aes(Layer[[X]]/10)) +
facet_wrap(~Layer[[Y]], ncol=3, scales="free")

Does anyone have any clues on how the syntax should look to get facet_wrap to work?


texts: Clipping raster using polyline or polygon using ArcGIS Spatial Analyst?

texts:
I have a large raster created from values that cover the whole United States and Alaska which are both polygon and polyline features, this raster was caluclated using Spatial Analyst IDW Interpolation.  How can I get the raster to clip with the either the polyline of the United States border or the polygon shape of the United States?

I followed the advice of @DerekSwingley and set the cell size environment setting to match the original raster.  It gives it a smooth border but instead of having smooth grades of differentiating colors now i have large blocks of differentiating colors.

Good

Bad

Any ideas on how I can fix this and get the Bad raster to look like the orignial (good) raster but keep the bad rasters border?


texts: Using GRASS GIS 7.3 for Atmospheric Correction of Sentinel-Data

texts:
I read that GRASS GIS version 7.3 will be able to atmospheric correct Sentinel-2 Data. After downloading this version, which is still in development, I imported the Data and put all needed information in a text file. After launching the program i.atcorr with all necessary data, a warning message appears: "WARNING: The sun is not raised"

All information in the text file hopefully are as they needed. Has anyone alreade made experience with atmospheric correction of Sentinel-2 Date in GRASS GIS?


texts:
After clipping I am receiving these errors. I was clipping some satellite images from Sentinel 2 via SCP plugin but QGIS is not processing it. I already posted my problem with Scipy where in QGIS 3.8 is not recognizing the python 3.9 I installed and is still using the python 3.6 which I already removed from my MAC.
I am using QGIS 3.8. I don't really understand that multiprocessing problem.


texts:
The subject that I’m working on is the deployment of optical fiber .
I have a shapefile which represents a country containing homes.
And I have another shapefile representing an infrastructure plan of an operator containing boxes (fiber optical room).

I want to cut this country into many areas respecting these rules:

Every area must contain a maximum of 648 houses.
All houses belonging to an area must be accessible from the boxes through the cable ducts of the infrastructure plan.
We can attach a house not touching the optical fiber line with a house that touches it

The result must be something like that.

I have succeed to achieve the first rule. I get this result.
But I couldn't realize tne second one.