texts: Administrative boundaries topology rules

texts:
I have administrative boundaries at 4 levels, each one is a feature class in a file geodatabse.

Each level should be a split of the layer above.
In other words one polygon from level 2 should be filled with one or more complete polygons from level 3.

I am using ArcGIS Topology to validate. What are the topology rules that I need to select?


texts:
This is the code to highlight and pop up a country area on mouse click:

window.onload = function () {
    map = new OpenLayers.Map('map', {sphericalMercator: true});
    var osm = new OpenLayers.Layer.OSM({sphericalMercator: true});

    var shpLayer = new OpenLayers.Layer.Vector({projection: new OpenLayers.Projection('EPSG:4326')});
    map.addLayers([osm, shpLayer]);
    map.setCenter(new OpenLayers.LonLat(0, 0), 1);

    // Interaction; not needed for initial display.
    selectControl = new OpenLayers.Control.SelectFeature(shpLayer);
    map.addControl(selectControl);
    selectControl.activate();
    shpLayer.events.on({
        'featureselected': onFeatureSelect,
        'featureunselected': onFeatureUnselect
    });

    // load the shapefile
    var theUrl = 'naturalearthdata/cultural/110m-admin-0-countries';
    getOpenLayersFeatures(theUrl, function (fs) {
    // reproject features
    // this is ordinarily done by the format object, but since we're adding features manually we have to do it.
    var fsLen = fs.length;
    var inProj = new OpenLayers.Projection('EPSG:4326');
    var outProj = new OpenLayers.Projection('EPSG:3857');
    for (var i = 0; i < fsLen; i++) {
        fs[i].geometry = fs[i].geometry.transform(inProj, outProj);
    }
    shpLayer.addFeatures(fs);
    });
}

    // Needed only for interaction, not for the display.
    function onPopupClose(evt) {
        // 'this' is the popup.
        var feature = this.feature;
        if (feature.layer) { // The feature is not destroyed
        selectControl.unselect(feature);
        } else { // After "moveend" or "refresh" events on POIs layer all
        //     features have been destroyed by the Strategy.BBOX
        this.destroy();
        }
    }
    function onFeatureSelect(evt) {
        feature = evt.feature;
        var table = '
';
        for (var attr in feature.attributes.values) {
        table += '
' + attr + '	' + feature.attributes.values[attr] + '
';
        }
        table += '
';
        popup = new OpenLayers.Popup.FramedCloud("featurePopup",
                             feature.geometry.getBounds().getCenterLonLat(),
                             new OpenLayers.Size(100,100), table, null, true, onPopupClose);
        feature.popup = popup;
        popup.feature = feature;
        map.addPopup(popup, true);
    }
    function onFeatureUnselect(evt) {
        feature = evt.feature;
        if (feature.popup) {
        popup.feature = null;
        map.removePopup(feature.popup);
        feature.popup.destroy();
        feature.popup = null;
        }
    }

I want to highlight a city area using its latitude and longitude - how can I do this?

For Example: I want to highlight California on window load


texts: Inner workings of GRASS reclassification on QGIS. How does it work with multiband rasters?

texts: I have a 33 band raster layer I need to reclassify. I tried using GRASS module inside QGIS and it worked out fine, except that it outputs a raster with a single band. Which band of the 33 did it reclassified? Does it make some kind of merging to take into account all 33 bands? I tried reading GRASS' reclass manual, but it doesn't mention multiband rasters or any similar situation.


texts: Create mean raster stack from raster stack

texts:
I am trying to create a mean stack image(1 panchromatic layer) from a stacked tiff image(13 panchromatic layers). How do i go about doing this?

My code:

from osgeo import gdal
import numpy as np

#read ndvi_stack layers
ndvi_stack = gdal.Open('NDVI_stack/ndvi_stack.tif')
bands = ndvi_stack.RasterCount

for band in range(1, bands+1):
    data = ndvi_stack.GetRasterBand(band).ReadAsArray().astype('float')
    mean = np.mean(data[data != 0]) #calculate mean without value 0


texts: Photogrametry using e-Foto (Stereo Plotter, Stereo Imagery, Stereoscopy)

texts:
I'm learning Photogrametry with Stereo-view using e-Foto (http://www.efoto.eng.uerj.br/en). At the end, I want to know how to do Photogrametry via Stereoscopy in PC with Anaglyph Glasses.

I could replicate successfully e-Foto's tutorial, however trying to work with my own data like Imagery taken by an UAV, Digital Camera, UTM Ground Control Points, Camera Calibration with AgiSoft Lens, I couldn't get a good Stereo-View to work with.

1. My opinion.

I think this is happening because the camera calibration parameters.

Do the parameters need to be converted?

I also think that is because the orientation-rotation. I dont know what this is.

And finally the Fiducial Marks

My imagery doesn't have fiducial marks, so for this I made manually the marks in the images.

What other software do you recomend for stereo-imagery or stereo-view? (I hope use an Anaglyph 3D Glasses)


texts: Setting layer visibility in QGIS Python API?

texts:
Im having a problem toggling the visibility of map layers in pyqgis (2.4). I have a simple interface, just the map canvas, no legend, and basically a radio button that toggles between two raster layers.

Layers are initialized like this:

lyr = QgsRasterLayer(os.getcwd() + '/data/rasters/dof/dof.vrt', 'Ortophoto')
QgsMapLayerRegistry.instance().addMapLayer(lyr, False)
ml = QgsMapCanvasLayer(lyr)
ml.setVisible(False)

lyr2 = QgsRasterLayer(os.getcwd() + '/data/rasters/chm/chm.vrt', 'Chm')
QgsMapLayerRegistry.instance().addMapLayer(lyr2, False)
ml2 = QgsMapCanvasLayer(lyr2)
ml2.setVisible(True)

When i want to toggle their visibility i do:

ml.setVisible(True)
ml2.setVisible(False)

But unfortunately, the rendering on the map canvas stays the same.

What am i missing?


texts: Dropping all fields exept one with QGIS 3 Model Builder

texts:
Is there a way to remove all but one of the fields in the attribute table (with a pre-known name) using QGIS 3 Model Builder?
The problem is, I don't want to include all of these names in the model builder (besides they might not be known beforehand).


texts: Make background color transparent in MapTiler?

texts:
I have MapTiler 5.0.3 pro, and hundreds of 8-bit per band .tiff files with transparent channel Alpha. After batch conversion:

"C:\Program Files\Maptiler Pro\maptiler.exe" -f png -nodata 0 0 0 -zoom 10 18 -store mbtiles -o "outfiles.mbtiles" input_files.tif

and so on...

For some reason some of the output mbtiles have transparency (when load in QGIS, I see only image, without black and\or bounding box). Some output mbtiles have both white and black background (white bg is much bigger), and some only white. What is that?

The input .tiff files have the same metadata (4 bands, no compression). Documentation to MapTiler says, that if input files doesn't have transparency (e.g. jpeg or 3 band png (png8, png24)), MapTiler adds default White background. But, I have alpha band in input. I also tried make bgcolor black (0 0 0 ) and set nodata value (0 0 0), selected different -f options: png, png32, png24, png8, jpeg the result is the same. I also tried tileserver.php but nothing helps. Nodata parameter works fine in most cases, but the white giant  background - I don't need it.


texts: Why would arcpy ApplySymbologyFromLayer tool not work in ArcGIS Pro Python window?

texts:
I can't seem to make ApplySymbologyFromLayer work in ArcGIS Pro console, or in a scripting tool.

It works in an independant python script run in IDLE, but not when executed inside ArcGIS Pro GUI, neither from a scripting tool or the console.

For a basic "single symbol" symbology, this snippet should do the trick :

aprx = arcpy.mp.ArcGISProject("CURRENT")
map = aprx.listMaps(your_map_name)[0]

map.addDataFromPath(path/to/your/layer)

in_layer = map.listLayers(your_layer_name)[0]
symbology_layer = path/to/your/layer.lyrx
arcpy.management.ApplySymbologyFromLayer(in_layer, symbology_layer, None)

It doesn't matter weither I try to apply the symbology from a lyrx file or from a data layer present inside map. Changing arcpy.management.ApplySymbologyFromLayer to arcpy.ApplySymbologyFromLayer_management doesn't change the output either. Neither does replacing the "None" as a last argument of the function by #, or nothing.

When executed, these code blocks say they worked, but the symbology is unchanged on the map. However, it appears to be changed in the symbology tab.

Here, I have tried to apply a grey, 7 pt symbology to this layer currently in pink. It is still pink in the map, and in the TOC, but if you go in the symbology panel for this layer, you see the symbology it is supposed to have.

It looks like a bug to me for the symbology displayed on the map does not match the one you can see on the symbology panel.

Is there a way to make this command work when executed from within ArcGIS Pro console or from a scripting tool ?

I tried this at first using ArcGIS Pro 2.0. I recently updated to ArcGIS Pro 2.1, but the problem is still there. Also tried it on different computers and for different datas, and data types (points, lines, polygons)


texts: Raster size Out of Range error even if within Capabilities' Lower and Upper Corner

texts:
I'm fairly new at WCS so bear with me if this is a really stupid question, but I can't really understand the problem. I'm trying to retrieve some DEM data from the Italian Ministry of Environment's freely available datasets, which are structured as a Web Coverage Service.

In the capabilities of the dataset I'm interested into (http://wms.pcn.minambiente.it/wcs/dtm_75m?service=WCS&request=GetCapabilities)
it clearly states that the Bounding Box limits are the following:

<ows:WGS84BoundingBox dimensions="2">
<ows:LowerCorner>6.53911455679 35.2154346311</ows:LowerCorner>
<ows:UpperCorner>19.6147808261 47.1364182243</ows:UpperCorner>
</ows:WGS84BoundingBox>

But if I try to send the following GetCoverage request:

http://wms.pcn.minambiente.it/cgi-bin/mapserv.exe?map=/ms_ogc/wcs/dtm_75m.map&service=WCS&version=1.1.2&request=GetCoverage&Identifier=EL.DTM.75M&BoundingBox=10.27,41.80,13.36,43.72,urn:ogc:def:crs:EPSG::4326&format=image/tiff

The server keeps returning me this error:

<ows:ExceptionReport xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"    xmlns:ows="http://www.opengis.net/ows/1.1" version="1.1.2" xml:lang="en-US"  xsi:schemaLocation="http://www.opengis.net/ows/1.1  http://schemas.opengis.net/ows/1.1.0/owsExceptionReport.xsd">
<ows:Exception exceptionCode="width/height" locator="InvalidParameterValue">
<ows:ExceptionText>
 msWCSGetCoverage(): WCS server error. Raster size out of range, width and height of resulting coverage must be no more than MAXSIZE=2048.
</ows:ExceptionText>
</ows:Exception>
</ows:ExceptionReport>

Why is that? Where's the error?

I tried also to swap Lon and Lat duplets or give them in the Cartesian coordinates of the dataset, but to no avail...


texts: CartoDB createLayer function not working

texts:
I'm using cartodb.js on my website and it were working fine. But couple of days ago I found that cartodb layer is not loading properly. When I check the console in Firefox, I found this error

ReferenceError: functionreturn_cdbc_cartodb is not defined

I search on google but there exist no suggestions. Clicking on error show

functionreturn_cdbc_cartodb.uniqueCallbackNameJSON.stringifyt.toJSON({
"layergroupid": "an_id_written_here",
"cdn_url": {
    "http": "ashbu.cartocdn.com",
    "https": "cartocdn-ashbu.global.ssl.fastly.net"
},
"last_updated": "2014-10-29T05:03:33.820Z"
});

This is how I'm adding my layer

 cartodb.createLayer(map, url)
 .addTo(map).on('done', function(layer) {
 def=layer})

Another thing is, I'm using same way to add layers on other pages of my website and this function is working properly on those pages. I'm using google maps for base maps.

Any suggestions would be warmly welcome.


texts: ArcMap - Definition query with subquery using a date returns wrong result

texts:
In ArcMap 10.6, I am using the following definition query on a featureClass to get the features that are underway but the results is not what I expected. Data are stored in an Enterprise geodatabase hosted on PostgreSQL.

numeroprojet IN (SELECT numeroprojet FROM FeatureClass_P WHERE datedebut <= CURRENT_TIMESTAMP AND datefin >= CURRENT_TIMESTAMP)

It seems that the subquery does not select all the features where CURRENT_TIMESTAMP is between datedebut and datefin. This happens only in subqueries.

If I add the same query to the definition expression of the featureClass FeatureClass_P, like below, the result is correct

datedebut <= CURRENT_TIMESTAMP AND datefin >= CURRENT_TIMESTAMP

Here the results from the query above in arcmap. It returns 18 records (those in the red rectangle are missing from the same query in pgadmin)

If I do the same query on FeatureClass_P in pgadmin, I got 14 records (just like the subquery in the first statement)

It seems that the subquery in the first statement returns the same result than pgadmin, which is 4 record less than a direct definition query in arcmap.

At the moment of this post the CURRENT_TIMESTAMP from pgadmin is the following:

Any idea why there is such difference? The correct results are those from direct definition query in arcmap (18 records)


texts: What is ArcGIS for Desktop (and Geodatabase) equivalent to Regions of ArcInfo Workstation Coverages?

texts:
Here's one for those who've been around the block a few times.

ArcInfo Workstation from days of yore had an intelligent data model called Regions or Region sub coverages. With this model one could select one or more polygons and create new virtual polygons from the selected set. Furthermore one could dissolve within the set, creating a result akin to merged polygons, or have disparate polygons, possibly overlapping, that identify themselves as one unit (multipart polygons in shapefile parlance). These virtual polygons are called regions. When the regions are loaded into Arcmap (or ArcEdit/View back then) they behave like an independent shapefile or feature class in terms of querying, symbolizing, etc.

The beauty of this model and what is missing from feature classes and feature datasets in the ArcGIS world, so far as I've been able to determine, is that in regions the geometry is stored only once. Imagine, for example, the continent of North America split into polygons depicting countries, islands, states & territories, counties, municipalities, etc. down to whatever level you feel like stopping at, all stored in a single feature class. (ugh, what a mess) Yet, each level has convenient handles (country, province, whatever) that can be grabbed and used as if they were independent feature classes.

Coverages and Regions have their problems, they weren't left behind lightly, but they weren't all bad. Sooo, what is the ArcGIS answer to ArcInfo Regions?


texts: Creating ISO 19139 xml metadata from spreadsheet?

texts: I have metadata that was collected according to ISO categories and is currently housed in a spreadsheet where each row houses metadata for a GIS dataset.  I am interested in automating the creation of XML metadata from this spreadsheet.  A similar question (Writing ISO 19139 compliant XML metadata from spreadsheet?) was asked a few years ago, but I'm hoping that someone may have found a better way! A solution that uses Python or R scripting to create the XML files would be great!


texts: How to add ol.layer.Vector source dynamically?

texts:
How to set source dynamically to a layer vector in OpenLayers 3

like this

var vectorLayer=  new ol.layer.Vector({
    source: new ol.source.Vector()
});

vectorLayer.setSource(another_source);

Is it possible?


texts:
I would like to use Hammer's projection in QGIS. I am providing the follow proj string:

+proj=hammer +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs

But QGIS complains that it is invalid:

According to the proj documentation all parameters are optional, I would expected it to function like this.

I have QGIS 2.18.20 and proj 4.9.3.


texts: Strange "inFeatures" value for Intersect tool

texts:
I am rather new to Python but I am trying to define a function that will take the output of one of my prior functions (5 buffered stream files) and intersect each unique streamBuffer file with an LULC shapefile. So, I want to create 5 unique intersection shapefiles.

My problem comes when trying to figure out how to input this into a variable that will fit inside the arcpy Intersection_analysis tool.

My code currently looks like this:

    def identHabitat(self):
         '''Creates a polygon shapefile by intersecting the shapefiles
         resulting from calcBuffer, identCover, and identDev. The result
         shows areas within a specific watershed where habitat suitable for
         the desired species could be found.'''
         workspace = arcpy.env.workspace = (sys.argv[7])
         arcpy.env.overwriteOutput = True
         buffList = arcpy.ListFeatureClasses("*buff_{0}*".format(self.abbrv))
         unionFile = arcpy.ListFeatureClasses("*{0}_LULC*".format(self.abbrv))
         inFeatures = [unionFile]
         for fc in buffList:
              outFeatures = "{0}_Habitat".format(self.abbrv)
              inFeatures.extend(fc)
              try:
                  intLayers = arcpy.Intersect_analysis(inFeatures, outFeatures)
                  message = "Habitat successfully identified!"
              except:
                  message = "Habitat could not be identified, please try again."
         return message

I put a breakpoint on the line defining the workspace and watched the "inFeatures" variable. This was its value:

Obviously, this is not what I'm trying to do. I want my variable to have a value of: [[u'nopa_LULC.shp'], fc]

And again, I'm new to Python so I'm sure some of this is not the best coding you've ever seen.


texts:
Is it possible to load the services in https://staging.gcm.telstra.com.au/gisserver/rest/services/Hosted into QGIS.
I have tried a Virtual Layer, Vector Tile and ArcGIS Map/Feature service with no luck.
In https://crowdsupport.telstra.com.au/t5/network-roaming/gis-layer-of-telstra-mobile-coverage/td-p/893241 a poster mentions that QGIS 3.14 can't load the format (PBF) but is it possible in newer versions and if so how?


texts: Seeking outline/boundary (shapefile) of Ganges river basin?

texts:
I'm looking for an outline/boundary (shapefile) of the Ganges river basin.

Where can I find this?


texts: Does Re-Projecting ShapeFile only update content of *.prj file?

texts: If I wish to change the projection of some shapefiles from .. let's say NAD83 (4269) to WGS84 (4326) ... and by using some 3rd party app .... does this mean the .prj file is only updated .. and not the content in the .shp file or any other associated files?


texts: How is the Point on surface tool calculated in QGIS

texts:
I am interested to know how the Point on surface tool works in QGIS 3. I am choosing to use this rather than the Centroids tool which works of the "center of mass" methodology because this creates a point outside the area of the polygon which is not suitable for my needs.

The QGIS documentation does not seem to explain the calculation used but I would like to know so I can test whether it is suitable for my analysis and so I can incorporate an error with my calculations.


texts: Determining average RGB values of raster layer image inside selected polygon

texts:
I have a layer which contains a number of polygons(around 10,000). My work is to check the imagery layer colour inside of that polygon and write a colour in the 'Ceiling_ID' field based on its average rgb values.

I need to populate rgb values only on few selected polygons.

How can I achieve this?

I do have a license for ArcInfo, but not able to run zonal statistics.

Imagery layer is an online wms service.

This attached image can explain my problem a little bit more clearly.

Here I have selected a polygon and by analyzing its image I have to give its 'Ceiling_ID' field a value based on its roof colour. Which is light grey in this example.


texts:
After the Route 360 QGIS plugin stopped working, I developed my own script to handle CSV to GeoJSON travel time polygons.

Up until recently, I had been using this page I created to turn CSV files into travel time polygons outputted as a GeoJSON file. This would just use PHP to make a cURL request and give the response to the user. Right now it's in a debugging mode to output as much info as possible.

I have been in discussion with the developers and they can't seem to get to the bottom of the issue.

Does anyone know of any resources for getting travel time polygons for Targomo (formerly Route 360) outputted as a GeoJSON file instead of just being displayed in the browser? All of the examples on the Targomo website are to output to the Web browser instead of a file.


texts: QGIS - Changing the length of Shortest_Line for CAD-like slope pattern?

texts:
I haven't found a tool/plugin that could do this, but I am almost there. Here is what I want to do:

I want to create a slope like the ones from CAD programms, but I don't have and don't need a DEM. It's just for the visualization of the slope.
What I did so far:

Create 2 line layers for the upper slope (blue) and the lower slope (green)
The Tool "Points along geometry" is creating 2 point layers. 1 layer has points every 6 meters, the other one has points every 6 meters with an offset of 3.
With
shortest_line($geometry, aggregate('testline', 'collect', $geometry))

on both point layers, I am able to create lines from one point to the lower slope line, so it looks like a connection of both slope lines.
But I need every second line to be half as long. How can I do that?
Can I change the length of the "shortest_line"? I tried to multiply with 0.5 and other things, but nothing worked so far.
Is there an alternative function? Should I try something with "make_line", some points and angles?
I am not really good when it comes to coding in Qgis.
So any help is very much apprechiated!!
Thanks a lot !


texts: How to manipulate OpenStreetMap schema to be useful in ArcGIS applications?

texts:
I have a shapefile of OSM data. I plan to use the data for an ArcGIS Engine application, and I was wondering if there are any existing solutions to getting the data's schema into your typical GIS format. What I'm talking about is the "tags" attribute. It contains a comma-separated list of categories and their values, whereas typically you would see these values as separate fields in a feature class.

Are there any tools available that separate the tag values out into separate fields? It would probably be straightforward to create a Python script that separates these out, but I'd like to know if a solution already exists before I create one myself.


texts: Modifying feature attributes of a shapefile in geotools

texts:
I have a shapefile with an attribute I want to compute programmatically in Java for each feature. I have found a lot of documentation on how to read and create shapefiles with geotools (http://docs.geotools.org/stable/userguide/library/data/shape.html
http://osgeo-org.1560.x6.nabble.com/Add-an-attribute-to-all-the-features-of-a-shape-file-td4318371.html), but no documentation on how to modify an attribute.

So far, I am loading the shapefile data into a FeatureCollection, and modifying each feature attribute with:

myFeature.setAttribute("ATTRIBUTE", myAttributeValue);

but these modifications are not stored into the shapefile (only in memory I guess).

Is there a way to make these modification persistent into the shapefile?
Maybe there is a more standard and simple way to update feature attribute with geotools. The only thing I would like to avoid is to re-create a new shapefile.


texts:
The following wikipedia page says that 3d mode in bing maps was discontinued

http://en.wikipedia.org/wiki/Comparison_of_web_map_services

Whats the reason for that?
Is there any official announcement or page regarding this from Microsoft?


texts: Creating Data Dictionary File (.ddf) from a Geodatabase

texts:
At my workplace we use Trimble GeoXT GPS devices or similar with CartoPac installed for GPS feature data collection.  We have not used Trimble Pathfinder for several years.

For a few projects, we have been hiring 3rd party contractors to do our GPS/surveying work.  One company we have hired uses software called BlueSky, and we were able to give them a copy of our SDE geodatabase schema for them to generate their GPS data collection forms to match our data collection standards.  Another company we have hired seems to be using Trimble Pathfinder and requested that we provide them with a .ddf (data dictionary file) so that they could generate the .fxl (feature code library file) needed for their GPS data collection forms.

Each company will be providing deliverables in the the form of a file geodatabase at the end of the projects.

My question is:  Is there a way to generate a data dictionary file (.ddf) or feature code library file (.fxl) from an existing geodatabase in ArcGIS?  The third party company that requested this file indicated that it would be too much work to create these manually based on the geodatabase schema that we sent them.


texts: leaflet: how to show draw control in custom (own) div?

texts:
I'm trying to put a Control into an existing div but I don't really know where or how I can force the map.addControl method to show the control (it is a draw control by the way) within an already existing div on the map. I'm using the leaflet draw plugin by the way.

My html looks something like this:

<div class="tooldiv" ng-controller="ClientState">
   ...
</div>

tooldiv is where the control should be placed.

This is my leaflet config:

var drawnItems = new L.FeatureGroup();
map.addLayer(drawnItems);

var drawControl = new L.Control.Draw({
    position: 'topleft',
    draw: {
        polyline: false,
        polygon: {
            title: 'Draw a sexy polygon!',
            allowIntersection: false,
            drawError: {
                color: '#b00b00',
                timeout: 1000
            },
            shapeOptions: {
                color: '#bada55'
            },
            showArea: true
        },
        circle: false,
        rectangle: false,
        marker: false
    },
    edit: false
});

// Add and remove DrawControl menu when layer is selected/unselected
this.toggle_layer_edit = function(edit_polygon) {
    if (edit_polygon === true) {
        if (draw_control_check === null) {
            draw_control_check = map.addControl(drawControl);
        }
    } else {
        if (draw_control_check !== null) {
            map.removeControl(drawControl);
            draw_control_check = null;
        }
    }
}

While searching for an answer I got the idea that it might not even be possible?


texts:
When I began to start the OTP server using the command

java -Xmx2G -jar otp-0.19.0-shaded.jar --build /home/username/otp --inMemory

The program executes normally and says: Grizzly server running.

But, when accessing
    http://localhost:8080/
it gives me "Internal Server Error"

Jar file was not found
 1: org.glassfish.grizzly.http.server.CLStaticHttpHandler.getJarFile(CLStaticHttpHandler.java:349)
 2: org.glassfish.grizzly.http.server.CLStaticHttpHandler.handle(CLStaticHttpHandler.java:288)
 3: org.glassfish.grizzly.http.server.StaticHttpHandlerBase.service(StaticHttpHandlerBase.java:188)
 4: org.glassfish.grizzly.http.server.HttpHandler$1.run(HttpHandler.java:219)
 5: org.glassfish.grizzly.threadpool.AbstractThreadPool$Worker.doWork(AbstractThreadPool.java:565)
 6: org.glassfish.grizzly.threadpool.AbstractThreadPool$Worker.run(AbstractThreadPool.java:545)
 7: java.lang.Thread.run(Unknown Source)

Any Help with this error?


texts: Unable to select layers for intersect in QGIS

texts:
I have a series of buffers around points that I want to intersect with FEMA's National Flood Hazard Layer. My goal is to find out which of the buffered points fall within 1% and 0.2% flood chance zones.
Unfortunately, when I open the Intersection tool, I cannot select the Flood Hazard Zones layer as an input -- it just isn't showing up. This layer comes from an ArcGIS REST connection to FEMA's NFHL server that provides the Flood Hazard Zones. Here are the Properties of the layer:

And here is the screenshot showing the Intersection dialog and the dropdown that doesn't show the Flood Hazard Zones layer:

Why can't I select Flood Hazard Zones as a layer to intersect on? I have a suspicion that there's something about this layer that makes it different from the other layers I'm working with, but what?
I've already tried: checking that the projections are all the same, converting the layer from raster to vector (got a nondescript error), and closing/reopening everything.


texts: GeoServer+PostGIS performance for very detailed data

texts:
I have a layer on GeoServer, taking very detailed data from a PostGIS table.
I can see the size of the table is quite big
Size 977MB, External Size 977Mb
becoming the performance really low. GIST indexes and normal index to unique id have been applied.
As the user zooms in, the speed becomes worse (normal, more detailed data appears)
I have implemented JVM optimization
https://docs.geoserver.org/stable/en/user/production/container.html#optimize-your-jvm
as well as postGIS
https://postgis.net/workshops/postgis-intro/tuning.html
I can cache the initial layer view but then the user has the possibility to apply feature filtering, so at some poing caching is useless.
The geometry types are MULTIPOLYGONS: that is, very few rows on the table but each of these with huge data covering a lot of area.
I was thinking maybe splitting these multipolygons to POLYGONS could help on the performance? so each time user zooms in, GeoServer will find faster which polygons are in the view port and has not to deal with so big (MULTI)polygons.
Another possibility could be to at some zoom level show a simplified layer/table, and 'closer' we are, switch to a more detailed layer/table.
Does some of these options make sense? other advises to take into account?
Moreover, my users are supposed to be on areas with low internet connectivity..
UPDATE:
ST_Subdivide() helps a lot on performance! but newly created polygons border cannot be hidden from Geoserver (even setting polygon border to 0 or transparent), they are slightly visible


texts: Spatially constrained multivariate clustering in ArcPy gives AttributeError?

texts:
i want to use the Spatially constrained multivariate clustering tool for a point shapefile. I'm using PyScripter. The shapefile is named "LITO2016_area_modelli_20190312_v2" and the attribute fields i want to use in the analysis are "Rvmedian", "Rv_int_ran" and "zero". Prova is the output field. Arcpy give me this error:

  AttributeError: 'module' object has no attribute
  'SpatiallyConstrainedMultivariateClustering_stats'

I don't know where i'm wrong.

import arcpy
    from arcpy import env
    from arcpy.sa import *
    #Check out the ArcGIS Spatial Analyst extension license
    arcpy.CheckOutExtension("spatial")
    arcpy.env.workspace = r'C:\python'
    arcpy.SpatiallyConstrainedMultivariateClustering_stats("LITO2016_area_modelli_20190312_v2", "prova", None, "Rvmedian", "Rv_int_ran", "zero", None, None, None,"CONTIGUITY_EDGES_CORNERS")

Can anybody help me?


texts: QGIS 3.18.2: ModuleNotFoundError: No module named 'supy'

texts:
starting QGIS I get the error:

Konnte Erweiterung 'processing_umep' nicht laden aufgrund eines
Fehlers beim Aufruf der classFactory() Methode
ModuleNotFoundError: No module named 'supy'
Traceback (most recent call last):   File
"/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py",
line 335, in _startPlugin
plugins[packageName] = package.classFactory(iface)
:
:
etc

But when I try to install supy

sudo -H pip install supy

I am informed that supy is already installed and that all "Requirement already satisfied".
What can I do?

Python-Version: 3.8.7 (default, Feb 10 2021, 09:04:08) [Clang 12.0.0 (clang-1200.0.32.29)]
QGIS-Version: 3.18.2-Zürich Zürich, 44d7dfabbe
OS: macOS Big Sur 11.3.1


texts: Calculating expected number of points within polygons of different size under hypothesis of random distribution?

texts:
Let's imagine we have some polygons that cover a whole study plot (see attached picture). Further, we have some points (red dots) lying within the polygons.
For each polygon we want to calculate the expected number of points, against which compare the observed number of points, and to calculate the significance (i.e., p-value) of the observed number of points under the Null Hypothesis of a random distribution of points within the study area as a whole.

We actually know (a) the area of each polygon, (b) the total area of the study plot (=sum of the polygons' area), (c) the total number of points, and (d) the number of points falling inside each polygon.

For the p-value, my best guess is to use the binomial distribution, but I can't figure out how to exactly do that. In particular, where I am stuck is in working out the value of p, which in the case of polygons of equal size (i.e., quadrats) should be the reciprocal of the number of quadrats into which the study area is divided. But I do not know how to calculate p in case of polygons of unequal size.

In case of quadrats (equal sized polygons), I have read that we should use a binomial distribution with:

p=1/x (where x is the number of quadrats)

n=number of events in the pattern (i.e., total number of points)

k=number of events in a quadrat.

In the specific case I described, would p be equal to the fraction of each polygon's area relative to the whole study plot (i.e., p=polygon area/sum of all the polygons area)?


texts: Setting X/Y shadow offset values to 1 and -1 in advanced text window using ArcPy?

texts:
I am trying to use a label expression with ArcGIS formatting tags to create a shadow on my text. I have done everything I need except assigning a value of 1 to X offset, and a value of -1 to Y offset.

Below is my code:

import arcpy
mxd = arcpy.mapping.MapDocument(r"CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "Layers") [0]
updateLayer = arcpy.mapping.ListLayers(mxd, "OCParcelsAPN selection", df) [0]
sourceLayer = arcpy.mapping.Layer(r"L:\GISLibrary\Scripts\Python\ReferenceLayers\OCParcels_selection.lyr")
arcpy.mapping.UpdateLayer(df, updateLayer, sourceLayer, True)

del sourceLayer

mxd = arcpy.mapping.MapDocument(r"CURRENT")
layer = arcpy.mapping.ListLayers(mxd) [0]

if layer.supports("LABELCLASSES"):
    for lblclass in layer.labelClasses:
        lblclass.className = "ASSESSMENT"
        lblclass.expression = '"<FNT name=""Futura Md BT"" size=""6""> <CLR red=""255"" green=""255"" blue=""0"">" & [ASSESSMENT] & "</CLR> </FNT>"'
        lblclass.showClassLabels = True

        layer.showLabels = True

        arcpy.RefreshActiveView()


texts: Why don't my polygons from MSSQL2008 display in QGIS?

texts:
I have imported a polygon shape file in MSSQL2008 using Shape2SQL loader without any issue. When I import the loaded polygon file in QGIS, I can't see it ! When I import other point or line features (previously loaded in MSSQL2008), I can see a few polygons from my file but not all of them.

At first I thought that the geometry of the poly file was not complete, so I repaired the geometry using ArcGIS 10.1 and reloaded the ""fixed geometry" file in MSSQL2008, and still am running in the same problem when I import the layer in QGIS.

One thing I have noticed in QGIS: when I connect to my MSQL2008 database, the points and line features have the geometry and SRID already determined but the polygon file requires the user to define the geometry type and SRID. I have tried selecting polygon and multipolygon type and set the SRID to 0 (like the other features in the DB) and still am running into problems viewing my poly file in QGIS. My original poly file has 1285 polygons in it.

I have tried using a subset of my original file and still am running into problems in QGIS.


texts:
I would like to know why I am getting no data pixels in my soil wetness raster, since none of the source rasters used to calculate the SWI have no data values.

method I used: soil data from HWSD FAO used to derive depth and calculated soil hydraulic conductivty by SoilPAR program.
First made TWI: SRTM DEM-watershed clip-flow direction-flow accumulation
DEM-slope-slope radian-remove zero values from radian raster
then calculate TI with raster calculator: this gives a nice raster without No Data values:

For the STI I calculated the transmissivity (=soil hydraulic conductivity*soil dpeth*0.000864), again without No DATA values
Then I used raster calculator for
Ln[(average of transmissivity)(flowaccu)/(transmissivityslope)]
And suddenly everything is pixelated and there are No Data pixels everywhere. I understand that the lakes are NoData pixels but I don't think this should be the same for other places.

Why do I get them and how to prevent them?


texts: True GIS clip of polygons + attribute table

texts:
How can a true intersection/clipping of polygons be done in R?

Problem description:
Shapefile1 contains only one feautre without attributes (or - at least - attibutes don't matter).
Shapefile2 contains several features, some of them fall entirely, some of the partly in the feature contained in Shapefile1. Other features fall totally outside.

I need to perform a true clip or intersection (just like in QGIS Vector > Geoprocessing Tools > Intersect). So far, the solution is described here. However, in addition to the cited I also DO need the clip to inherit attributes of Shapefile2.

At the end, the result needs to be exported to an ESRI Shapefile.

Any answer, suggestion is welcome!


texts: Trouble in Natural Earth Shapefile "jumps" read using GeoTools

texts:
I'm using GeoTools 12.0 to read in Natural Earth data.

I'm really happy with the performance, and ease of implementation.

But, I'm having trouble reading in very large polygon data.  The 110M Ocean dataset contains two polygons: one for the Caspian Sea, and the other for the rest of the World (ROW).

But, in the ROW polygon I'm getting lots of "jumps" in the data.  For example, record 3665 is for this location in North America: (-93.61275590694041, 74.97999726022442, NaN).

The next record, 3666 is for this location in the UK: (-3.0050048486351955, 58.635000108466244, NaN)

Does anybody know if the data includes some other specifier that indicates that the polygon breaks between these two points?


texts:
I am using OSM for my project, I have downloaded data from http://download.geofabrik.de/ in the OSM file and imported the data in POSTGIS, In the planet_osm_roads, I am not able to see any key related to lanes.

Can you please help me get the information about the number of lanes from OSM?

I have visited Getting lane data from OSM

Number of lanes in road data


texts: Clean way to identify color scheme in Python QGIS

texts:
When I do in the Python QGIS console QgsColorSchemeRegistry.instance().schemes(), it returns the following

[<qgis._core.QgsColorScheme object at 0x7f6cc5d8ccc8>,
 <qgis._core.QgsColorScheme object at 0x7f6cc5d8c938>,
 <qgis._core.QgsColorScheme object at 0x7f6cc5f54c30>]

If I loop the schemes with [scheme.schemeName() for scheme in QgsColorSchemeRegistry.instance().schemes()], I get:

[u'Recent colors',
 u'Standard colors',
 u'Project colors']

So, it seems that all schemes have been through a cast so the only way to make the difference is via their schemeName but it's not generic when you use other languages.

In my understanding, I would expect a list with returned content like below:

[<qgis._core.QgsRecentColorScheme object at 0x7f6cc5d8ccc8>,
 <qgis._core.QgsGplColorScheme object at 0x7f6cc5d8c938>,
 <qgis._core.QgsProjectColorScheme object at 0x7f6cc5f54c30>]

First, my guess is right?

Secondly, what is the "best" way to manage this case properly, if any?


texts: dbf to shp within BASH (ogr2ogr would be ideal)

texts:
I have a CSV with many points and attributes that I'm trying to convert into a shapefile but I'm having no luck. I have installed gdal-bin and can use my terminal environment to create a shapefile with ogr2ogr, but I'm not getting the desired outcome. I'm running the following:

ogr2ogr -f "ESRI Shapefile" 5_2834k2-25_FAST.shp 5_2834K2-25_DataCompiled_K2.csv -progress -skipfailures -a_srs EPSG:4326 -nlt POINT

And I get:

"5_2834k2-25_FAST.prj" 1 KB
"5_2834k2-25_FAST.shp" 8 KB
"5_2834k2-25_FAST.shx" 6 KB
"5_2834k2-25_FAST.dbf" 7,625 KB

but despite having all my objects and all the features, my SHP has no geometry. It's extent is 0 and recalculating gives me NaN errors in Arc.
Now my source CSV does have a lat and lon column (two separate columns) and I've tried to glean from other answers how to get ogr2ogr to see and interpret these values but I cannot.

For example, perhaps using this bit of sqlite to choose the columns:

-dialect sqlite -sql "SELECT easting_col as long, northing_col as lat, geometry FROM input"

and "SELECT lon as long [...]" from input where input is my CSV... but nope. I tried a two-step cheat, thinking I could provide the dbf or shp from the incomplete output of the first as that input in the second but no...

I've tried as well GEOMETRY=AS_XY and I'm reading the man pages, the --long-usage help file as well as the GDAL ogr2ogr page but I just can't come up with the proper syntax.

How can I do this? I'm very frustrated because my columns don't have elaborate names ('Lat' and 'Lon') so I figured it would be very straightforward.


texts: Copying rasters which are intersecting with polygons using ArcPy?

texts:
I'm new to ArcPy and to Python in general.

I have a huge raster dataset of all over the World (~1000 tiles .tif) and i'd like to automatically copy the tiles which are intersecting with 24 polygons (one shapefile containing polygons for 24 cities in Africa) in order to keep only the tiles that partially or completely contains my cities. I try using the script below. But when I do, he gives me all the tiles that are in between those cities (so pretty much all tiles in Africa). Probably because he uses the total extent of my shapefile, and not the polygons that are in that shapefile.

Does anyone knows how I could correct this so I get as result only the tiles that are intersecting with the 24 polygons inside that shapefile?

import os, sys, arcpy

AOI     =
"C:/Users/jbihin/Desktop/testGufPython/To_Jeremie/Cities24_50km.shp"  # the
area of interest to find the rasters for
rFolder = "D:/GUF/GUF04" # the folder containing the rasters to search from
oFolder = "C:/Users/jbihin/Desktop/testGufPython/To_Jeremie/gufSelect" # the
folder to copy to

desc = arcpy.Describe(AOI) # get the extent of the AOI
sExt = desc.extent

arcpy.env.workspace = rFolder

for ThisRas in arcpy.ListRasters():
    rDesc = arcpy.Describe(ThisRas)
    rExt  = rDesc.extent
# check if this extent is related spatially
# by using not disjoint
    if sExt.disjoint(rExt):
        arcpy.AddMessage("Raster %s is outside" % (ThisRas))
else:
    arcpy.AddMessage("Raster %s overlaps" % (ThisRas))
    outFile = os.path.join(oFolder,ThisRas)
    arcpy.Copy_management(os.path.join(rFolder,ThisRas),outFile)


texts: Transferring raster band features to another raster in QGIS

texts:
I have one raster layer containing a band that indicates the land type of a random sample of pixels and another raster containing NDVI, VV, and VH values. The first raster has fewer pixels than the second, hence there are areas where the single bands of this layer would show "no data".
How would I go about transferring band values from pixels from the second raster to the first where there are actually pixels in the first raster (the single band has a value)?
A similar question to mine, Extracting value from one raster to another raster using QGIS, suggested using the raster calculator with the expression "layer1"*"layer2", but this only returns one band that seems to have just multiplied the two bands from the input rasters. I am a beginner when it comes to using QGIS, and especially to using the raster calculator, hence I am having trouble finding a solution to this.


texts: ArcGIS Raster Calculator output raster has one band but cannot build attribute table?

texts:
I built a raster using the raster calculator of ArcGIS 10.1 for Desktop.

The properties of the output raster ("usle_1") indicate that this raster has one band.

There is no attribute table but I want to reclassify the values in the raster, so I open the "Build attribute table" tool.

When I input "usle_1", it gives me an error:

  Only single band integer raster dataset is a valid input.

I tried saving the raster as a .tif, to no avail.

Any ideas for how I can reclassify this raster in order to use it in another raster calculator equation?


texts: Assigning latitude and longitude to a particular pixel in satellite image

texts: Pixel of a satellite images can be assigned latitude and longitude, but how this is done as each pixel represent certain area depending on the resolution of image that can be too high to assign latitude and longitude?


texts:
We are having a major issue. We are using ArcGIS Enterprise for our workgroup projects for multi-user editing. Long story short we had a project on a service and this morning it was gone. No trace of it whatsoever. I can't even find a file with the corresponding service name.

I'm trying to figure out what happened. Is there a way to check a log of activity? Such as seeing what happened to different services, changes, publishes, deleted services and so on?

This is a major problem in our office and we need to find out what happened and our periodic back-up does not cut it.


texts: Cannot .saveACopy() When Certain Featureclass in Original MXD

texts:
I'm running a script (ArcMap 10.5) that loops through a list of mxds, and saves a copy with a new name showing the year and quarter. This works for all mxds, unless they have an annotation feature class in them.

When there is an annotation featureclass, I get this error:

  AttributeError: MapDocObject: Unable to save.  Check to make sure you
  have write access to the specified file and that there is enough space
  on the storage device to hold your document.

The script is very simple.

mapList = [r"C:\Projects\Figure 1.mxd", r"Z:\Projects\Figure 2.mxd"]

#make updates
for mapDoc in mapList:
    mxd = arcpy.mapping.MapDocument(mapDoc)

    #generate year and quarter
    some code...
    Q_Y = "Q" + quarter + " " + year

    #generate names/paths
    bn = os.path.basename(mapDoc).replace(".mxd", "") + " " + Q_Y
    #path
    fullPathMXD = os.path.join(C:\Projects, bn + ".mxd")

    #save as new map doc
    mxd.saveACopy(fullPathMXD)

    #create new map doc object to export EPS from
    mxd2 = arcpy.mapping.MapDocument(fullPathMXD)

    #update title/eps name
    mxd2.title = ....

    #export to eps
    arcpy.mapping.ExportToEPS(mxd2, locationPDF + "\\" + bn + ".eps",
                              "Page_Layout",640,480,200,"BEST","RGB",3,
                              "ADAPTIVE","RASTERIZE_BITMAP",True,False)

Has anyone encountered this problem, or can anyone tell me what's going on here?

The mxds that do not throw this error do not have annotation featureclasses, but do have featureclasses that live in the same .gdb.

UPDATE

I used disutils to copy the mxd to a new location, under a different name, and this solved the first error. However, it now throws this error:

  AttributeError: 'NoneType' object has no attribute 'exportToEPS'

Again, this only happens when there are annotation layers in the map. Without them, the original, and the "disutils", version work as expected.


texts: QGIS: Use 'replace' expression for 2 characters in one field calculation

texts:
I have a field of currency values in an attribute table. As an example, I'm attempting to replace £10,000.00 to 10000.00 so QGIS can run graduated styles on the new field. The following replaces the £ sign with a blank character.

replace(  "field_name", '£', '')

I'm struggling to get rid of the comma in the same expression. I have tried:

replace(  "field_name", '£', '') || replace(  "field_name", ',', '')

but this returns 10,000.00£10000.00

The 'replace' expression doesn't seem to cater for more than one character to target as a replacement. Do I need to do this over 2 expressions?


texts: Adding Excel file with degrees to map with CRS in QGIS?

texts:
How can I add a layer with X,Y in degrees (from an excelfile) to a map with CRS?

I am using QGIS 2.8.2


texts: Enabling QGIS Official Plugin Repository in the 3.4.1 Version?

texts:
I am using QGIS version 3.4.1.
Recently, I cannot install new plugins. I have looked in the plugin setting. It seems that the link is blocked. The address is "https://plugins.qgis.org/plugins/plugins.xml" and it is directed to the "QGIS Python Plugins" website.

In the log Messages I have this warning "WARNING: Failed to load C:/PROGRA~1/QGIS3~1.4/apps/qgis/plugins/grassplugin7.dll (Reason: Cannot load library C".

I have a proxy from my university, disabling the automatic configuration script did not solve the problem.

Also, re-installing QGIS did not help.
What else should I look into in order to enable the plugin link?


texts: Copy a .prj-file from another shapefile to a new one when the projection is known?

texts:
I have a bunch of shapefiles of which I know the projection (EPSG:28992). However, they did not come with a prj-file (.prj).

Rather than defining the projection and re-saving the files, I wonder whether it would be okay to copy (and rename) the prj-file from another shapefile which is in the same projection?


texts: Export properties of a featureCollection along with features of it in Google Earth Engine

texts:
I want to read an imageCollection, reduce all images in the collection by the region. and save the result.

However, when I save the result, the properties of the feature Collection is gone. It does not get saved.

Any suggestion?

// Load four 2012 NAIP quarter quads, different locations.
var naip2012 = ee.ImageCollection('USDA/NAIP/DOQQ')
                .filterBounds(ee.Geometry.Rectangle(-71.17965, 42.35125, -71.08824, 42.40584))
                .filterDate('2012-01-01', '2012-12-31');

var reduced = naip2012.map(function(image){
  return image.reduceRegions({
    collection:ee.Feature(ee.Geometry.Rectangle(-71.17965, 42.35125, -71.08824, 42.40584)),
    reducer:ee.Reducer.mean(),
    scale: 30
  });
});

reduced = reduced.flatten();
print(reduced);

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


texts: Converting a PostGIS raster back to vector

texts:
Is it possible to convert a PostGIS raster back to a line?

I'd like to generalize some vector data by converting it to raster and altering the cellsize to fix 'gaps' in the vector connectivity.

Once I've reached a good raster cellsize, I'd like to convert that raster back to a vector line, connecting all the cell centroids into continuous lines within a given tolerance (so as not to connect too many adjacent cells at intersections).

I figure the query would look something like this, but can't figure out what the outer function would be:

Convert To Vector (line) (
  Convert to Raster(cellsize)
  (
  Original Line (geom)
  )
)

One simple example would be to convert the raster 'line' back to a vector as in the ST_AsRaster documentation SQL:

ST_AsRaster(
        ST_Buffer(
            ST_GeomFromText('LINESTRING(50 50,150 150,150 50)'), 10,'join=bevel'),
            200,200,ARRAY['8BUI', '8BUI', '8BUI'], ARRAY[118,154,118], ARRAY[0,0,0]);


texts:
I am attemtping to call a dojo tab container within an ESRI pop-up. The code that I have (see below) is able to create a pop-up with a tab container and the contents, but also creates three 'dead' dropdown menus. Is there a way to create the pop-up without the dropdown menus?

  var map;
  var resizeTimer;
  var identifyTask,identifyParams;

  function init() {
    //setup the popup window
    var popup = new esri.dijit.Popup({
      fillSymbol: new esri.symbol.SimpleFillSymbol(esri.symbol.SimpleFillSymbol.STYLE_SOLID, new esri.symbol.SimpleLineSymbol(esri.symbol.SimpleLineSymbol.STYLE_SOLID, new dojo.Color([255,0,0]), 2), new dojo.Color([255,255,0,0.25]))
    }, dojo.create("div"));

    map = new esri.Map("map",{
      infoWindow:popup,
    });

    dojo.connect(map,"onLoad",mapReady);

    //Add the imagery layer to the map. View the ArcGIS Online site for services http://arcgisonline/home/search.html?t=content&f=typekeywords:service
    var basemap = new esri.layers.ArcGISTiledMapServiceLayer("http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer");
    map.addLayer(basemap);

    var landBaseLayer = new esri.layers.ArcGISDynamicMapServiceLayer("http://<someServer>/ArcGIS/rest/services/test/MapServer");
    map.addLayer(landBaseLayer)
  }

      function mapReady(map){

       dojo.connect(map,"onClick",executeIdentifyTask);

       //create identify tasks and setup parameters
       identifyTask = new esri.tasks.IdentifyTask("http://<someServer>/ArcGIS/rest/services/test/MapServer");

       identifyParams = new esri.tasks.IdentifyParameters();
       identifyParams.tolerance = 3;
       identifyParams.returnGeometry = true;
       identifyParams.layerIds = [0];
       identifyParams.layerOption = esri.tasks.IdentifyParameters.LAYER_OPTION_ALL;
       identifyParams.width  = map.width;
       identifyParams.height = map.height;

       //resize the map when the browser resizes
       dojo.connect(dijit.byId('map'), 'resize', map,map.resize);

      }

       function executeIdentifyTask(evt) {
            identifyParams.geometry = evt.mapPoint;
            identifyParams.mapExtent = map.extent;
    var deferred = identifyTask.execute(identifyParams);

    deferred.addCallback(function(response) {
      // response is an array of identify result objects
      // Let's return an array of features.
      return dojo.map(response, function(result) {
        var feature = result.feature;
        feature.attributes.layerName = result.layerName;
        console.log(feature.attributes.OBJECTID);
        var template = new esri.InfoTemplate();
        template.setTitle("Hello World");
        template.setContent(getWindowContent);
        feature.setInfoTemplate(template);
        return feature;
      });
    });

    // InfoWindow expects an array of features from each deferred
    // object that you pass. If the response from the task execution
    // above is not an array of features, then you need to add a callback
    // like the one above to post-process the response and return an
    // array of features.
    map.infoWindow.setFeatures([ deferred ]);
    map.infoWindow.show(evt.mapPoint);
  }

  function getWindowContent(graphic) {
    //make a tab container
    var tc = new dijit.layout.TabContainer({
    });

    return tc.domNode;
  }

  dojo.addOnLoad(init);


texts:
I have two grids, grid A representing observed extents of flooding and grid B representing simulated extents from a hydraulic model. The rasters have the same extent and grid size.

I want to be able to compare each pixel and create a new raster with 4 unique values (A/B/C/D equivalent in the table below) to show the difference between the two rasters:

If this was a shapefile field calculator problem I'd use a CASE logic like the one below but I'm struggling to identify a QGIS raster-based equivalent.  Can anyone suggest anything I could try?

CASE
    WHEN "Modelled" > 0 AND "Observed" > 0 THEN D
    WHEN "Modelled" = 0 AND "Observed" > 0 THEN C
    WHEN "Modelled" > 0 AND "Observed" = 0 THEN B
    WHEN "Modelled" = 0 AND "Observed" = 0 THEN A
END


texts: What Geometric model can I use for the old black and white aerial imagery?

texts:
I am working a lot of aerial photos that are very old and were scanned by a contractor. The flight were taken and snapshot over the national forest. The aerial photos are not georeference.

I am going to use them in ERDAS Imagine 2016 AutoSync to perform them with control points and run APM.

The only thing I have some concerns is do I need to have the pixel size to match from the input and the reference ?

For the aerial Photos, what geometic model do I need to use?


texts: Use Esri Portal Sign-in page for ASP.NET?

texts:
I have an ASP.NET application that is hosted on IIS that deals with map services and story maps.
I would like to federate it to an ArcGIS Portal that is federated to ADFS. I want to make it so whenever a user hits the site and is not signed into the portal, it will redirect to Portal to sign in and then back to the ASP.NET site with the user’s information. if they are signed in, then just redirect to the portal sign in and then back to ASP.NET with the security information.
I have followed the answer to the following question to achieve that, but the guide uses a custom View to sign in.
ArcGIS Portal use server-side obtained authorization in the client / JavaScript
My question is, is there a way to use the sign-in page that the Esri Portal already uses so that I don’t have to set up a new one? And can I use that to get the security information from the portal (I only need the username for records) using this method?


texts: Refresh linear referencing layer that's based on a database view, without closing ArcMap

texts:
Background:

I have a linear referencing event layer that was created using the Make Route Event Layer geo-processing tool.

The Input Route Features are in a standard feature class, but the
Input Event Table is a database view (both from the same Oracle
geodatabase).
The Event Type is set to
POINT, and the Measure Field uses the MIDPOINT field in the database view.

I've used the
   database view as the Input Event Table so that I can symbolize a point at the midpoint of my linear referencing events.

--ENG.EVENTS_MIDPOINT_VW
--Calculates the midpoint measure between the from and to measures.

SELECT
    ROAD_ID
    ,FROM_MEAS
    ,TO_MEAS
    ,(TO_MEAS-FROM_MEAS)/2 AS MIDPOINT
FROM
    ENG.EVENTS

The Problem:

When I make edits to the underlying table, they are not reflected in the linear referencing event layer, even after refreshing the map. The only way I've found to get the linear referencing event layer to refresh, is to restart ArcMap. This is not ideal.

Oddly, the attribute table of the database view (not the linear referencing event layer) does refresh just fine in ArcMap, after edits are made to the underlying table. But, for some reason, the linear referencing event layer doesn't reflect the edits.

The Question:

How can I get a linear referencing event layer that's based on a database view to reflect edits made in the underlying table?


texts: Bulk Download App 'Insuficient Previledges'

texts:
Ordered Landsat 7 data from USGS and wants to download it using the Bulk Download application. However, i get this error while running the software.

I have JRE Installed properly.
Does someone knows how to got about this Insufficient error problem.


texts: R: Reclassify a raster stack through the range of the quantiles of each layer. Generic quantiles function ; layer by layer

texts:
I have a raster stack that contains terrain bends extracted with different neighborhood windows "crosc_nxn" (e.g. 3x3, 5x5, etc ...).

crosc_nxn <- raster::stack(path.files.list)
> crosc_nxn
class       : RasterStack
dimensions  : 170, 172, 29240, 5  (nrow, ncol, ncell, nlayers)
resolution  : 30, 30  (x, y)
extent      : 794418.4, 799578.4, 7400299, 7405399  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=22 +south +ellps=WGS84 +units=m +no_defs
names       : crosc_3, crosc_5, crosc_9, crosc_11, crosc_15

I want to create a new stack containing only the extreme values of each layer using two quantile ranges ( <= 0.01 or 0.99 >=).

out_prob <- function(x, min= quantile(x, 0.01,na.rm=TRUE),  max= quantile(x, 0.99, na.rm=TRUE)){ifelse(x <= min | x>= max , x, NA)}

range_calc <- calc(crosc_nxn, fun= out_prob)
range_over <- overlay(crosc_nxn, fun= out_prob)

It appears that calc applies the function to a certain pixel and traverses all layers to return the desired value. However I want my function to sweep all the pixels of the layer and then move on to the next one.  I'm surprised that overlay has returned the same result.

On the other hand when specifying the quantile for a certain layer the result is consistent with my query, in that layer where I expressed the interval.

    # calc quantile
    > quantile(crosc_nxn, c(0.01, 0.99), na.rm=TRUE)
                       1%         99%
    crosc_3  -0.001324438 0.001358304
    crosc_5  -0.000933790 0.000937456
    crosc_9  -0.000481476 0.000486876
    crosc_11 -0.000356400 0.000362014
    crosc_15 -0.000223900 0.000228338

    # explicit range in function
    out_prob_ex <- function(x) {ifelse(x <= -0.001324438 | x >= 0.001358304, x, NA)}

    # generic
    out_prob <- function(x, min= quantile(x, 0.01,na.rm=TRUE),  max= quantile(x, 0.99, na.rm=TRUE)){ifelse(x <= min | x>= max , x, NA)}

    quantile_generic <- overlay(crosc_nxn, fun= out_prob)
    quantile_explicit <- overlay(crosc_nxn, fun= out_prob_ex)

    # comparison using only the first stack layer
    par(mfrow = c(2,2))
    plot(outlier_generic[[1]], main= "generic formula")
    plot(outlier_explicit[[1]], main= "explicit quantile")
    hist(outlier_generic[[1]])
    hist(outlier_explicit[[1]])

I need a generic function for each layer, without having to express the quantiles one by one. I thought overlay would do this and I can not figure it out. Someone has an explanation to help me see the way out.


texts: Library to produce pretty print maps including rasters

texts:
I've come across a variety of libraries for easily producing maps for printing, in a fully reproducible way using some sort of code or configuration file. For example, Kartograph looks lovely, Vincent can do some cool stuff, and basemap allows integration of geographic data with matplotlib.

However, I'm trying to produce maps which include raster data - and most of the options above only work with vector data. Are there any alternatives that will allow me to display rasters on my maps too?

Of course, I can produce the sort of maps that I want in ArcMap or QGIS fairly easily, but it is difficult to make the process automatable and reproducible. For example, if I produce one map, and then a few weeks later want to produce the same map (exactly the same view extents, exactly the same zoom, exactly the same scaling for the raster stretches) with some different data, it gets quite tricky. By doing this in some sort of automated/scripted way I'd remove a lot of these issues.

The key things I'd need are:

Produces good looking maps
Relatively concise (some of the basemap examples seem very long, and this is where something like Kartograph shines - see its documentation)
Can display rasters and vectors, with configurable symbology/display for both
Automatable

Does anyone know of a tool that will allow me to easily do this? If not, what is the best 'hacky' solution you've found? I can't be the first person to run into this problem, but I can't find a good solution online.


texts: When working with raster calculator there is an error 000539? The program failed to create a common raster map. I don't know what to do. I did the same in school without any problems. But now, on my computer the raster calculator will not make a map.


texts: Getting ALL demand points covered in location allocation analysis if facility capacities are sufficient using ArcGIS Pro

texts:
I want to match existing facilities with points of demand on the basis of individual capacity/demand and shortest distance. There is more capacity available than there is demand. The desired outcome would be to have ALL points of demand covered, I do not care how many facilities are chosen, I only need to know which ones supply which demand points.
I have tried location allocation maximise capacitated coverage, but that appears to minimise the number of facilities while never including all demand points. I've tried changing the FacilityType to 'required' or 'chosen' but the solver still rendered them 'candidates' in the process (despite close demand points not being served). For 'facilities to find' I used the total available number of facilities.
Am I doing something wrong or is this not possible, or is there a different tool to do this with?
I'm using ArcGIS Pro 2.5.2.


texts: Why does ArcGIS Online project my files right where they should be but, ArcGIS for desktop does not

texts:
When loading shapefiles into Arc Map 10.2.2 I get the message saying "unknown spatial reference", then my files are not projected where they should be and are usually too small to see regardless how much I zoom in. When I load the same files into ArcGIS Online they load fine and are projected exactly as they should.

If the spatial reference is unknown to ArcMap then how can ArcGIS Online know how to reference them?


texts: Openlayers Vector Layer Polygons not displaying properly with OSM

texts:
I have an Openlayers map displaying Acetate tiles in OSM, and I'm drawing polygon regions on a vector layer to represent the countries of the world.  This works fine for the most part, however, there seems to be a problem with where the vector layer ends and it not repeating the regions.  For example, here is the zoomed out view:

As you can see, the polygon ends/begins on the east side of Russia near Alaska.  This looks okay in this view, but if you start navigating around, it starts to get confusing.  As you can see in these two screenshots, parts of Europe or North America aren't showing even though I've drawn polygons over those regions.  As you scroll to the left and right these regions will hide and show:

This is preventing me from being able to use OSM layers in general.


texts:
I work with ArcGIS. In the picture you can see several pipelines, more over pipeID's A; 1; 2; 3
I would like to have as an output the following columns and rows:

A - 1;
A - 2;
A - 3;

As an inputfile I have a polyline file. I have taken a buffer around it of 0.5 meters and then I tried to spatially join this buffer file with the original shapeline file, using "Within". Nevertheless this does not deliver me the requested matches as I hoped.
Does anyone have an idea on how I should do this?


texts: Updating vertices Z in polyline from points

texts:
I'm trying to snap a line's vertices to several points in height. The XY have already been snaped, but the snap tool doesn't seem to get the Z value. There is one vertex for each point. Is there a way to update the vertices of the line with the Z value of the points?

I've had a similar problem before, but that time it was between points which shared the same ID. I'm not sure the vertices are able to store an ID the same way. I have acces to Arcgis Pro 2.2 and ArcMap 10.6.

Update: I tried to edit BERAs code to edit an existing fc instead of creating a new one.

with arcpy.da.UpdateCursor(existingFC, ['SHAPE@']) as cursor:
for row in cursor:
    linelist = []
    for part in row[0]:
        partlist = []
        for pnt in part:
            if (pnt.X, pnt.Y) in pointzdict:
                print("in list")
                row[0] = arcpy.Point(pnt.X,pnt.Y,pointzdict[(pnt.X,pnt.Y)])
            else:
                row[0] = pnt
    #row[0] = arcpy.Polyline(arcpy.Array(linelist), arcpy.Describe(polylinefc2).spatialReference, True)
    cursor.updateRow(row)

It gives me this error:

Traceback (most recent call last):
  File "<string>", line 13, in <module>
SystemError: <built-in method updateRow of da.UpdateCursor object at 0x00000195B48ECDF0> returned NULL without setting an error


texts: Creating a list of Rasters from another list of Rasters based on file name?

texts:
in Arcpy trying to create a subset of rasters from a list of rasters.

Basically the idea, run arcpy.ListRasters on a geodatabase, check all raster names for a user input value, then add those rasters that match to a new list. ideally, the user would be able to add the values they want as a list of values and the code iterates through each one. And actually I'm trying to get it so it will refine the new list further. but I only really need the code to do the first one.

Code so far

 listvals = arcpy.GetParameterAsText(0)
 listrasters = arcpy.ListRasters
 for raster in listrasters:
    rastname = str(raster)
    for val in listvals
       a = rastname.rfind(val,7,10)
       if a > 0:
          #add raster to newlist
       else:
          #do nothing


texts: Copying values from virtual field to non-virtual field using PyQGIS?

texts:
Using the python console, I am trying to copy values from a virtual field into another non-virtual field. Both fields are integer types but I get NULL in the non-virtual field. I followed the answer in this post: Is it possible to programmatically add calculated fields?

The following code works if both fields are non-virtual but doesn't work if copying from a virtual field. Here is the code:

from PyQt4.QtCore import QVariant
from qgis.core import QgsField, QgsExpression, QgsFeature

vl = QgsVectorLayer("C:\Users\Me\Desktop\Test\\Layer.shp", "Layer", "ogr")
vl.startEditing()

idx = vl.fieldNameIndex( 'Field_1' )             # Field_1 is the non-virtual field

e = QgsExpression( ' "Field_2" ' )               # Field_2 is the virtual field
e.prepare( vl.pendingFields() )

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

vl.commitChanges()

This works manually through the Field Calculator.


texts: How to obtain the maximum/minimum longitude and latitude for a polygon?

texts:
I have the MSA/CMSA shapefile which contains 274 metropolitan areas. How do I have to obtain the maximum/minimum longitude and latitude for an area?

So, it's like drawing a box right outside and getting the 4 corner points in the lon/lat format.

Can we do this using Python to iterate the process for all the cities and to export the coordinates to a text file as follows?:

msa     min_lon   max_lon   min_lat  max_lat
6162   -76.XXXX  -74.XXXX   38.XXXX  40.XXXX
9094   -97.XXXX  -96.XXXX   37.XXXX  38.XXXX
....
....


texts: how can i plot news stories to show how many stories include each country on a map of the world

texts:
We have a table that shows the number of news stories that include particular countries over a set period of time

there are three news channels that we have this data for and we would like to show/ plot on a map of the world so that it will easily show the concentration of the stories per country whereby the marker/ number etc for each country shows clearly the concentration of stories.

Such that a country with 100 stories per period is shown as a higher concentration than one with just 10 stories per period.

We have our own data as a spread sheet

Is this possible using QGIS?


texts:
How can I open a KML Network Link in QGIS?
I have googled this question and searched for similar questions on Stack Exchange. I am not asking how to open a KMZ in QGIS.
I don't see any references to trying to open a KML Network Link in QGIS.


texts:
Do you know how we can show Google Maps as a base map in ArcObject stand alone map?

For example is there any WMS address from Google Maps that we can use it to connect to it or something like that?

I have used open layers before and it has ability to set Google Maps as base map... But I want to use Google Maps in my stand alone application without using GeoServer and ...

Do you have any suggestions?


texts: How do I calculate the distance between multiple associated points in QGIS?

texts:
I have two large lists of points. List A contains the location of a farmers farmstead and list B contains locations of pastures. Both lists have a column that will identify the farmer that the pasture or farmstead belongs to. There is only one farmer per farmstead but there may be several pastures linked to each farmer. Each farmer does not have a pasture.

So to my question, how do I calculate the distance between each pasture and its farmers farmstead? I know how to calculate the distance between two given points but since these are rather large sets it is not possible to go through the lists point by point. Is there a way to set this up without making a matrix of all possible combinations of points? The lists contain ~65000 and ~2000 points which means that such a matrix would become quite cumbersome to handle.

I'm using QGIS 2.18 and my projection (EPSG:3006) uses meters as the map unit.

This question is similar but using a different software.


texts: Does QGIS DB Manager plugin lock PostGIS transactions?

texts:
I do most of my geospatial analyses in PostGIS 2.0 using PgAdmin III. I use QGIS 1.8 for visualizing my data but also for querying my database using the DB Manager plugin.

Sometimes (quite often actually...) when I import tables from PostGIS into QGIS using DB Manager and then I run queries on these layers in pgAdmin, they keep on running forever. The only way to bring the queries to stop is to shut QGIS down. Not sure if this is a bug in the plugin or if my PostGIS settings are wrong but I must say I'm well annoyed by this behaviour.

Many thanks in advance.


texts: I am working on writing basic expressions to control the Marker Properties. My previous question on the syntax for color was graciously answered , but I soon realized that I need to know the syntax for all the properties like , size, outline, angle. Is there a place I can find them all


texts: Can feature layer in MXD with joined data from non-geodatabase be published as feature service?

texts:
I have a table stored in an SDE geodatabase for all point features with an additional GEOMETRY_ID. A few non-spatial data tables are stored in an ordinary MS SQL Server database. They can be joined by the common GEOMETRY_ID and the layer is stored in a MXD.

The MXD file couldn't be published as feature service, so Add, Update and Delete operations cannot be applied to this layer. Therefore, I am planning to move the non-spatial database to SDE geodatabase and published as feature service.

I have questions as below:

Is that true that a feature layer in MXD that has joined data from non-geo database cannot be published as feature service in ArcGIS Server?
In order to allow operations to those non-spatial tables, say Add, Update and Delete records, is it a good practise to store them into SDE geodatabase and published as feature service?


texts: Digitizing and converting PDF to shapefile using QGIS

texts:
Is there a way to convert the pdf to shapefile in QGIS?
I have plenty of PDF files (I think they are exported from AutoCAD) and I need to georeference them and convert them to shapefile. It would be easier if I do it automatically (Not digitize everything by hand at least), since we have plenty of files.


texts: How to apply custom OSM map style in OpenLayers?

texts:
I have a small web mapping application based on PostGIS/GeoServer/OpenLayers.

I'm overlaying few WMS layers served by GeoServer and using OpenStreetMap as a background.

In OpenLayers OSM is simply called as:

var osm = new OpenLayers.Layer.OSM("OSM");
map.addLayer(osm);

How could I display OSM basemap in shades of grey to accentuate my WMS layers?

Is there a way to use existing styles from CloudMade's Style Editor for that purpose?


texts: Problem editing, clipping, and saving large shapefile QGIS 3.4

texts:
I'm working with a very large shapefile in QGIS 3.4. I only need one section of the layer (near the main city, not the areas outside it) and I'm trying to delete the large sections which I don't need. Unfortunately, QGIS keeps crashing (it's fine with all other operations). I got the following error message on the last crash:

Stack Trace

RaiseException :
_TypeMatch :
RtlCaptureContext :
QgsEllipseSymbolLayer::properties :
QgsLineString::points :
QgsCurvePolygon::coordinateSequence :
QgsGeometry::convertPolygon :
QgsGeometry::asMultiPolygon :
QgsGeometryValidator::run :
QThread::start :
BaseThreadInitThunk :
RtlUserThreadStart :

QGIS Info
QGIS Version: 3.2.3-Bonn
QGIS code revision: 9b176802e5
Compiled against Qt: 5.9.2
Running against Qt: 5.9.2
Compiled against GDAL: 2.2.4
Running against GDAL: 2.2.4

System Info
CPU Type: x86_64
Kernel Type: winnt
Kernel Version: 10.0.17134

So far I've tried the following:

1) Deleting vertices in small batches
2) Deleting vertices in large batches
3) Exporting to GeoJSON and reloading the GeoJSON and repeating the above
4) Creating a new shapefile of a bounding polygon for the area I need and using the "clip" geoprocessing tool to isolate the section of the original shapefile I need in a new layer
5) Setting the canvas to the area I need and saving as a new layer using the "map view extent" option under "Layer" > "Save As"

None of the above have worked. The first three all result in QGIS crashing, the fourth does not create a correctly clipped new layer (to get this to run I have had to change "Invalid Features Filtering" to "Ignore features with invalid geometries"), and the fifth doesn't save the map view extent but instead saves the whole of the original shapefile.

I'm thinking option five is the best, but I cannot get it to save just the "map view extent". I have searched the forums for "map view" "map view extent" "canvas extent" but couldn't find anything about the issue.  I've also tried this on two different computers but the results are the same.

Does anyone have any ideas?


texts: QGIS - Pop-up menu not displayed when right click on a layer in layer tree view

texts:
I'm developing a custom QGIS application using the QGIS API. I've been using the API provided with QGIS Brighton (Version 2.6.1) and have been developing the application in Python in a Windows 7 environment. I have provided a basic map canvas for my application with an accompanying layer tree view. While most of the functionality appears to be working fine, I've been unable to get a pop-up menu to display when right-clicking on a layer in the layer tree view. I have followed the scheme suggested by Martin Dobias in his "QGIS Layer Tree API (Part 3)" post on the Jan 30 2015. I've also looked through the qgisapp.cpp source code, and although written in C++, my Python application seems to be following the same approach as the qgis code. I have included contextMenuEvent event handler with is called when I right click on the layer tree. However the createContextMenu() function in the Menu Provider is not being called. Do I need to explicitly call it from contextMenuEvent event handler. This doesn't appear to be what is required from Martin's post.

I've included a couple of code snippets so that you can see what I am doing.

class FeweMenuProvider(QgsLayerTreeViewMenuProvider):

    def __init__(self, view):
        QgsLayerTreeViewMenuProvider.__init__(self)
        print "FeweMenuProvider.init()"
        self.menuProviderView = view

    def createContextMenu(self):
        print "FeweMenuProvider.createContextMenu()"
        if not self.menuProviderView.currentLayer():

            print "No Layer - Got in here at least"
            return None

        contextMenu = QMenu()
        contextMenu.addAction("Show Extent", self.showExtent)
        print "Layer clicked."
        return contextMenu

    def showExtent(self):
        print "FeweMenuProvider.showExtent()"
        r = self.menuProviderView.currentLayer().extent()
        QMessageBox.information(None, "Extent", r.toString())

class FeweMainWindow(QMainWindow):

    def __init__(self, parent=None):
        super(FeweMainWindow, self).__init__(parent)
.
.
.

            self.rootNode = QgsProject.instance().layerTreeRoot()
            self.model = QgsLayerTreeModel(self.rootNode, self)
            self.model.setFlag(QgsLayerTreeModel.AllowNodeReorder, True)
            self.model.setFlag(QgsLayerTreeModel.AllowNodeRename, True)
            self.model.setFlag(QgsLayerTreeModel.AllowNodeChangeVisibility, True)
            self.model.setAutoCollapseLegendNodes(10)

            self.feweLayerTreeView = QgsLayerTreeView(self.centralWidget)
            self.feweLayerTreeView.setGeometry(10, 40, 171, 521)
            self.feweLayerTreeView.setModel(self.model)

            # Create popup menu for layer tree
            #
            contextMenuProvider = FeweMenuProvider(self.feweLayerTreeView)
            self.feweLayerTreeView.setMenuProvider(contextMenuProvider)

.
.
.
    def contextMenuEvent(self, event):

        print "Got into contextMenuEvent()"
        QMainWindow.contextMenuEvent(self, event)
        #menu = self.feweLayerTreeView.menuProvider().createContextMenu()
        #print menu
        #if menu is not None:
        #    print "About to show menu."
        #    menu.show()
        #event.accept()


texts: Filtering Small Areas of Shapefiles Using R

texts:
I'm trying to use the method proposed by the @ialm in the Question "How to speed up the plotting of polygons in R?"

There, he loads the maps using getData from the raster package, I tested and it works fine.

I downloaded files from the Brazilian Institute of Geography and Statistics - IBGE and one of the ShapeFiles could be this one, then I use readOGR from rgdal.

I'm facing an error on method 2

  # Get the main polygons, will determine by area.
getSmallPolys <- function(poly, minarea=0.01) {
  # Get the areas
  areas <- lapply(poly@polygons,
                  function(x) sapply(x@Polygons, function(y) y@area))

  # Quick summary of the areas
  print(quantile(unlist(areas)))

  # Which are the big polygons?
  bigpolys <- lapply(areas, function(x) which(x > minarea))
  length(unlist(bigpolys))

  # Get only the big polygons
  for(i in 1:length(bigpolys)){
    if(length(bigpolys[[i]]) >= 1 && bigpolys[[i]][4] >= 1){
      poly@polygons[[i]]@Polygons <- poly@polygons[[i]]@Polygons[bigpolys[[i]]]
      poly@polygons[[i]]@plotOrder <- 1:length(poly@polygons[[i]]@Polygons)
    }
  }
  return(poly)
}

Particularly here if(length(bigpolys[[i]]) >= 1 && bigpolys[[i]][4] >= 1)

Error in if (length(bigpolys[[i]]) >= 1 && bigpolys[[i]][4] >= 1) { :
  missing value where TRUE/FALSE needed

As for me bigpolys[[i]][4] is always NA, so I started to suspect the structure of files are not the same.

So far, as starter in this area, I have tried a combination of query methods ( slots, head, str, attributes, class, unclass ) trying to uncover the differences that would justify the error in my case.

> class(mex)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
> class(mapa)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

> slotNames(mapa)
[1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"
> slotNames(mex)
[1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"

I need help in 2 directions, what are the best methods to compare 2 shapefiles like in the example above and detect their structural differences, and if possible an appropriate solution for clipping the small polygons.


texts: ArcGIS - where is the zip & unzip toolbox available for download?

texts: ESRI mentions a zip & unzip toolbox available for download (http://webhelp.esri.com/arcgisserver/9.3/java/index.htm#geoprocessing/python_2126079299.htm) Unfortunately the direct link provided is out of date & I can't seem to find it searching through ESRI's sites or anywhere else.  Does anyone know where ESRI hid it?


texts: Reprojecting a shapefile with ogr2ogr when a line is crossing the international date meridian

texts:
I have a line shapefile with a local transverse mercator projection
("+proj=tmerc +lat_0=50 +lon_0=172.5 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs") and I want to reproject it into EPSG:4326 with ogr2ogr.
The line feature is crossing the international date meridian, as you can see in the image.

This particular situation seems to be a problem for ogr2ogr. In fact, the reprojected output line is crossing the entire world, and this is clearly wrong.

This is the command I used for the reprojection:

ogr2ogr -f "ESRI Shapefile" track_linedate_pr.shp track_linedate.shp -t_srs EPSG:4326


texts: ArcGIS Pro Project Template using Absolute Paths

texts:
i'm trying to create a project template with ArcGIS Pro that uses a shapefile stored on my companies shared network drives. This file contains basic parcel data for our county and is around 500MB in size.

The problem is that when I package this into a project template, instead of referencing the share drive location of the shapefile, it packages the shapefile with the project template file. Not only does this make the project file balloon to about 150MB in size, but once used it extracts out a copy of the full shapefile into the new project folder.

Is there any way to force Pro to save out the absolute path for the shapefile and not simply copy it to the new project template?


texts:
How do you get the text information to pop out with the information extracted from the Attribute Table when you hover over it or a right-click on it or something akin to that?  Example is in the picture supplied.

Are you able to extract certain columns from the AT or all the columns are automaticaly included?

When doing this, are you limited with options on how the pop-up is done?

Thank You, Zoran


texts: Showing overland flow route on DEM in QGIS?

texts:
I've been asked to produce a risk map at work which will show the possible flood route from sewers that surcharge (flood) during heavy rainfall. This is basically where I know the x-y coordinate (OS National Grid Ref. as I'm in the UK) of a number of manholes we know sometimes flood, and I have a 2 metre x 2 metre Digital Elevation Model (Raster) from the UK Environment Agency that covers the area in question.

Basically I need a method that looks at the elevation at the point each manhole occurs at, and then looks at the elevation of each raster square around that point. If a square is at the same or lower elevation then assume that square is at risk of flooding, and for each of THOSE squares, look around again and identify further squares at the same or lower elevation. Keep going until there are no more directly connected squares at the same or lower elevation, and this is the total area that could be at risk.

It doesn't need to account for volumes or the speed at which this might happen, just literally all connected raster squares at the same or lower elevation as the origin point.

I know how to use QGIS fairly well and it's all I have access to at work.


texts: Algorithm or command line tool to decimate point cloud of terrain points?

texts:
I need to take a larger (more dense) than needed list of lidar survey points (longitude, latitude, and elevation) for terrain definition and decimate it based on a 2 dimensional grid. The idea would be to end up with points based on a NxN (i.e. 1 meter x 1 meter) dimension grid using the longitude, latitude (x,y) values, therefore eliminating the points that are more than are needed. The goal is to determine what the elevation is at each point in the grid after the decimation, not use elevation as part of the decimation rule itself.

I would like to either code/script this myself in a scripting or programming language if I can start with a decimation algorithm or use a command line tool from a project that may already exist that can do this that can run on Ubuntu and called from our application as system call. The approach should not require using a GUI based type of software or tool to solve this. It needs to be part of an automated set of steps.

The data currently exists in a tab separated values file but I could load the data into a sqlite database file if using an database/sql query driven algorithm would be better/faster. The ideal scripting language would be ruby or python but can be any really and if there exists C/C++/C# libraries for this already then we could wrap those for our needs.

Ideas?


texts: In CartoDB I want to import geo-referenced points based on where I took photos and then, I'd like to see the photos / photo URL when I hover over the points.


texts: RemoveMarkers() is not working for PruneCluster markers in Leaflet

texts:
I have a set of markers clustered with PruneCluster plugin. The data of these markers is rendered from GeoJSON layers, the way to do it is like this:
    var propsAmp = amphibiens.features
    var propsMam = mammiferes.features
    var propsOis = oiseaux.features
    var propsPois = poissons.features
    var propsRep = reptiles.features
    var propsSb = sibes.features
    var propsPrc = parcs.features
    var propsZh = zonesHumides.features
    var propsArb = arbres.features

        var propsEspeces = []
        propsEspeces.push(propsAmp)
        propsEspeces.push(propsMam)
        propsEspeces.push(propsOis)
        propsEspeces.push(propsPois)
        propsEspeces.push(propsRep)
        propsEspeces.push(propsSb)
        propsEspeces.push(propsPrc)
        propsEspeces.push(propsZh)
        propsEspeces.push(propsArb)

        //console.log(propsEspeces)

        for (i in propsEspeces) {
            for (j in propsEspeces[i]) {
                var marker = new PruneCluster.Marker(propsEspeces[i][j].geometry.coordinates[1], propsEspeces[i][j].geometry.coordinates[0], {
                    data: propsEspeces[i][j].properties,
                    icon: L.icon({
                        iconUrl: "markers/" + propsEspeces[i][j].properties.category_image,
                        iconSize: [60, 60], // width and height of the image in pixels
                        shadowSize: [35, 20], // width, height of optional shadow image
                        iconAnchor: [12, 12], // point of the icon which will correspond to marker's location
                        shadowAnchor: [12, 6], // anchor point of the shadow. should be offset
                        popupAnchor: [0, 0]
                    })
                })
                switch (propsEspeces[i][j].properties.GrandeCat) {
                    case "AP":
                        marker.category = AIRES_PROTEGEES;
                        break;
                    case "Veget":
                        marker.category = VEGET;
                        break;
                    case "Faune":
                        marker.category = FAUNE;
                        break;

                }
                pruneCluster.RegisterMarker(marker);
            }
        }

        pruneCluster.PrepareLeafletMarker = function(leafletMarker, data) {
            leafletMarker.on('click', function(evt) {

                console.log('clicked')
            });
            leafletMarker.data = data
            leafletMarker.setIcon(data.icon);
        };

        pruneCluster.ProcessView();
        map.addLayer(pruneCluster);
        console.log(pruneCluster)

        var markers = []
        var remove_marker = document.getElementById('remove_marker');
        remove_marker.addEventListener('click', function() {
            for (var i in pruneCluster.Cluster._markers) {
                if (pruneCluster.Cluster._markers[i].data.data.category == 'Arbres') {
                    markers.push(pruneCluster.Cluster._markers[i])
                }

            }
            console.log(markers)
            pruneCluster.RemoveMarkers(markers)
        })

I'm getting the properties of each GeoJSON layer then pushing it in a global array, then assign to each marker it's proper data and icon.
I added a button to remove markers from map fulfilling a certain condition. So I'm storing the markers result of the condition in an array then passing it as parameter to the RemoveMarkers() function yet I'm getting nothing not even an error. I would like to know what am I missing and how can I handle this?


texts: Esri JS API 4.7 - Changing the view spatial reference in place

texts:
Asked this already in GeonNet, but asking here as well as I can get valuable help from both sites.

I'd like to change the view.spatialReference property in place to be able to switch between basemaps with different spatial reference.
All the basemaps I use should be of type TileLayer.

I can actually do what I want if initialize a new view with the spatial reference I want and then add the new basemap.

However, I'd like to do it witout recreating the view from scratch if possible.

Something like doing view.spatialReference: 32632 or similar with the extent property.

The problem is, when I do it, the console raises an viewpoint error and I don't know how to handle it.

The error is:

[esri.views.MapView] #viewpoint incompatible spatialReference {"latestWkid":3857,"wkid":102100} with view's spatialReference {"wkid":32632}

Is there a way to reinizialize the viewpoint when the view changes? I don't see in the docs any reference to setting the spatialreference for the viewpoint, I guess this is set internally when the view is inizialized and this gets the app confused when changing the spatial reference of the view, while its viewpoint "remembers" it as the old one...

I am aware that the viewpoint overrides the extent property in the view class.

Here is the codepen example to hopefully make you understand better my problematic.

Just, be aware that the above example may work or not depending on your browser (I am using FireFox and it works for instance).


texts: Creating several new fields at once with QGIS field calculator in batch mode in the same layer/file

texts:
How is it possible to create several new fields at once in the same layer, using QGIS field calculator in batch mode? There is a similar question, but using Python - I'm interested in a solution without Python.
In QGIS, using the field calculator in batch mode (from processing), it is possible to calculate several fields at once. Let's suppose a very simple case just for demonstration purpose here: an attribute table with a field named value. I now want to create 5 new fields (red box in the screenshot):

output_1: "value"+1
output_2: "value"+2
etc.

This works easy in batch mode of field calculator, using the increasing row-number: each row added will create a new field. And all input parameters can be set with Autofill... / Calculate by Expression… accordingly, writing expressions that generate the necessary values needed for each column (field calculator input parameters) in the batch dialog.
The problem now is: used in this way, QGIS will create a separate file for each row / each new field calculated: output_1, output_2 etc. If I define the same output path for each row (last column in the batch dialog window, see screenshot), than QGIS will iterate over all rows and each time overwrite the output of the previous row. Thus the result will contain only the last output, in this example with 5 rows/5 fields to create: output_5: "value"+5. When definining separate output filenames (like using increasing numbers, generated by the variable @row_number: output_1.gpkg, output_2.gpkg etc.), I get the desired result. But I get aseparate file for each field that is created.
The question: How is it possible to batch create new fields in the same layer or the same file - ideally adding the newly created fields to the attribute table of the input layer?
Screenshot: Field calculator in batch mode: field name, expression to calculate and output path are generated by an expression with Autofill... / Calculate by Expression…:


texts: I am using ArcGIS 10.1.  I inserted a legend however the items in the legend are too large because of the size i set them in the data view.  How do i scale the items in the legend so they will be the same size in the layout tab as they are in the data view tab.


texts: QWC2 - How to calculate WMTS Resolutions

texts:
I am trying to add this WMTS as a background layer in QWC2.

https://geoegl.msp.gouv.qc.ca/carto/wmts?SERVICE=WMTS&REQUEST=GetCapabilities&VERSION=1.3.0

In the GetCapabilities results I see 22 scale denominators that I am trying to convert into resolution values for my themesConfig.json file in QWC2.

I have been searching on internet for days now but I can't find any info or formula on how to do this.

Does such a formula actually exist ?  Or is there a way I haven't seen to make QWC2 figure out the resolutions itself ?


texts: Closing out a 'add polygon' operation in the interactive MapWinGIS ShapeEditor

texts:
I am creating a small app that uses a MapWinGIS control and allows a polygon layer to be edited. Users can add/delete/edit polygons. All is working fine, but the operation to finish a new polygon (Ctrl + Left Mouse Button) is not from intuitive. We would like to implement this functionality by using the more intuitive "left click on the first vertex" which most drawing & geospatial packages use.

How can I do this? I should be able to trap the left mouse event, but need to compare the coord with the start (which is stored where?) and then initiate the "finish this polygon" operation. I see the StartEdit() operation in the ShapeEditor documentation, but no equivalent FinishEdit()!?


texts: Custom cuts using R and `choropleth` in package `GISTools`

texts:
I am plotting some data on a shapefile using R and choropleth in packageGISTools
I would like to modify my code below in order to assign custom breaks.

Using the World Borders Dataset, download from: http://thematicmapping.org/downloads/TM_WORLD_BORDERS_SIMPL-0.3.zip

library(GISTools)
library(rgdal)

shape_file <- readOGR(dsn ="/home/hamid/TM_WORLD_BORDERS_SIMPL-0.3",
                      layer = "TM_WORLD_BORDERS_SIMPL-0.3")

shades <- auto.shading(shape_file@data$POP2005, cutter = quantileCuts,
                       n=5, cols = brewer.pal(5, "Blues"))
choropleth(shape_file, shape_file@data$POP2005,
           shading = shades)
choro.legend(70.00000, 32.50000, shades)

Can anyone help please?

UPDATE

No errors in the code I have posted. What I would like to do is to use custom cuts instead of quantile cuts.
I have trie the following:

cuts= c(0.1, 0.2, 0.3, 0.4, 0.6)
shades <- shading(shape_file@data$POP2005, breaks = cuts)
choropleth(shape_file, shape_file@data$POP2005, shading = shades)

which give the following error:

Error in findInterval(v, shading$breaks) :
vec' must be sorted non-decreasingly`

Note that the values for the breaks are specific to my case which is not using the same data on the example.


texts: Splitting OSM files into multiple smaller pieces using QGIS?

texts:
I have a large osm file of a city. I want to split this file into multiple smaller pieces (e.g 50) and retain all smaller pieces for further use rather than just a portion of it. Is anyone aware of anything out there on the Windows platform that can do this?

I'm expecting I'll have to just get the max lat/lon and  min lat long and just divide the numbers to get the coordinates for manual piece downloads from OSM, but this seems a shame if I have the full city osm already.

I am using QGIS.


texts: Resources to integrate ESRI gis to SAP ?

texts: I heard to this term to integrate GIS with SAP.What does this actually mean? Where can i find resources to integrate GIS to SAP. I have no idea on where to start and to proceed further?