texts: How to open image file saved with JAI in Geotools?

texts:
I need to find NDVI and other vegetation indices and then find changes in time with tif images.

I've been using Geotools for this. I took 2 tiff images, then using QGIS I changed them to PNG because original ones are huge. Then taking those 2 PNGs I created NDVI image using JAI.

The problem is that I couldn't show this file with Geotools. You can look at the code here:

private void constructNDVI(File band3File, File band4File)
        throws IOException {
    PlanarImage iNIR = JAI.create("fileload", band3File.getAbsolutePath());
    PlanarImage iVIS = JAI.create("fileload", band4File.getAbsolutePath());
    // The pixels on those images must be processed as floating-point
    // values!
    ParameterBlock pbConvert = new ParameterBlock();
    pbConvert.addSource(iNIR);
    pbConvert.add(DataBuffer.TYPE_DOUBLE);
    PlanarImage NIR = JAI.create("format", pbConvert);
    pbConvert = new ParameterBlock();
    pbConvert.addSource(iVIS);
    pbConvert.add(DataBuffer.TYPE_DOUBLE);
    PlanarImage VIS = JAI.create("format", pbConvert);
    // Calculate the denominator (NIR-VIS).
    ParameterBlock pbd = new ParameterBlock();
    pbd.addSource(NIR);
    pbd.addSource(VIS);
    PlanarImage denominator = JAI.create("subtract", pbd);
    // Calculate the numerator (NIR+VIS).
    ParameterBlock pbn = new ParameterBlock();
    pbn.addSource(NIR);
    pbn.addSource(VIS);
    PlanarImage numerator = JAI.create("add", pbd);
    // Calculate the NDVI.
    ParameterBlock pbNDVI = new ParameterBlock();
    pbNDVI.addSource(denominator);
    pbNDVI.addSource(numerator);
    PlanarImage ndvi = JAI.create("divide", pbNDVI);
    DisplayNBImage disp = new DisplayNBImage(ndvi);
    String filename = new Date().getTime() + ".png";
    JAI.create("filestore", disp.getSurrogateImage(), filename, "PNG");
    File file = new File(filename);
    showNDVI(file);
}

You can find DisplayNBImage here.

private void showNDVI(File file) throws IOException {
    AbstractGridFormat format = GridFormatFinder.findFormat(file);
    AbstractGridCoverage2DReader reader = format.getReader(file);
    ...
    final MapContent map = new MapContent();
    map.setTitle("ImageLab");
    map.addLayer(rasterLayer);
    ...
    JMapFrame frame = new JMapFrame(map);
    ...

I can open the file I've created in constructNDVI method in Windows Photo Viewer or QGis. But 2nd line in showNDVI gives

Exception in thread "main" java.lang.UnsupportedOperationException: Trying to get a reader from an unknown format.

I've tried bmp, tiff but no luck. How can I overcome this?

In fact, proposing non-Jai way of the same thing would be awesome.


texts: Where does the output of a grid volume run in the Saga tool/ Qgis end up?

texts:
I have tried to do a grid volume run with the Saga grid calculus - grid volume function. After the run the window closes and nothing seems to happen, no message, no new layer or information appear.
I use a raster .tif as input.

After a run, where does the output normally appear?


texts: How to remove legend items with PyQGIS without removing the layer?

texts:
I am trying to generate maps in pdf files by using a pre-defined .qgz layout. I am already able to load layers and adjust color scale, etc. Since the legend is synchronized with the map layers, all the layers I have loaded are listed in the legend.

Using PrintComposer in QGIS I can turn off synchronization and remove the unwanted items in legend. Is there a way to do it automatically in Python? I tried to use QgsLegendModel, but it also removes the layers from map canvas.

.

projectInstance= QgsProject.instance()
projectLayoutManager = projectInstance.layoutManager()

for comp in projectLayoutManager.printLayouts():
    for item in comp.items():
        if(isinstance(item, QgsLayoutItemLegend)):
            model = item.model()
            while model.rowCount()>1:
                 model.removeRow(0)

            break


texts: Accumulate values from raster using ArcPy

texts:
I wrote the script:

import arcpy
if arcpy.CheckExtension("Spatial")=="Available":
    arcpy.CheckOutExtension("Spatial")
    print("Spatial")

from arcpy import env
from arcpy.sa import *
env.workspace = "D:\Semestr 2 MGR\GeoinformatykaB\Cwiczenia\C07\DaneSkryptGISv02\Dane\Orto"
arcpy.gp.ZonalStatisticsAsTable_sa("punkt","OBJECTID","CentrumGdanska.tif","C:/Users/ilona/Documents/ArcGIS/Default.gdb/ZonalSt_punkt2","DATA","MAXIMUM")

arcpy.CheckInExtension("Spatial")
print "end"

I have a problem. This script is failed.
This is error:

ERROR 000824: The tool is not licensed.
Failed to execute (ZonalStatisticsAsTable).

Could anyone check if it works?

The content of the task:

Create a script to accumulate values from the raster. The tool comes with:

raster
the point from which the zone will rise
barrier value of the maximum sum of values obtained from raster cells.


texts: Why has my ArcGIS Enterprise username changed case?

texts:
Something rather peculiar has occurred on our federated, on premesis ArcGIS Enterprise instance.

Last week I had logged in with my username which is case sensitive containing caps and lower case letters as a named user. This week that username does not work unless I only use lower case letters.

From what I've read so far, you cannot change your username- but mine has changed and I'm curious as to why/how it happened. Also it seems to have only affected me and the change has been carried through to groups and content so that the lower case version of me still owns everything Mr uppercase did and is still owner and member of the appropriate groups.

I'm concerned that this is a symptom of some deeper problem with the server. On the face of it the server logs don't show anything useful (no warnings) but I will need to do a deep dive to see if anything crops up. Has anybody had a similar experience?


texts: Combine raster layers with same spatial extent but different data to multivariable value attribute table

texts: Does QGIS allow to combine several raster layers having the same spatial extent, but containing data of different dates, to one single raster with an attribute table of several columns (as with the COMBINE tool of the Spatial Analyst in ArcMap)?


texts: Need to set centerAt to true before calling locate() method

texts:
I need to allow the user to change the extent while the Locator Widget is tracking.  I defined the locator widget with the following code

var locateButtonLayer = new GraphicsLayer();
app.map.addLayer(locateButtonLayer);
var geoLocate = new LocateButton({
    map: app.map,
    graphicsLayer: locateButtonLayer,
    useTracking: true,
    setScale: false,
    centerAt: false
}, "LocateButton");
geoLocate.startup();
app.geoLocate = geoLocate;

If I don't set setScale and centerAt to false some browsers will reset the extent every few seconds.  If you user wants to re-center on their location the following code is called.

function customLocate() {
   require(["dojo/Deferred"], function(Deferred) {
       app.map.setScale(9000).then(function() {
           app.geoLocate.centerAt = true;
           app.geoLocate.locate().then(function(result) {
               app.geoLocate.centerAt = false;
           });
       });
   });

}

My issue is the locate() method will only center on the user's extent if I set a breakpoint on the app.geoLocate.centerAt = false; line.  Setting a breakpoint on any other line does not help.  I assume this has to do with the locate() returning a promise but I can't figure out how to fix the issue.  The code appeared to  run fine when I put it in the Locator Widget example sandbox at Locator Widget Sandbox

EDIT:  The code works fine if I add an alert as shown below:

function customLocate() {
   require(["dojo/Deferred"], function(Deferred) {
       app.map.setScale(9000).then(function() {
           app.geoLocate.centerAt = true;
           app.geoLocate.locate().then(function(result) {
               alert("here");
               app.geoLocate.centerAt = false;
           });
       });
   });

}

EDIT: I tried the following:

  app.map.setScale(9000).then(function() {
   app.geoLocate.centerAt = true;
    on(app.map, "update-end", function(){
        app.geoLocate.locate().then(function() {
            on(map, "update-end", function(){
                app.geoLocate.centerAt = false;
            });
        });
    });

});

This solved the issue with the locate() method executing before the map was done updated but the app.geoLocate.centerAt = false line never gets executed so the map keeps re-centering every few seconds.  Evidently the locate() method prevents the update-end event from occurring when it is in tracking mode.


texts: Append dynamic stats on attribute table in QGIS Composer

texts:
I have an attribute table containing records like Town, inhabitants, N° male, N° Female like this:

When I compose a map in QGIS I use "Show only visible feature" command resulting in a view with only few records.
I need to add a last dynamic row or separate custom table with statistics (based only on visible features) to summarize my table with:

Inhabitans sum,
n°male sum
n°female sum
male/female ratio

How can I do in QGIS Map Composer?


texts: Creating Groups of points based on proximity in QGIS?

texts:
I have a dataset with around 900 points and I want to cluster them into groups of 44 based on proximity.

Is this possible in QGIS?

I have tried to use the convex hull tool with no success.


texts: How to get distance between a point and a list of highways?

texts:
I'm new to GIS. I'm trying to determine whether a list of GPS points are within a highway/road (in terms of meters). I have downloaded a OSM file and populated my database using osm2pgsql.

I have the following tables:

                     List of relations
 Schema |            Name            |   Type   |  Owner
--------+----------------------------+----------+----------
 public | geography_columns          | view     | map_user
 public | geometry_columns           | view     | map_user
 public | planet_osm_line            | table    | map_user
 public | planet_osm_point           | table    | map_user
 public | planet_osm_polygon         | table    | map_user
 public | planet_osm_roads           | table    | map_user
 public | raster_columns             | view     | map_user
 public | raster_overviews           | view     | map_user
 public | spatial_ref_sys            | table    | map_user
 public | gps_coordinates            | table    | map_user
 public | gps_coordinates_id_seq     | sequence | map_user

My coordinates table looks like so:

id | longitude  | latitude  |                        geom

I populated my table with a few longitude/latitude points, then ran the following function to generate a geometry value:

UPDATE vehicle_coordinates SET geom = ST_SetSRID(ST_MakePoint(longitude, latitude), 4326);

I initially tried using the PostGIS function ST_DWITHIN like so:

SELECT gps_coordinates.id, highway.name
FROM gps_coordinates, planet_osm_line highway
WHERE highway.highway IS NOT NULL
AND ST_DWITHIN(gps_coordinates.geom, highway.way, 5);

But no results were returned. I don't think I can just use meters in this function as I'm using latitude/longitude as my base unit. I'm stuck and not sure where to go from here. Can I use St_Distance in my WHERE clause? It states that it returns the distance in meters. But I'm not getting any rows returned:

SELECT gps_coordinates.id, highway.name
FROM gps_coordinates, planet_osm_line highway
WHERE highway.highway IS NOT NULL
ST_DistanceSphere(vehicle_coordinates.geom, highway.way) <= 5;

It could be that my data is bad (I got it from GEOFABRIK, so should be good) or that I'm confusing my units here.


texts:
I would like to embed the live versions of BA infographics into story maps created on the newer StoryMap interface.  I can export infographics from BA directly into the old story maps, and I can export images or pdfs of the infographics into the new version, but I do not seem to be able to get the live version of the BA infographics into the new story maps.
The only solution I have thought of is to export the graphics as HTML, host them, and then embed that in the storymap.
Is there a simpler way that I'm not seeing?


texts: Recursive filter and clipping of ee.ImageCollection

texts:
I'm trying to source average NDVI raster images from July of every year from 2008 through 2018 for the USA. I need to export this data to link with PHI offline. Exporting the images for the entire country isn't allowed so I need to break it up into 4 regions. This will result in 44 raster images exported from GEE. I'm having a really hard time automating this process since for loops are not allowed.
Below is an example I'm trying with just two regions and one date. My approach was to try and create a ee.List composed of two ee.Dictionaries then map the drawDateRegion function over that ee.List. It's not working. I think that some of the issue might be in the fact that I can't seem to reference elements of the dictionaries contained in the flatRegionsAndDates2 list.
I'm willing to try another approach if somebody has a suggestion. Essentially I need to clip and filter an ee.Image collection over all combinations of regions 1-4 and Julys 2008-2018.
//Getting MODIS Image Collection
var modis = ee.ImageCollection('MODIS/006/MOD13Q1');

//New program to get all regions and years
//Import all 4 US regions as FeatureCollections
//Importing region
var region1 = ee.FeatureCollection("users/fossaal/us_region1_10km_buffer");
var region2 = ee.FeatureCollection("users/fossaal/us_region2_10km_buffer");
//var region3 = ee.FeatureCollection("users/fossaal/us_region3_10km_buffer");
//var region4 = ee.FeatureCollection("users/fossaal/us_region4_10km_buffer");

//Creating date ranges
var startDate = ee.Date('2016-07-01');
var dateRange = ee.DateRange(startDate,startDate.advance(1,'month'));

//Create arrays
var allRegions = [region1,region2]
var allDates = [dateRange]
//Functions for mapping
var arrayOfArrays = function(region){
  var regionDatePair = function(date){
    return {
      region:region,
      date:date
    }
  }
  return allDates.map(regionDatePair)
}
var regionsAndDates = allRegions.map(arrayOfArrays)

print(regionsAndDates)

var flatRegionsAndDates = ee.List(regionsAndDates).flatten()

print(flatRegionsAndDates)

//Make elements of list into Dictionaries
var toDict = function(a){
  return ee.Dictionary(a)
}

var flatRegionsAndDates2 = flatRegionsAndDates.map(toDict)
print(flatRegionsAndDates2)

//Creating function to clip and filter
var drawDateRegion = function(dateRegionObject){
  print(dateRegionObject)
  var dataset = modis.filter(ee.Filter.date(dateRegionObject.get('date')));
  //Reducing MODIS images to averages for all bands
  var redDataset = dataset.mean();
  // Clip to to US borders.
  var clippedRedDataset = redDataset.clipToCollection(dateRegionObject.get('region'));
  // Select only NDVI band
  var ndvi = clippedRedDataset.select('NDVI');
  return ndvi
}

var output = flatRegionsAndDates2.map(drawDateRegion)


texts:
I am using a feature service to show data on the map like so:

var featureTableLayer = new FeatureLayer("http://12.122.105.741/arcgis/rest/services/test/FeatureServer/1", {
    mode: FeatureLayer.MODE_ONDEMAND,
    showLabels: true,
    visible: true,
    syncSelection: true,
    zoomToSelection: true,
    outFields: [ "objectid", "mesto_sifra", "mesto_ime", "ulica_ime", "ime_na_ulica_opstina", "novo_ime_na_ulica", "zabeleska_opstina"],
                    id: "streetsLayer"  });
    map.addLayers([featureTableLayer]);

I also have a map service where I have labels

var mapLayer = new FeatureLayer("http://12.122.105.741/arcgis/rest/services/test/MapServer/1", {
    mode: FeatureLayer.MODE_ONDEMAND,
    showLabels: true,
    visible: true,
    syncSelection: true,
    zoomToSelection: true,
    outFields: [ "objectid", "mesto_sifra", "mesto_ime", "ulica_ime", "ime_na_ulica_opstina", "novo_ime_na_ulica", "zabeleska_opstina"],
                    id: "streetsLayer"  });
    map.addLayers([mapLayer]);

I am adding both layers to the map, both the feature and map services but I am not seeing the labels on the map but when I open the map service via the ArcGIS Online Map Viewer it loads the map with the labels. What could be the issue?


texts:
I am a novice to QGIS tasked with updating other developers' plugins from QGIS 2 to QGIS 3.

Here's the code I need to replace:

legend = self.iface.legendInterface()
index = legend.addGroup('Quadrant Shapes', False)
self.quad_group = index
self.iface.mapCanvas().setRenderFlag(False)
mlr = QgsMapLayerRegistry.instance()

According to GitHub, iface.legendInterface() has been replaced with iface.layerTreeView(), so I can replace the first line with legend = self.iface.layerTreeView() but the next line fails because the QgsLayerTreeView object does not have an addGroup() attribute.

The general idea is to add a group to the interface and then append some layers to it.

UPDATE:

With the input from Fran Raga and the discovery that moveLayer() is no longer available, the following solution worked for me:

root = QgsProject.instance().layerTreeRoot()
self.quad_group = root.addGroup('Quadrant Shapes')

self.north_layer = QgsVectorLayer([shp], "North_Quad", "ogr")
self.quad_group.insertChildNode(0, QgsLayerTreeLayer(self.north_layer))

self.south_layer = QgsVectorLayer(os.path.join([shp], "South_Quad", "ogr")
self.quad_group.insertChildNode(0, QgsLayerTreeLayer(self.south_layer))
...etc


texts: Splitting polygons at narrowest part using R?

texts:
I have a dataset of polygons like this one:

I would like to split these polygons in separate parts at their most narrow location (if they have one). For example the two small polygons should not be split. So I would need to identify narrow locations and then split the polygon there.

How could this be done using R?


texts: Control flow plug in of GeoServer 2.18.2 doesn't show on GeoServer Tools

texts: I already downloaded and added to Lib folder as all other extensions of GeoServer that I needed but I am having issue loading the control flow extension which I confirmed that it is same version as my GeoServer 2.18.2 (war file version installed on Java 64 bit) .


texts: I've loaded several shapefiles into QGIS and I'd like to increase the transparency across all of those layers while keeping them their original color.  Can this be done without changing each layer individually?


texts: Simulate water runoff to get water height over each pixel of DEM

texts:
Here is my problem that I am trying to deal with.

I have given a height map of around 512X512 resolution(It is Geo-referenced and I have the meta-data for this). And I have given that this area receives some x mm/hr rainfall for an hour(so I know the total volume of water that fell). Now I am assuming that the surface is concrete, so there is water loss to the ground, and all the water will reside inside this area(assume that the area is enclosed by infinity walls). Further there is no loss of water from evaporation(atmosphere has 100% humidity).

Now I want to find out at 5 or 10 minutes time interval, how much of water will each pixel on the height map will have. It's like snapshots of water levels(height of water), at various time intervals.

I have looked into GRASS r.topmodel[1], but could really understand if what I want is actually achievable from this function. Further there is a topmodel implementation in R, and I am not sure about this as well, if this can help me, in achieving what I want(all I see is a lot of parameters of flowing water).

Can someone help me figuring out how can I use these above functions to help my cause. If not, help by pointing out some (easy to implement) algorithm or some other implementation or library that I can use.

[1] http://grass.osgeo.org/grass70/manuals/r.topmodel.html


texts: How can I get the first selection of a layer in a QgsMapLayerComboBox?

texts:
I'm developing a very simple QGIS plugin to calculate raster average by using GDAL. It uses a QgsMapLayerComboBox and a filter to select only raster layers. However, it needs to select twice a layer at the Combo Box to produce the expected results. The snipped code is:

.
.
.
def initGui(self):
        """Create the menu entries and toolbar icons inside the QGIS GUI."""

        icon_path = ':/plugins/AverageRaster/AverageRaster.png'
        self.add_action(
            icon_path,
            text=self.tr(u'To calculate raster average'),
            callback=self.run,
            parent=self.iface.mainWindow())

        self.wcb = QgsMapLayerComboBox(self.dlg)
        self.wcb.setFixedWidth(220)
        self.wcb.move(140,40)
        self.wcb.setFilters( QgsMapLayerProxyModel.RasterLayer )
.
.
.

def run(self):
    """Run method that performs all the real work"""
    #My code starts here
    layer = self.wcb.currentLayer()
.
.
.
    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result:
        # Do something useful here - delete the line containing pass and
        # substitute with your code.
        #My code starts here
        message = "the raster average of " + filename + " is = " + str(average)
        QMessageBox.information(None, "Raster Average", message)
        #My code ends here

At the next image, my Plugin is determining the average of ndvi.tif raster. If I want the average of utah_demUTM2.tif, the expected value is obtained only at the second chance. How can I get the average at the first selection of a layer at the QgsMapLayerComboBox.

Afterward that I tried out a lot of possibilities, I found a possible solution but, I don't like it very much because I have to click on icon plugin again and not on a OK button.

At the 'Property Editor' of Qt Designer, I unmarked the 'OK' button and I modified the label "Click on OK to get raster average" for this one: "Click on plugin icon again to get raster average".

After reloaded the plugin, I observed that the order in the Combo Box was alphabetic.

When I selected anyone of these layers (see above image) and I clicked on icon plugin, its behavior was always as expected and the principal Dialog never closed.

I changed the QMessageBox for a QTextBrowser and put again the OK button at the widget Dialog. In this case, I could observe that the QgsMapLayerComboBox has the expected behavior when I select the raster layer. The problem is when I clicked on the OK buttom. The Dialog was closed and I need to launch the plugin again to see the results. It is possible to keep the window Dialog open when the OK button is clicked?


texts:
I have a shape file with nature preservation areas (Biosphärenreservate) in Germany. I think the projection  must be something like gauss-krüger zone 4 Pulkovo or something like that. When I reproject it to UTM-33N WGS84 it reprojects but is not exactly accurate. The borders don't quite match my other shape (the other file)

http://www.projektidee.org/images/extern/biosph-res-2010.prj.zip
http://www.projektidee.org/images/extern/gebiet_UTM33.dbf.zip

Could anyone here help me with that problem. That would be nice!

best regards
desputin


texts: QGIS hide the field in the feature properties form if the field's value is empty

texts:
When clicking on a feature with the  tool and obtaining its properties, in the form that is shown, I would like to view only the fields that are filled (not empty). In some layers I have 50 fields but only 4 are filled.

Is it possible to do it using Edit Widgets?


texts: Shapefiles from Nationalgrid without a .prj file - what projection are they in?

texts:
Im working with these shapefiles:

http://www2.nationalgrid.com/uk/services/land-and-development/planning-authority/shape-files/

They show pipelines & electricity grid of the UK. They do not have a .prj file included in the shapefile.

They display fine in QGIS when layered on top of each other, but I need to know which projection they are in so I can display them on a world map.

The x and y values are very large so I tried the Universal transverse Mercator projection "WGS84/UTM zone 30N ESPG:32630" because the UK is in zone 30.
This placed the shapefiles at the right scale, but just above the equator by Ghana (48 degrees too far south and 3 degrees too far west).

This seems to imply that the shapefiles' coordinates are zeroed half way between UTM zone 30U and 31U.
http://en.wikipedia.org/wiki/File:Utm-zones.jpg

Does anyone know which projection I should be using?

Thanks for your help.


texts: Which Projection system to use for area calculation in Bangladesh (in South Asia)?

texts:
I am vectorizing a raster layer using ArcScan. I need to calculate the area of the vectorized feature after the vectorization process.

What projection should I use for the vectorized polygon layer to calculate the area of the polygon?

Geographic location for the layer is: Bangladesh (in South Asia)


texts: Are Landsat 8 data from ESPA corrected and calibrated?

texts:
I am trying to figure out Landsat 8 data corrections and calibrations. Actually, I don't really know how this can be done and which data needs what kind of correction and calibration for further processing. Basically new for these processing. So, I prefer to download already corrected and calibrated Landsat scenes.

I have downloaded Landsat 8 OLI/TIRS C1 Level-1 data from https://glovis.usgs.gov/. It seems that there is no correction and calibration. I got following data in the package.

LC08_L1TP_129026_20140712_20170421 = *

        *_01_T1_b1 to b11
        *_01_T1_bqa
        *_01_T1_MTL
        *_01_T1_ANG

Also I downloaded same scene from https://espa.cr.usgs.gov/. It supposed to be corrected and calibrated. I guess it is. Because the value range of band data are different. Following are the additional data plus the above package, cause I ticked TOA calibrated, and SR, ST and angle corrected.

*_01_T1_sr_band1 to band7
*_01_T1_toa_band1 to band9
*_01_T1_pixel_qa
*_01_T1_radsat_qa
*_01_T1_sr_aerosol
*_01_T1_sensor_azimuth_band4
*_01_T1_sensor_zenith_band4
*_01_T1_solar_azimuth_band4
*_01_T1_solar_zenith_band4

Are the data downloaded from https://espa.cr.usgs.gov/ corrected and calibrated to be used for further analysis directly?

I am wondering, if the band data already corrected and calibrated, what are the additional data for (i.e. sr_band1-band7, _01_T1_toa_band1-band9, _01_T1_sensor_azimuth_band4 and _01_T1_sensor_zenith_band4)?


texts: Getting minimum distance of points with scalable algorithm

texts:
I have a number of geographical coordinates with latitudes/longitudes:

id  | name          | lat       | lon
1   | New York      | 40.730610 |  74.0060
2   | Philadelphia  | 75.1652   | 39.9526
3   | Washington DC | 77.0369   | 38.9072
... | ...           | ...       | ...

I could calculate all distances between my points using standard tools, but I'm only interested in the minimum distance of one point to all others. But that would require N*N (N: number points) calculations and is too slow.

A heuristic/approximate approach would be better.

Any ideas?


texts: Control Points Disappear from ArcMap Link Table

texts:
I georeferenced a scan map and updated my georeferencing after adding the control points. But after I started digitising, and saving edits, I soon realised that when I open the link table on the Georeference toolbar, the table is empty and my control points gone, however, my work is still spatially correct.

What could be the reason and how do I solve this?


texts:
I have a set of points with latitude/longitude coordinates and some value that I want to map using QGIS software.

When I import my points and apply a graduated style I obtain the following map. On the latter all the points are visible and the overall effet is not nice from my point of vue.

Is there a way to improve that ? I'd like to have a smooth and continuous field.


texts: In cartodb, how to automatically update the_geom with new longitude and latitude?

texts:
I am currently posting new table row using the cartodb sql api with the following commands:
update {table_name} set geom = ST_Transform(ST_SetSRID(ST_MakePoint(latitude, longitude),4326)
latitude and longitude are taken from numeric columns in my table. However the_geom is not updated, and instead an error is received.

I am aware I can use 'Georeference your dataset' within cartodb to update the the_geom...but I don't want to do this manually whenever I add new coordinates.

Is there a way to automatically update the_geom when I add new new longitude and latitude coordinates?


texts: Why can't I add most plugins under Ubuntu 12.04?

texts:
I loaded QGIS into 12.04 with much trouble but it works. Although I do not seem to be able to load most plugins.  I go to fetch plugins and most are not loaded so I try to load some and they will not load including many that I used with 11.10.  I click on fetch and then plugins and then highlight a plugin and like Table Manager and click Install Plugin and the box appears:

Failed to unzip the plugin package. Probably it's broken or missing from the repository. You may also want to make sure that you have write permission to the plugin directory:
/home/hebert/.qgis//python/plugins

This is true for 90% of the plugins.  Help!


texts: How to extract road centreline from delaunay triangulation using open source QGIS?

texts: From some months I have been solving a problem but I cannot finish. My question is, "How can I extract a road centreline from delaunay triangulation using open-source QGIS?" Actually, I am using Chugiak.


texts: Making lookup table for performing Lookup(inRaster, lookupField)?

texts:
I want to replace the Digital number (DN) values in a raster file with some other values which are present in a text file;in the text file, I have a table which it's first column are all possible DN and on the next column I have corresponding values to be replaced.

First of all, I want to know which ways make it possible to use this table to replace DN values with Python specially arcpy?

I want to use Lookup(inRaster, lookupField) in ArcMap but I do not know how can I introduce my separate text file to this statement as a lookup table?


texts:
I want to compute the shortest path between two points
origin=(34.052299,-118.1812657)
destination=(34.141132,-117.853179).

Below is my workflow.

Step 1: Find the nearest nodes.

origin: the following returns 590279764

SELECT id
FROM nodes
ORDER BY st_distance(the_geom, st_setsrid(st_makepoint(-118.1812657, 34.052299), 4326))
LIMIT 1

destination: the following returns 123007528

SELECT id
FROM nodes
ORDER BY st_distance(the_geom, st_setsrid(st_makepoint(-117.853179, 34.141132), 4326))
LIMIT 1

Step 2: Compute shortest path between origin and destination.

SELECT *
FROM pgr_dijkstra(
    'SELECT gid as id, source, target, length AS cost FROM ways',
    590279764, 123007528, false, false
);

However when I execute this I get the following error.

ERROR:  Start vertex was not found.

********** Error **********

ERROR: Start vertex was not found.
SQL state: XX000

Why is this happening?


texts: Problem creating feature layer with synchronized edits on both ArcGIS Online and ArcMap

texts:
I am working on a project where I have created a shapefile of building footprints on ArcMap 10.4, and I need the features and their attribute tables to be editable inside of ArcGIS Online for those using it in the organization. I also need the edits to transfer over to ArcMap. I have the data stored in the organization's Server, but I have not been able to export the data into ArcGIS Online with editing and have the changes made in AGOL transfer over to the shapefile in ArcMap successfully.

If anyone knows a good way to make this possible, that would be fantastic. I apologize if this is very simple as I am new to the industry.


texts: Definitive list of Nominatim types and categories?

texts:
Is there a definitive list somewhere on the interwebs of the category and type values in Nominatim?

I've been doing some Googling but found nothing as yet.


texts: Adding WMTS from Geoserver getcapabilities via OpenLayers?

texts:
I am new to web mapping and web in general.

I'm trying to add the default WMTS Italy mosaic that can be found in Geoserver (nurc:mosaic) to my map.

I think there is a way to do it by providing get capabilities wmts link, not xml, not Rest.. (This one is with xml and works for me : http://openlayers.org/en/latest/examples/wmts-layer-from-capabilities.html)

(This one is with Rest which I don't want to implement : http://openlayers.org/en/latest/examples/wmts.html)

I get the wmts link simply right clicking on Geoserver's main page WMTS 1.0.0 and copy the link location. Here it is: http://localhost:8081/geoserver/gwc/service/wmts?REQUEST=GetCapabilities

And here is my naive code:

<!DOCTYPE html>

    <title>WMTS</title>
    <link rel="stylesheet" href="https://openlayers.org/en/v4.6.4/css/ol.css" type="text/css">

I get this message from Fiddler:

<ExceptionReport version="1.1.0"
xsi:schemaLocation="http://www.opengis.net/ows/1.1
http://geowebcache.org/schema/ows/1.1.0/owsExceptionReport.xsd"><Exception
exceptionCode="InvalidParameterValue"
locator="TILEMATRIX"><ExceptionText>Unknown TILEMATRIX
10</ExceptionText></Exception></ExceptionReport>

I tried to add 'EPSG: ' to my matrixIds [z] value, (Saw it on : wmts layer in Openlayers from Geoserver) it didn't work.

I tried to change matrix Id's new Array values up to 21, didn't work.

In my humble understanding GeoServer goes into the columns rows that don't exist for some reason. For example say I have columns ranging from 70 to 80 and rows ranging from 100 to 110, when I provide wmts get capabilities link and run the code when I check from the Fiddler, link goes to tilecolumn: 250 and tilerow 300something etc. These files don't exist so there is a mismatch.

What should I learn to get over these kind of issues ?


texts: QGIS "transform error" after downloading from OpenStreetMap

texts:
Beginner's question. I am working on this tutorial:
http://www.qgis.org/en/docs/training_manual/vector_analysis/basic_analysis.html

After selecting WGS84 UTM zone 34S in Project Properties, I receive several of these transform error messages:

Simplify transform error caught: forward transform of
(0.366535, 1.570812)
PROJ.4: +proj=longlat +datum=WGS84 +no_defs +to +proj=utm +zone=34 +south +datum=WGS84 +units=m +no_defs
Error: latitude or longitude exceeded limits

Can anyone guide me as to next steps? From other questions I've read on this site I assume it's a long/lat issue. Many thanks.


texts:
I am trying to upload a photo into a folder in ArcCatalog, yet I haven't been able to find where and how to do that.

I am using ArcGIS 10.1 (ArcGIS for Desktop) through my school's Citrix receiver.

I have a MacBook which apparently does not work with ArcGIS very well?


texts: Authenticate to geoserver using openlayers 2

texts:
I want load layers from geoserver using openlayers in php page.
How I should do it? Maybe can give ane code examples?
Layers not public, need authentification.


texts: Showing only 1 layer at a time in CartoDB

texts: Is there a way in the free CartoDB to set the layer selector so that you can only select 1 layer at a time?  When one layer is selected, deselect all other layers.


texts: Categorized symbols are not working in QGIS 3.4

texts:
I'm trying to generate categorized symbols for a polygon layer in QGIS 3.4. The code below should generate random colour polygons for each record of 'run' layer and I'm not getting any error message. The layers toolbar shows the coloured categories created, but the 'run' layer is not visible at all.

    ls = QgsProject.instance().layerStore()
    r_layer = ls.mapLayersByName('run')[0]

    # get unique values
    fni = r_layer.fields().indexFromName('id')
    unique_ids = r_layer.dataProvider().uniqueValues(fni)
    QgsMessageLog.logMessage("sstyle for run layer..." + str(unique_ids))
    # define categories
    categories = []
    for unique_id in unique_ids:
        # initialize the default symbol for this geometry type
        symbol = QgsSymbol.defaultSymbol(r_layer.geometryType())
        symbol.setOpacity(0.5)

        layer_style = {}
        layer_style['color'] = '%d, %d, %d' % (randrange(0, 256), randrange(0, 256), randrange(0, 256))
        layer_style['outline'] = '#000000'
        symbolLayer = QgsSimpleFillSymbolLayer.create(layer_style)

        if symbolLayer is not None:
            symbol.changeSymbolLayer(0, symbolLayer)
        category = QgsRendererCategory(unique_ids, symbol, str(unique_id))
        categories.append(category)

    renderer = QgsCategorizedSymbolRenderer('id', categories)
    # assign the created renderer to the layer
    if renderer is not None:
        r_layer.setRenderer(renderer)
    r_layer.triggerRepaint()


texts: Specify a fields datatype while creating a new field in geodataframe with expression?

texts:
I would like to know how to set the field type of ['denominator'] in a geodataframe manually in this expression to "String":

geo_df['denominator'] = df[["basalareap","basalareas","basalaread"]].sum(axis=1)

Where the fields "basalareap","basalareas","basalaread" contain floating point numbers

Automatic it assign the data-type "Double" (when check in QGIS), I tried to add the datatype in the first bracket, geo_df['denominator', type= "String"] but this throws an error.

I'm using Windows 10, Python 3.7


texts: Link QGIS attribute table with Web Image?

texts:
I have some layers, and I would like to link the field photograph from the attributes table to an image (Web URL).

How can I do that?  [Qgis 2.18 - Windows 10]

Something like...

I finally found the answer


texts: SAGA Raster Calculator, ifelse doesn't return specified value exactly

texts:
I am using the SAGA Raster Calculator in QGIS 2.6 to change NaN values to zero and that part is working great, however the specified 'if' value is being altered slightly. I have developed a work around but it seems a bit clumsy.

For example, I am using the following formula:

ifelse(eq(a,10), 10, 0)

yet when the the temporary file is created in QGIS the result ('Result' in the image below) appears as a raster with values from 0 to 9.99.

Through trial and error I have found that by using the following formula I can achieve the result I am looking for but it seems a bit clumsy.

ifelse(eq(a,10), 10.01001, 0)


texts: Visually clipping a search area (MCE) (QGIS)

texts:
I have some search areas to refine using a variety of vector inputs. I can clip these to get my end product but it is a messy process. Is it possible to use blending modes to visually "knock out" unsuitable areas without actually doing the clipping?
Inputs are vector polygons, system is QGIS 2.8 Mac.


texts: Nearest Neighbor - return object ID in R

texts:
I have a csv file that contains the GPS coordinates of about 100 villages and their unique ID. I have an other csv file with the GPS coordinates of some potential destinations and their unique ID.
I found a way to measure the shortest distance between each village and the nearest destination using get.knnx. However the get.knnx function only provides me with the distance, it does not provide me with the ID of the village nor the ID of the destination. For now I tryed this:

g1 = get.knnx(coordinates(points_period1$X), coordinates(villages$hid), k=1)

Here is my code :

library(sp)
villages <- read.csv("/path/villages.csv", header=TRUE)

head(villages,15)
     hid exact_lat exact_lon
1      1  16.97587  7.992230
2      2  16.97371  7.998500
3      3  16.96942  7.984560
4      4  16.95704  7.965890
5      5  18.74667  7.400229
6      6  18.73254  7.401260
7      7  18.71110  7.336340
8   7020  18.71110  7.336340
9   7060  18.71110  7.336340
10     8  16.31698  8.013590
11     9  15.62224  7.892040
12 10040  15.52151  7.034985
13 10050  15.52151  7.034985
14 10060  15.52151  7.034985
15 10070  15.52151  7.034985

coordinates(villages) <- c('exact_lon', 'exact_lat')
proj4string(villages) <- CRS("+init=epsg:32633")
villages = spTransform(villages, CRS("+init=epsg:32632"))

head(villages,15)
     hid
1      1
2      2
3      3
4      4
5      5
6      6
7      7
8   7020
9   7060
10     8
11     9
12 10040
13 10050
14 10060
15 10070

names(villages)
[1] "hid"

points_period1 <- read.csv("/path/points_p1.csv", header=TRUE)

head(points_period1,15)
     X layer        x        y optional
1    1     1 6.625000 17.56250     TRUE
2    2     1 6.633929 17.56250     TRUE
3   79     1 6.616071 15.59821     TRUE
4   80     1 6.625000 15.59821     TRUE
5   81     1 6.633929 15.59821     TRUE
6   84     1 6.633929 15.58929     TRUE
7  117     1 5.294643 15.50893     TRUE
8  128     1 5.116071 15.50000     TRUE
9  129     1 5.125000 15.50000     TRUE
10 130     1 5.133929 15.50000     TRUE
11 133     1 5.285714 15.50000     TRUE
12 136     1 6.160714 15.50000     TRUE
13 144     1 5.107143 15.49107     TRUE
14 145     1 5.116071 15.49107     TRUE
15 146     1 5.125000 15.49107     TRUE

coordinates(points_period1) <- c('x', 'y')
proj4string(points_period1) <- CRS("+init=epsg:32633")
points_period1 = spTransform(points_period1, CRS("+init=epsg:32632"))

head(points_period1,15)
     X layer optional
1    1     1     TRUE
2    2     1     TRUE
3   79     1     TRUE
4   80     1     TRUE
5   81     1     TRUE
6   84     1     TRUE
7  117     1     TRUE
8  128     1     TRUE
9  129     1     TRUE
10 130     1     TRUE
11 133     1     TRUE
12 136     1     TRUE
13 144     1     TRUE
14 145     1     TRUE
15 146     1     TRUE

names(points_period1)
[1] "X" "layer" "optional"

require(FNN)
g1 = get.knnx(coordinates(points_period1), coordinates(villages), k=1)
[1] "nn.index" "nn.dist"

write.csv(g1, file = "/path/dist_2011.csv")


texts: Creating space at intersection of two links

texts:
I have a shapefile consisting of two links which are crossing over each other.

I want to create an empty space at the point of intersection of the two links.

As of now, I have been manually doing this like this:

Split the line at two points near the intersection. And then remove that feature.
Similarly carry out this for the other line.

So by doing this, I create a empty space at the point of intersection.

I want to know whether there is any better way to achieve this using a plugin or any other way?

Can a plugin can be created specifically for this task ? (I'm ready to write one, if this is achievable.)

Below is an image where one link touches other link. Now, I want to create an empty space there.

This is the output I'm expecting:


texts: TopoJSON: how to calculate scale and translate?

texts:
Context:

I'm using the topojson file format along with d3. Since I could not find a library in Python that handles this format, I'm writing one. I've got the topojson to other formats working.

Question:

When converting other formats to topojson, how do I calculate the ideal value for the topojson transforms scale and translate?

In order to save space (among other tricks), the format uses relative integer offsets from the previous point. For example, this is the topojson file for Aruba:

{
  "type": "Topology",
  "transform": {
    "scale": [0.036003600360036005, 0.017361589674592462],
    "translate": [-180, -89.99892578124998]
  },
  "objects": {
    "aruba": {
      "type": "Polygon",
      "arcs": [[0]],
      "id": 533
    }
  },
  "arcs": [
    [[3058, 5901], [0, -2], [-2, 1], [-1, 3], [-2, 3], [0, 3], [1, 1],
     [1, -3], [2, -5], [1, -1]]
  ]
}

The same information as a GeoJSON feature looks like this:

{
    "type" : "Feature",
    "id" : 533,
    "properties" : {},
    "geometry" : {
        "type" : "Polygon",
        "coordinates" : [[
                 [-69.9009900990099, 12.451814888520133],
                 [-69.9009900990099, 12.417091709170947],
                 [-69.97299729972997, 12.43445329884554],
                 [-70.00900090009, 12.486538067869319],
                 [-70.08100810081008, 12.538622836893097],
                 [-70.08100810081008, 12.590707605916876],
                 [-70.04500450045005, 12.608069195591469],
                 [-70.00900090009, 12.55598442656769],
                 [-69.93699369936994, 12.469176478194726],
                 [-69.9009900990099, 12.451814888520133]
            ]]
    }
}

It is easy to see how the former is more compact than the later. In order to obtain the absolute coordinates, the function is:

def arc_to_coordinates(topology, arc):
    scale = topology['transform']['scale']
    translate = topology['transform']['translate']

    x = 0
    y = 0
    coordinates = []

    for point in arc:
        x += point[0]
        y += point[1]
        coordinates.append([
            x * scale[0] + translate[0],
            y * scale[1] + translate[1]
        ])

    return coordinates

Reading topojson is easy enough:

>>> arc_to_coordinates(topology, topology['arcs'][0])
[[-69.9009900990099, 12.451814888520133],
 [-69.9009900990099, 12.417091709170947],
 [-69.97299729972997, 12.43445329884554],
 [-70.00900090009, 12.486538067869319],
 [-70.08100810081008, 12.538622836893097],
 [-70.08100810081008, 12.590707605916876],
 [-70.04500450045005, 12.608069195591469],
 [-70.00900090009, 12.55598442656769],
 [-69.93699369936994, 12.469176478194726],
 [-69.9009900990099, 12.451814888520133]]

Update:

Reading Mike Bostock's code, I saw this:

kx = ((Q - 1) / (x1 - x0)) if x1 - x0 else 1
ky = ((Q - 1) / (y1 - y0)) if y1 - y0 else 1
...
scale = [1/kx, 1/ky]

In the example, given the bounding box in the form (x0, y0, x1, y1), Q seems to be around 6 and 12 for the values [0.036003600360036005, 0.017361589674592462].

>>> x0, y0, x1, y1
(-70.08100810081008, 12.417091709170947, -69.9009900990099, 12.608069195591469)
>>> Q = 6
>>> 1 / ((Q - 1) / (x1 - x0))
0.036003600360035644
>>> Q = 12
>>> 1/ ((Q - 1) / (y1 - y0))
0.017361589674592892

I know the Q has to do with quantization, any idea where the Q values near 6 and 12 came from?


texts: ST_AsPNG query has wrong position on Mapbox

texts:
I am trying to import GeoTIFF files into a PostgreSQL database using raster2pgsql.
I use the following query to import the data into my database:
raster2pgsql -s 4326 -c -P -I -C -t 256x256 *.tif public.my_raster_table > out.sql

The import seems ok. Then, to work with Mapbox with a z/x/y format, I build a request.
For example, for 18/133274/88180 (Lille, France) I create the following request
SELECT ST_AsPNG(ST_Union(rast)) As new_rast FROM my_raster_table WHERE ST_Intersects(rast, ST_Transform(ST_Segmentize(ST_MakeEnvelope( 336628.6725679152, 6556921.160358975, 336781.5466244854, 6557074.034415545, 3857), 38.21851414255798), 4326))

The data is retreived, but my raster file has a wrong position, and it is repeated multiple times (I try to replace a road).
I know my ST_Intersects is right as it works well for vector tiles.
Where did I go wrong ?


texts: ArcGIS Collector with ArcGIS Portal webmap failing on attachments

texts:
I have a web map that has been created on ArcGIS Portal. Our Portal is federated. The feature layer (service) and web map have all the required settings based on the documentation of ArcGIS Collector app. When logged in to the ArcGIS Collector app on iPhone (iOS), we can see the webmap and collect points, the problem is we can't attach pictures. When we take a picture for the collected point and click "Submit", the following error occurs: "Update Failed: 1 attachment could not be sent".

We are using the Classic Collector app, however, I also tried downloading the new version, the same error occurs.

I've tried to simulate the same workflow over Blue Stacks (android) Collector on PC, and everything seems to be working and attachments uploading to the web map.

I also have tried my ArcGIS Online account with another web map where I tried to attach a photo using my iPhone and Collector app, it also worked perfectly without giving me any errors.

However, we must get it to work on ArcGIS Portal for my team, as we don't have access to additional ArcGIS Online user accounts right now, and probably, there should be no issue like this working wherever...

Anyone had this problem or has any idea what could have gone wrong?

Attached you will find an error screenshot from iPhone, and also what the ArcGIS Server Manager puts as a Warning when we try to attach a photo.


texts: Which coordinate system does MapTiler use for its "Photoshop Overlay for Google Maps" page?

texts:
I am attempting to use this MapTiler webpage to create a background image that I can position my image over and then use GDAL to generate tiles. My process is:

Position the map on the MapTiler page over the area I'm adding an overlay to.
Save the static map image and the four coordinate values representing west, south, east, and, north.
Bring the saved image into Photoshop and use it as a background to create my overlay image. The new image always has the same dimensions as the background.
I then use gdal_translate to create a GeoTIFF from the overlay image and the saved coordinates. Here is an example of the parameters used. This command runs without error.

gdal_translate -of GTiff -a_srs EPSG:900913 -gcp 0 0 west north -gcp width_of_your_input_map_file_in_pixels 0 east north -gcp width_of_your_input_map_file_in_pixels height_of_your_input_map_file_in_pixels east south your_input_file_name.png your_output_file_name.png
I then use gdal2tiles.py to generate the tile set. Here is an example of that command:

gdal2tiles.py -p mercator -z zoom_range name_of_output_file.png

Which returns the error:

ERROR 1: Too many points (2704 out of 2704) failed to transform,
unable to compute output bounds.

I believe the problem is that I'm substituting in the coordinates from the MapTiler webpage which are in this form I'm not familiar with.

-4891969.810251281 -26369557.159749157
95295571.90369494 73817984.55419706

So, what coordinate system am I looking at here? How can it be converted into something gdal_translate understands?


texts: I installed both postgresql and postgis using brew in mac. How can I now install pgrouting. I couldn't install it using brew. Any suggestions?


texts:
I've been looking all over the place and I can't find the download!

Where do I find the ArcGIS 10.1 SP1 for Server Web Adaptor?


texts:
I can't seem to find a tool to merge files from the active TOC into one file...

The datamanagement --> Merge Shape files only allows you to choose a directory of shapefiles to merge which is really inconvenient as you then have to place only the required files in one directory etc or use the "select by layers" and then choose files which is a pain if the files are in multiple directories.

The Geoprocessing Tools --> Difference on the other hand just lets you choose files from a drop down of files in the TOC. This is so much easier.

Is there a version of merge in QGIS that works like "Difference"?


texts: Changing the scale of the screen display in QGIS using PyQGIS

texts:
I'm trying to answer my own question asked here, but which has not found favour: In QGIS save a view like a bookmark based on center and scale

I thought I might be able to work backwards, by defining screen center (centre) and screen scale, then working backwards to get the geographic extents. Note I'm a Python novice and have never tackled PyQGIS before.

Screen centre turned out to be trivial (so long as points were in correct units for the CRS):

canvas = qgis.utils.iface.mapCanvas()
center = QgsPointXY(44, 42)
canvas.setCenter(center)

I naively thought I could do a similar trick with setScale, but this gives the error 'QgsMapCanvas' object has no attribute 'setScale' no matter what I try.

I know I can change the display scale using lines like this:

canvas.resize(QSize(3920,4080))

However, this changes the scale by changing the window extent, which is not what I want. How can I do this?


texts: SLD - Is it possible to create a Rule that stylise only the first and last point of a Line?

texts:
I've read the SLD doc, but I can't find any clue on how to create an SLD file to add a graphic symboliser on the first and last point of a line.

Does someone know how to do that ?

Thanks.


texts: How to plot millions of points in web maps with python, postgis, leaflet, and etc.?

texts:
I am building a web map for plotting millions and billions of point data. It seems impossible because millions of point data might blow the memory and browser. I am wondering whether it is possible to plot that much data on a web map using for example leaflet or google earth engine.

My idea right now is convert all of these data into a raster image, so the plotting might become possible. But I have two questions here, how to convert the point data into raster, and whether it is still possible to do gis with the point data?

If there are other ways to plot the points (probably point cloud?), please provide me some clues. I am totally stuck right now.


texts: Editing PostgreSQL geodatabase table using ArcPy?

texts:
I'm trying to update a table using a python script and it's not working.  The table resides in a postgres geodatabase and is versioned.  I have also verified my user permissions (I'm the admin).  The version is public (i.e. not protected).  The table properties show that the current user has select, insert, update, and delete privileges.

Here's my command

with arcpy.da.Editor(arcpy.env.workspace) as edit:
    with arcpy.da.InsertCursor("service", ["service_id"]) as cur:
        cur.insertRow([0])

The error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 3, in <module>
RuntimeError: The requested operation is invalid on a closed state [ugdb.sde.service]


texts:
I wonder if there is a way to combine fill and line style for a polygon seperately.

I tried with simple fill + fillstyle or line pattern fill but both resulted in lines without no background coulor.


texts: Retain band names in mosaic dataset Extract Bands

texts:
I have a collection of multi-band Sentinel-2 rasters. From these, I apply the Extract Band Raster Function. On the resulting mosaic dataset, I need to retain the initial band names. For example, if I extract B4 (red) and B8 (ir) bands, in the resulting mosaic dataset, I will have B1 and B2 as names, but they can be easily misinterpreted as B1 (coastal aerosol) and B2 (blue) bands.

Is there a way in ArcMap or Python to retain the original band names?


texts: How to get back closed panes in QGIS 1.8?

texts:
I was using the map composer in QGIS 1.8 in Mac OS X 10.6.8 today and closed a couple of the object attribute panes on the right side of the map composer. However, I can't find a way to get those panes back? Can someone tell me how to do it? It is probably pretty obvious but I also can't do this for panes I close in the main window, too....

Thanks!


texts:
maybe you can help

I have a string:

"width"=>"6","name:be"=>"Калінавы завулак","name:ru"=>"Калиновый переулок","surface"=>"asphalt","int_name"=>"Kalinavy zavulak","maxspeed"=>"40"

necessary as a distinct attribute : Kalinavy zavulak

my test:

Split( [Your_Field_Name] , "int_name"=>",)(1)

But something goes wrong


texts: polygon Self intersection find in arcobjects

texts:
how to get self intersection polygons in arcobjects could you please tell me how to get already i got self intersection for polyline features as per below code but it is not working on polygon features

if (intF.Shape.GeometryType == esriGeometryType.esriGeometryPolyline)
{
    IPolyline pline = intF.Shape as IPolyline;
    ITopologicalOperator3 ptopo = pline as ITopologicalOperator3;
    esriNonSimpleReasonEnum reasen = esriNonSimpleReasonEnum.esriNonSimpleSelfIntersections;
    ptopo.IsKnownSimple_2 = false;
    if (!ptopo.get_IsSimpleEx(out reasen))
    {
    if (reasen == esriNonSimpleReasonEnum.esriNonSimpleSelfIntersections)
    {
        selfinterfind = true;
    }
    }
    if (selfinterfind == true)
    {
    dt.Rows.Add(intF.OID, intF.Class.AliasName, "", "Self Intersection Found");
    }
}


texts: Count the number of raster tiles in each polygon in ArcGIS

texts: I have a raster from the TRMM 3B43 product that contains a grid of 0.25 x 0.25 tiles (in degrees). I also have a shapefile with a number of districts. I want to know how many raster tiles are included in each district (polygon). Ideally, I'd like to know it as a continuous measure (i.e. polygon A contains 2.43 tiles) but if it only counts complete tiles its also OK.


texts: Pyqgis: Style Layer After Processing with qml file

texts:
I'm doing a spatial join, but the styling doesn't work. What needs to be fixed?

##Style=name
##drivetimepoints=vector
##gridwgs84layer=vector
##DriveTimeGrid=output vector
outputs_QGISJOINATTRIBUTESBYLOCATION_1=processing.runalg('qgis:joinattributesbylocation', gridwgs84layer,drivetimepoints,['intersects'],0.0,0,'sum,mean,min,max,median',0,DriveTimeGrid)
layer = processing.getObject(outputs_QGISJOINATTRIBUTESBYLOCATION_1['OUTPUT'])
layer.loadNamedStyle("C:/QGIS Style File/Drive Time Grids.qml")
layer.triggerRepaint()


texts: Configure Google Maps to display travel time by bike?

texts:
I move around by bicycle but didn't find how to configure Google Maps so that it displays travel time by bike: "x min ride home" would be more useful to me than "x min drive home".

Does Google Maps supports this, and if yes, how to configure it accordingly?


texts: Crediting imagery provider after delineating polygon using satellite imagery

texts:
I would like to create a polygon layer of wetlands for example. If I use satellite imagery to delineate the features, does that mean that the vector polygon is a derivative/modification  of the imagery?

I would like to create and distribute a vector dataset with a very lax licence, but I want to make sure I am truly permitted to distribute it under a highly open license.


texts: Adding attribute data together within multiple polygons in QGIS?

texts:
In QGIS I have two vector layers, the first is a set of polygons and the second is a series of points.  Each point has a numerical attribute associated with it and is located within one of the polygons.  What I'd like to do is to run a points in a polygon function to collect how many points fall within each polygon but in addition, I would like to sum each of the numerical attributes found in each point so that the end result is a total of the numerical attributes for each polygon.

For example: Say I want to generate a map with the total population of each city in the greater Los Angeles region (County). I have point data representing the centroid of each census block group. I then table join the centroid and the population data associated with each block group.  My attribute table for the centroids then has a column listing the population.  I would like to add the population of each centroid together with the population of other centroids that fall into each city's polygon. I've tried the instructions given in Total of attribute field of point layer per polygon layer but I am still am not getting the results that I need.


texts: ArcGIS Server map service layer placed in QGIS?

texts:
I try to load this ArcGIS Server map service layer https://gwa.israntique.org.il/mobidic/rest/services/landuse/MapServer/0 into QGIS 2.18.

When I do it with "ArcGIS REST API Connector" plugin, with the option of "only request feature overlapping the current view extent" the computer stuck:

and I must restart the QGIS, whereas in ArcGIS desktop the layer works fine.


texts: Set expression in QGIS edit widget relation reference with PyQGIS

texts:
I've setup an edit widget for a field as a relation reference like shown below and save the configuration as a QGIS layer style file (*.qml):

With PyQGIS I do the following steps

load the vector layer necessary for the relations
create the relations
load the style files

But after this, the expression dropdown is empty because it is not stored in the qml file:

<edittype widgetv2type="RelationReference" name="nv_bezeichnung">
   <widgetv2config OrderByValue="0" fieldEditable="1" ShowForm="0" Relation="rv-nv" ReadOnly="0" MapIdentification="1" labelOnTop="1" AllowNULL="0"/>
</edittype>

The edit widgets can be configured with PyQGIS like

layer.setEditorWidgetV2Config(fieldId, {'key':'value'})

with keys according to the xml snippet above, but there seems to be no key to store the expression in, therefore my question:

How to store or save the 'Display Expression' in a relation reference or set it up using PyQGIS?


texts: SQL expression error for Select Layer By Attribute

texts:
arcpy.MakeFeatureLayer_management("N05-12_Station2_JGD.shp","station_12")

year=2012

arcpy.SelectLayerByAttribute_management("station_12","#",'"closing">=%d AND "opening"<=%d' % year)

  Traceback (most recent call last):File "", line 1, in

      arcpy.SelectLayerByAttribute_management("station_12","#",'"closing">=%d
  AND "opening"<=%d' % year) TypeError: not enough arguments for format
  string

I want to select records that closing>=2012 and opening<=2012.

Above script, it occurs an error.

How can I write right script?


texts: run a make query table based on a selection

texts: I have a feature class of land parcels. I do a select by location to select the land parcels that intersect a line feature. Now I want to  use the Make Query Table tool in model builder. I include the Shape field to ensure my output is still a feature class however the table is created with all the data and not just the selected features. How do I run a make query table on a selection and not the whole feature class?


texts: ArcGIS not finding newly installed toolbar (.dll, .tlb files?)

texts:
I'm trying to use a toolbar that was developed by the EPA for BMP siting (see below link). The download seems like it should be fairly straightforward, an executable (BMPSitingTool_Setup.msi) is run and the toolbar is downloaded to the Program Files x86 folder within my C-drive. The BMP siting tool comes with a host of .dll, .tlb, .olb, etc files. According to the installation instructions, once I open a new ArcGIS map, it should find this new toolbar automatically and it should show up under Customize > Customize Mode. It doesn't. I've tried doing "add from file" and manually selecting the .tlb file but I get a message of "can't load type library from specified file" followed by "no new objects added". Any idea what's causing this? I'm on ArcGIS 10.3 and the tool was developed within ArcGIS 10.1, but I didn't think there were any compatibility issues with later versions of ArcGIS.

Edit: I've verified that I have administrator privileges on my computer.

https://www.epa.gov/water-research/best-management-practices-bmps-siting-tool


texts: Locking shapefiles from editing by ArcGIS for Desktop?

texts:
How we can lock or protect a shapefile from editing in arcgis?

I want another person to be able to view the data/attributes in the shapefile, but don't want them to be able to change or edit the data.


texts: Mapping simple function to print date and time stamp on Google Earth Engine Python API

texts:
I am new to the Google Earth Engine Python API, and I am trying to create functions that when mapped over a collection will print some metadata (e.g. time stamp, band names) but I always get the same error:

  EEException: ValueNode is empty

Here is the code:

import ee
ee.Authenticate()
ee.Initialize()

# Getting the Sentinel 1 image over the study site
collection = ee.ImageCollection('COPERNICUS/S1_GRD').filterBounds(ee.Geometry.Point(121.1737372,24)).filterDate('2019-01-01','2019-01-30')
collectionList = collection.toList(collection.size())
sentinel1 = ee.Image(collectionList.get(1)).select('VH')

##Function to print the dates of a collection
def printDate(image):
  date = ee.Date(image.get('system:time_start'))
  print('Time stamp:'+str(date.format('Y-M-d').getInfo()))

collection.map(printDate)

def collectionBands(img):
  bands = img.bandNames()
  print('band names:'+str(bands.getInfo()))

collection.map(collectionBands)


texts: Seeking tool like Merge of ArcGIS Desktop in QGIS?

texts:
I'm looking for a tool that might do the same thing as "Merge" in ArcGis. https://desktop.arcgis.com/en/arcmap/10.4/tools/data-management-toolbox/merge.htm
This tool allows you to select several entities, merge them and put them on another table in one manipulation without using another table. Except I need to use it in Qgis.

Do you know any plugin that would do the job?


texts:
When try to preview a layer, got following error related to blocked “HTTP” request:
Mix Content: The page at 'http wms service=WMS&vers..tion/openlayers:321s://dev.geoserver .dmp.azure.cloud.dfo-mpo.gc.ca/geoserver/wms servic 58B1&width=768&hei ght=539&srs=EPSG%3A26713&Forma pplication/openlayers’ was loaded over HTTPS, but requested an insecure frame ‘http://dev.geoserver.dm
p.azure.cloud.dfo-mpo.gc.ca/geoserver /wms ?SERVICE=WM..3659. 3321252958%2C49176
48. 712699675%2C597515 .080901723%2C4921504.461476102". This request has been
blocked; the content must be served over HTTPS.


texts: Calculating average slope along line using ArcGIS Desktop?

texts: I have a slope layer and I want to calculate the average slope along a line. This line isn't necessarily straight (e.g its a route, I would like to be able to input geographic coordinates along that line and have a script that outputs the average slope). Is this possible?


texts: Is there any possibility to make a drop down list of a table in PostGIS using QGIS

texts:
I am trying to create a droplist table filter by id class in PostGIS so that when I start editing in QGIS it will appear as a droplist table with the values or descriptions as I put them in the table.
I can accomplish this task through the QGIS widget attributes form, but this definitely forces me to have the table or tables within table of content, which I do not like it.
Does anyone know how to accomplish this task or has experience with how to accomplish it before?
.


texts: Creating a subclass that inherits from an ArcPy superclass?

texts:
I would like to create a subclass that inherits methods from the arcpy Describe class, as in:

import arcpy

class my_sublass(arcpy.Describe):
    pass

However, arcpy.Describe is a function (that returns a class object, according to ESRI), not a class in and of itself. How can I code my_subclass such that it inherits from the Describe class?

For reference, I am in Python 3, using a cloned environment from Pro's default conda distribution.


texts:
I am using QField v1.8.4 on a Samsung Galaxy Tab S4 running Android 10.
Sometimes when editing an existing feature, the text in a text box widget appears to be duplicated and slightly offset, making it appear blurry. If you then try to edit the text, it appears that the existing text stays and your new text overwrites it. You therefore cannot see properly what you are editing. After saving, the text appears correctly, with edits applied.
Is this a known issue, and if so, is there a solution?


texts: PostGIS st_transform not working as expected

texts:
I am trying to write script about transform coordinates in PyQGIS. I am working on geopackage. When I create geometry with st_transform, coordinates will move right side with same distance. Not transformed geometry's SRID and desired SRID can be same. Here an example:
AsGPB(st_transform(GeomFromEWKT('SRID=5253;POINT(-2217885.89475928433239460 -45973.46579360473697307 0.00000000000000000)'),5253))


texts: PostGIS Street Intersection Autocomplete with full text search

texts:
What I'm attempting to do is create a custom autocomplete for addresses, road intersections, POI's, etc.

For now we are going to using Mapbox GL (we'll be creating a custom map at a later phase), which currently only supports addresses and POI's. Additionally, we will be adding and extracting data from the shape files below.

http://rlisdiscovery.oregonmetro.gov/?action=viewDetail&layerID=556#

http://www.civicapps.org/datasets/street-centerlines

What I've been struggling with is creating a ST_Intersection self join query for the intersecting streets portion of the autocomplete.

The goals include:

1) Street intersection using " / " instead of the customary "&".

2) Partial matching with rank (i.e. "7th Ave & Stark St" would include both SE "7th Ave / SE Stark St" and "SW 7th Ave / SW Stark St").

3) Use of the custom tables and have either TIGER or OSM tables as a fall back for addresses or intersecting streets that outside the area.

I already have the tables created with the appropriate SRID's for each. If needed, I can explode the lines for point data and/or create a new topology edge table, if needed.


texts: Filtered images don't add up - Google Earth Engine

texts:
I came across something bizarre that I can't explain. When I run the following code,

var start = '1982-01-01'
var end = '2018-12-31'
var aoi = ee.Geometry.Polygon([[-92.56227555674951,40.858223918075275],
                               [-92.56227555674951,40.85666590660416],
                               [-92.56227555674951,40.85666590660416],
                               [-92.55961480540674,40.85828883442468]]);

var chirps7a = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY")
                .filterDate(start,end)
                .filter(ee.Filter.dayOfYear(1,9))
                .filterBounds(aoi)
    print(chirps7a.size())
//returns 333 values, an odd number

    var chirps7b = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY")
                .filterDate(start,end)
                .filter(ee.Filter.dayOfYear(1,10))
                .filterBounds(aoi)
    print(chirps7b.size())
//returns 370 values, an even number, which makes sense considering the timespan is 36 years.

selecting 9 days returns 37 years worth of data, while having 10 days returns the desired 36. I'm sure there's something simple I've missed, but can someone point it out?


texts: Standalone QGIS App keeps crashing. Issue with Context or Threads?

texts:
I'm trying to setup an application for my company, that generates a PDF file for specific user generated coordinates:

The user enters an address into the frontend (Flask generated), which is translated via Here Maps Geocoder API to coordinates.
The coordinates are sent to pyQGIS running on the company server.
PyQGIS connects to PostgreSQL, performing some spatial joins. This way suitable information for the coordinates is collected in a PDF file being returned to the user.

The good thing: It works! Except for when loading QGIS Command Prompt says:

libpng warning: iCCP: known incorrect sRGB profile
QObject::connect: invalid null parameter

But it seems to be no big issue, it works anyway.
The bad thing: It only works when all the PyQGIS code is in one single script, looking like this (shortened, loading only the PDF cover):
def print_pdf(new_address_dict):
    x = new_address_dict['x']
    y = new_address_dict['y']
    address = new_address_dict['address']

    path = os.path.dirname(os.path.abspath(__file__))

    QgsApplication.setPrefixPath("C:/Program Files/QGIS 3.4", False)

    # Create a reference to the QgsApplication.  Setting the
    # second argument to False disables the GUI.
    qgs = QgsApplication([], False)

    # Load providers
    qgs.initQgis()
    print('QGIS initialisiert, Brudi')

    # Load project
    project = QgsProject.instance()
    composition = QgsPrintLayout(project)
    document = QDomDocument()
    project.read(path + '/pyqgis/dynamischer_test.qgz')

    # Load layout
    layout = project.layoutManager().layoutByName("dynamisches_layout")

    # Set Headline
    headline = layout.itemById('00 Titel')
    headline.setText(address)
    headline_shadow = layout.itemById('00 Schatten')
    headline_shadow.setText(address)

    # Set extent
    refmap = layout.itemById('00 HG')
    xmin = x - 250
    ymin = y - 250
    xmax = x + 250
    ymax = y + 250
    bb = QgsRectangle(xmin, ymin, xmax, ymax)
    refmap.setExtent(bb)

    # Set background
    hintergrund = project.mapLayersByName('gruen')[0]
    project.layerTreeRoot().findLayer(hintergrund.id()).setItemVisibilityChecked(True)
    print('Cover done')

    # Print pdf
    pdf_path = os.path.join(path + '/pyqgis/', 'test.pdf')
    exporter = QgsLayoutExporter(layout)
    image_settings = QgsLayoutExporter.PdfExportSettings()
    image_settings.dpi = 300
    result = exporter.exportToPdf(pdf_path, image_settings)

    # Exit QGIS
    qgs.exitQgis()
    print('Done!')

This isn't very clean and also QGIS is being loaded every time the user presses the 'Generate PDF'-button, which ist very unefficient. So I tried to move all the startup code, like creating a QgsApplication reference, loading the project etc, into an own method, which is called right at the start of the script.
Also I wanted to have an own method for each of the pages in the PDF, like load_cover(), load_page1(), etc. I passed the project and layout variables into these new methods like this:
def load_cover(project, layout, address, x, y):
    headline = layout.itemById('00 Titel')
    headline.setText(address)
...

But when I things this way, there seems to be some issue with the context or threads or something. Python crashes at the very last line:
result = exporter.exportToPdf(pdf_path, image_settings)
And Command Prompt says the following:

QPrinter: Must construct a QCoreApplication before a QPrinter

I'm new to python and PyQGIS and I'm stuck with this. Please help! I'm using Python 3.8 and QGIS 3.4.


texts: Adding LAS files into AutoCAD map 3D with Point Cloud Manager tool

texts:
It's my first time to use AutoCAD map 3D for LiDAR. Based on my googling, it is said that I can bring .las file by creating new index and using the Point Cloud Manager tool (the image on the top - reference).

However, my AutoCAD does not show the same interface as shown in the video (the image at the bottom). I was not sure because it is a trial version. Any advice?


texts: libLAS georeference las file and update header

texts:
I am just starting to use libLAS and I had a small issue.

I have next code to generate a las file with georeferenced points:

std::ofstream ofs;
liblas::Header header;
liblas::SpatialReference srs;
srs.SetFromUserInput("EPSG:4326");
header.SetSRS(srs);
header.SetPointRecordsCount(1);
ofs.open("H:\\temp\\test2.las", std::ios::out | std::ios::binary);
liblas::Writer writer(ofs, header);

liblas::Point point(&header);
point.SetCoordinates(1.0, 2.0, 3.0);
writer.WritePoint(point);

point.SetCoordinates(5, 5, 5);
writer.WritePoint(point);

point.SetCoordinates(-1, -1, -1);
writer.WritePoint(point);

header.SetPointRecordsCount(3);

writer.SetHeader(header);

writer.WriteHeader();
ofs.Close();

When I run lasinfo over this file I get: It says that I have 3 points, but my spatial reference is none.

Same code above, but if I comment this line:

//writer.SetHeader(header);

I can see the spatial reference, but the number of points is 1.

Am I doing something wrong?


texts: Extract boundary cells between classes from raster

texts:
I want to extract from a RasterLayer all boundary cells between classes (possibly according to the 4-cell rule), excluding from this operation cells that neighbour NA cells

The raster::boundaries function detects boundaries as cells that have more than one class (including NA) cells in the 4 or 8 cells surrounding it , or, if classes argument is set to FALSE, cells with values and cells with NA.

What about if I do not want to consider as edge the NA border?


texts: How do I programmatically update/recreate an ArcSDE feature class from external sql data?

texts: I have an ArcSDE feature class that contains project points based on XY data from an external SQL database. (Not a geodb). I am looking for the best way to "recreate" the ArcSDE feature class from the external data on a daily basis. The external database is not not connected to ArcSDE and projects may be added or removed from the external database, which needs to be reflected in the ArcSDE feature class.


texts: Calculate min distance between points in PostGIS

texts:
I have point layer in PostGIS database. Help please with query which add for each points in layer distance in attribute table between this point and another nearest point. Thanks!

Update
Calculate distances between series of points in postgis there is similar question, but
there is in answer:

UPDATE my_table SET dist=ST_Distance(my_table.geom, b.geom)
FROM (SELECT geom FROM my_table WHERE **gid = 1**) b;

distance added for points - from gid = 1 point to another points.
For me each point object in attribute table should have distance from point to another nearest point.

Update

For example, I need add distance between point and another nearest point in attribute table for point:
  for point 1 - distance between 1 and 2;
  for point 2 - distance between 2 and 3;
  for point 3 - distance between 3 and 2;
  for point 4 - distance between 4 and 3;
  for point 5 - distance between 5 and 4;
  for point 6 - distance between 6 and 3.


texts: Replacing part of raster with another one using ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10.5 and I am unable to find how to replace a part of a raster with another raster.

I have a Raster for topography and for a part of it I need to recreate a past morphology which was done by creating a new raster from isoline but only for a part of it as briefly shown here:

How can I do it?


texts: How to cut unwanted area (polygon) using ArcGIS for Desktop?

texts:
I am using ArcGIS 10.2 for Desktop with a Basic level license.

There are two overlapping polygon. Polygon A (with line stripes) and Polygon B (in beige color). Polygon A and Polygon B is a separate feature class. Both feature class contains similar attribute column details.

Is it possible to remove all area in A and get the remaining area B?

I've been trying to use Union but it does not works.


texts: WMS web service - blocking users with wrong referrer?

texts:
I would like to to use the WMS url from https://www.alegrodialog.de/profiles/markaspot/modules/mark_a_spot/modules/markaspot_wms/markaspot_wms.js which is about the ALEGrO DC power line.

var alegro = L.tileLayer.wms("https://map-proxy.zebralog.de/maps.zebralog.de/cgi-bin/maps-amprion", {
  layers: 'zebralog',
  service: 'wms',
  version: '1.1.1',
  format: 'image/png',
  srs:"EPSG:3857",
  attribution: "Alegro Trasse",
  transparent: true
});

Shouldn't it simply be https://map-proxy.zebralog.de/maps.zebralog.de/cgi-bin/maps-amprion ?
Or more elaborate:

  wms:https://map-proxy.zebralog.de/maps.zebralog.de/cgi-bin/maps-amprion?FORMAT=null&VERSION=1.1.1&SERVICE=WMS&REQUEST=GetMap&LAYERS=zebralog&STYLES=&SRS={proj}&WIDTH={width}&HEIGHT={height}&BBOX={bbox}

But I've tried with QGIS and JOSM and the URL is dead.

Have I done a mistake or is the WMS endpoint really blocking me?


texts:
I'd like to merge linestrings which touch, but do not cross. The linestrings share the same SRID and come from a PostgreSQL 9.2/PostGIS 1.5 database. Within the table I query are several line strings, not just those I want to merge, but also others.

For example, consider lines from the OpenStreetMap database. I start with a temporary table that has all linestring which are tagged as being part of the power grid, together with the respective voltages. However, not all parts of the grid that belong together are also represented as a single (multi-) linestring, they are within the database the way they were inserted. So, I need to find out which lines touch, do not cross, and share the same voltages, and finally merge them.  The result will be another, shorter list containing merged lines and those which did not need to get merged with others.

I've come up with an algorithm in C++ (Qt), which works, but is slow:

for (QHash<qlonglong, MapItem *>::iterator current = gridItems.begin();
        current != gridItems.end(); current++) {
    for (QHash<qlonglong, MapItem *>::iterator other = current;
            other != gridItems.end(); other++) {
        if (other.value() == current.value()) {
            continue;
        }

        const OGRGeometry *firstGeom = current.value()->geometry();
        const OGRGeometry *itGeom = other.value()->geometry();

        if (current.value()->tag("voltages") == other.value()->tag("voltages")
                && firstGeom->Touches(itGeom)
                && !firstGeom->Crosses(itGeom)) {

            // Merge:

            MapItem *mergedItem = new MapItem(firstGeom->Union(itGeom));
            mergedItem->tags(current.value()->tags());

            *current    = mergedItem;
            *other      = mergedItem;
        }
    }
}

Basically, it does the following: Given a list of pointers to geometries, it creates two iterators within two loops. The second iterator always starts as a copy of the first; then both travel until the end of the list. If the geometries the two pointers point to touch, but do not cross (and also share some tags), they are united.

To illustrate what it does, I've drawn a very simple picture featuring five line strings which could be a possible result of my example query outlined above. In this picture, only the three grey ones are merged:

I now want to improve this approach by (a) possibly improving the algorithm and (b) putting the strain on the database server, not the slower client machine.

I've two options here, but I do not know which one is the best or which actually works:

Create a SQL query which does this. However, I do not know how. (Self-join? Aggregates?) If someone could suggest a query that does this, I'd be absolutely happy.
Create a PostgreSQL C extension that implements my algorithm. If so, I'm happy for any suggestions of how to improve the algorithm or improve it by enabling parallel processing.


texts: What expression should be used to return the total area of all polygons in a layer?

texts:
I'm presuming this must be possible, which may be an error, but I need an expression in Qgis 2.16 which calculates the total area of the polygons in a layer. I intend to use it to give a proportion of cover for each polygon in the layer.

I found this which seems to claim to create such an expression in the editor, but it doesn't work as expected. The preview in the field calculator shows the right answer, but the column in the table is populated with zeros and when I next open QGIS and then open the table it crashes.

I'm not interested in Statist or Groupstats (where similar questions seem to lead) because I need to use the value in my field calculator.

If anyone knows of an expression that would be perfect, but otherwise, can anyone could have a look at the script I've linked and see why it might no longer work (presuming it did once) and suggest any updates?


texts: Add Feature Class from Geodatabase Dataset to MXD using ArcPy Standalone Script?

texts:
I am trying to add a Feature Class from a Geodatabase Dataset to an MXD and save a copy like this but I am getting following errors on each methos

import arcpy
import os

arcpy.env.workspace = r"C:\arcgis\ArcTutor\AAA\MXDs"
MXD = arcpy.mapping.MapDocument(r"C:\arcgis\ArcTutor\AAA\MXDs\map.mxd")
DF = arcpy.mapping.ListDataFrames(MXD)[0]

lyr = arcpy.mapping.Layer("C:\arcgis\ArcTutor\AAA\GDBs\Delaware.gdb\Kent\Kent")
arcpy.mapping.AddLayer(DF, layer, "AUTO_ARRANGE")
MXD.saveACopy('newmxd.mxd')
del MXD
print(" Done")

and

import arcpy
import os
arcpy.env.workspace = r"C:\arcgis\ArcTutor\AAA\MXDs"
MXD = arcpy.mapping.MapDocument(r"C:\arcgis\ArcTutor\AAA\MXDs\map.mxd")
DF = arcpy.mapping.ListDataFrames(MXD)[0]

lyr = "C:\arcgis\ArcTutor\AAA\GDBs\Delaware.gdb\Kent\Kent"

arcpy.MakeFeatureLayer_management(lyr, 'shplyr')
arcpy.mapping.AddLayer(DF, lyr, "AUTO_ARRANGE")
print(" Step 4")
MXD.saveACopy('newmxd.mxd')
del MXD
print(" Done")


texts: Error "StopIteration: iteration not started" when updating row

texts:
I get the following error when trying to update the row on a for loop where I have assigned values to the row's attributes:

Traceback (most recent call last):
  File "D:\ArcGIS_Pro_Projects\Optimity_Alf\Arcpy_Scripts\LOS_Many_Points.py", line 39, in <module>
    cursor.updateRow(row)

StopIteration: iteration not started
Failed to execute (LOSManyPoints).

See below the code:

#Import geoprocessing.
import arcpy

#Import workspace.
arcpy.env.workspace = r'D:\ArcGIS_Pro_Projects\Optimity_Alf'

#Input parameters
buffer_distance = arcpy.GetParameterAsText(0)

#Input variables.
fc_Opt_All = r'D:\ArcGIS_Pro_Projects\Optimity_Alf\Optimity_Equipment.gdb\Optimity_All'
fc_Opt_Opp = r'D:\ArcGIS_Pro_Projects\Optimity_Alf\Optimity_Equipment.gdb\Optimity_Prospect'
fc_London_Bldngs = r'D:\ArcGIS_Pro_Projects\Optimity_Alf\London_Building.gdb\London_Buildings_Multipatch'
DSM_model = r'D:\Data4Optimimty\ESRI_UK_Bluesky_DSM\ESRI_UK_Bluesky_DSM.gdb\DSM_25cm_2015'

arcpy.MakeFeatureLayer_management(fc_Opt_Opp, 'fl_Opt_Opp')
field_names = ['LOCATION', 'ESTATE', 'ADDRESS', 'OFFSETA', 'OFFSETB', 'Z', 'OFFSETC', 'AZIMUTH1', 'AZIMUTH2', 'VERT1', 'VERT2', 'RADIUS1', 'RADIUS2', 'STATUS', 'CONNECTION']
observers = []
targets = []
#Calculating LOS for each point
with arcpy.da.UpdateCursor(fc_Opt_Opp,field_names) as cursor:
    for row in cursor:

        #Adding default values to observers' attributes
        row[1] = 'Prospect'
        row[3] = 1
        row[4] = 0
        arcpy.AddSurfaceInformation_3d('fl_Opt_Opp', DSM_model, 'Z', 'BILINEAR')
        row[6] = row[3] + row[5]
        row[7] = 0
        row[8] = 360
        row[9] = 1
        row[10] = -16
        row[11] = 0
        row[12] = 2000
        cursor.updateRow(row)


texts:
Exporting the results of a shortest path analysis through the Road Graph plugin in QGIS 2.0.1 results in an "empty" layer.

After running teh shortest path calculator, I selected export --> new temporary layer.  I select the coordinate system (WGS84 - UTM18N)

Note in the photo, that a path on my roads vector layer is selected (red) by Road Graph, but the attribute table is empty.  When I save this temporary layer permanently, the layer remains empty as well.

Please send instructions on how to make the export work/another way to save shortest paths as separate shapefiles.


texts: How to automatically place parcels on suitability map?

texts:
I need to make a map that randomly places parcels (groups of pixels or polygons) on a suitability map by first filling the most suitable areas, then the next most suitable and so on, until a pre-defined number of parcels have been placed on the map.

This is probably best explained with an example.  I have a shapefile with polygons on it (counties) and each polygon has a value (number of acres for a crop).  I also have a suitability map (0-9 suitability for that crop).  I need to place the number of acres defined in each polygon onto the suitably map.  I want to first fill the most suitable lands (9) and if there is still land left fill the next most suitable land (8) and so on until the number of acres for that county are distributed.  Each polygon has a different number of acres and a different range of suitability, e.g. some will have 0-9, others might only have 3-7.  Within each suitability layer the placing of the parcels can be random.

I have just started looking into Python scripting for ArcGIS.  It seems like I might be able to do this with Python but I'm not really sure where to start.  I would be interested in hearing any suggestion or advice on how to accomplish this.


texts:
I am trying to download raw radar data off of https://s3.amazonaws.com/noaa-nexrad-level2/index.html

The files do not have extensions. What format is the data in (grib, netcdf, etc.)? How do I extract it using 7-zip? How do I view it as an image?

I know that there are Python programs, but I was having trouble installing the cartopy module.

https://unidata.github.io/MetPy/latest/examples/formats/NEXRAD_Level_3_File.html#sphx-glr-examples-formats-nexrad-level-3-file-py


texts:
Currently, I'm mapping data from the Hudson's Bay Company Archives in the nineteenth century, and I am tracing the routes travelled by fur traders. I have gone through the documents, recorded the start and end points of each trip, and the decimal coordinates of each. I created point data on Google Earth, exported the data as .kml files, and imported it all into QGIS. My current map looks like this:

I would like to have lines drawn between points on the map indicating the travels made by traders. Ideally, there would be a line for each voyage to a specific point, meaning heavily travelled routes are clearly visible to the observer.


texts: Keeping properties in an Intersection - Google Earth Engine

texts:
I have a set of points and I need to classify them according to which polygon they belong. The natural way of doing this is creating the intersections between all points and each polygon. The points have properties I want to keep and the intersection function requires me to have the polygon on the left side and I lose the properties of the input of the right side. How should I go about this? If I use the points on the left side it doesn't work.

https://code.earthengine.google.com/f7b34e63c5e9ad70f96e1ad91620d48a

I want an output that I can download as a .csv with each point, it's characteristics, and an indicator of the polygon they belong.

Code:

var p1 = ee.Feature(p1, {id: 'id1'});
var p2 = ee.Feature(p2, {id: 'id2'});
var p3 = ee.Feature(p3, {id: 'id3'});
var p4 = ee.Feature(p4, {id: 'id4'});

var points=ee.FeatureCollection([p1,p2,p3,p4]);
print(points,' points');

var cc1 = poly1.intersection(points, ee.ErrorMargin(1));
var cc2 = poly2.intersection(points, ee.ErrorMargin(1));

Map.addLayer(cc1,{},'cc1');
Map.addLayer(cc2,{},'cc2');

var ccs=ee.FeatureCollection([cc1,cc2]);
print(ccs,'ccs');


texts: Export to PDF with Nepali (Preeti, Kalimati) fonts in ArcGIS 10.2.2

texts: I am using ArcGIS 10.2.2, created a map with Nepali (Preeti, Kalimati) fonts. Fonts not correctly displayed in exported PDF file. Any solution?


texts: Automatically snapping points to line ends on two separate layers?

texts:
I have US interstate ramp end location data points that I need to match road network data.  All the points are off by a few meters as they were collected in the field.  I am hoping there is a way to automatically match/snap all the points to the end nodes of the ramps if within 20 meters.  I can't do it manually as there are over 30,000 points.  Road network and ramp points are different layers.

EDIT: The blue lines are the on/off ramps of an interstate. The red points are the field collected location of the end of the off ramp.  I need the red points to auto snap on to end of the off ramp.

After thinking about this and the goal all I really need is the start and end points for each ramp and a way to differentiate between them.  I'll ask in a different question.  However for clarity sake I thought I would post what I originally was requesting.

I haven't found any information on how to automatically snap points to line ends in bulk.  Anything I've found that was similar didn't translate.


texts: How can you tell if a RDBMS is an ArcSDE Geodatabase?

texts:
I want to know if there is any way I can query the RDBMS to determine whether it is an ArcSDE Geodatabase, or if I have to do it through ArcObjects. I am wondering if it would be enough to check if the SDE schema exists on the target database, and then maybe check if the SDE.VERSION table also exists and then conclude that it's an ArcSDE Geodatabase.

The ArcSDE Geodatabases are Oracle, and I use ArcObjects SDK 10 .NET in C#.


texts:
We have a number of custom local coordinate systems at work.

I'm trying to create a well known text coordinate system for each of them however I'm stuck on how to define these systems they .

Some variables i'm aware of. centre of latitude longitdue -27 22 27.1954, 153 10 46.5722 (lat,long) GDA94 Zone56 centre 517760.037, 6972102.441 (E,N).

We set this coordinate to 3800.192,3129.572 then rotate the file about it by -57.975833 and scale it by 1.000400336 (this is a cad move rotate scale operation not a re-projection as such)

The only thing I've come close to getting a similar result has been a stereographic coordinate system ... but this scales in all directions? We require a cylinder type operation as our area of work is more rectangular and it will give greater accuracy along the centre of work. Can anyone help me define a WKT coordinate system for my proect?

updated wkt code ... im not sure on how to define the rest of the coordinate system help?

FITTED_CS["FIG",
      PARAM_MT["Affine",
               PARAMETER["elt_0_0", 3800.192],
                                 PARAMETER["elt_0_1", 3129.572],
                                 PARAMETER["elt_1_0", -0.53027691844514246],
                                 PARAMETER["elt_1_1", 0.84782450410702548],
                                 PARAMETER["elt_2_0", 0.84782450410702548],
                                 PARAMETER["elt_2_1", 0.53027691844514246]
                                 <need long lat of swing point somewhere?>
              ],

what is exported from my software ...

PROJCS["FIG",
 GEOGCS["FIG",
  DATUM["D_EPSG:6283",
     SPHEROID["EPSG_7019",6378137,298.257222096042]
    ],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]
 ],
 PROJECTION["CSMAP:54"],
 PARAMETER["Central_Meridian",0],
 PARAMETER["Latitude_Of_Origin",-27.3742209444444],
 PARAMETER["Scale_Factor",1.000400336],
 PARAMETER["False_Easting",0],
 PARAMETER["False_Northing",0],
 PARAMETER["A0",3800.192],
 PARAMETER["A1",-0.53027691844514246],
 PARAMETER["A2",0.84782450410702548],
 PARAMETER["B0",3129.572],
 PARAMETER["B1",0.84782450410702548],
 PARAMETER["B2",0.53027691844514246],
 UNIT["Meter",1]
]


texts: Automatically digitize polygons by most outlying points within cluster in ArcMap

texts:
I'm looking for a tool or script that will allow me to automate the digitization of polygons by the most outlying points in a cluster of points with the same key attribute. Below is a sample image displaying this done manually.

Selected points are all with the same ID of "MD_1001". I have drawn a polygon around these points, and angled the polygon via the most outlying points in the cluster.

I have thousands of these clusters to digitize polygons over, and need some sort of tool or script to get this done quickly. I am using ArcMap 10.2.2.


texts: Show updated attribute values during export animation (QGIS Temporal Controller)

texts:
I have a vector layer that shows the values for a specific timestep in a NetCDF.
If the user changes the timeframe, values for that timestep are read from the NetCDF and written to the vector layer. I do this by connecting the read/write method to the Temporal Controllers' updateTemporalRange signal.
This works fine when interacting directly with the temporal controller.
But if I export the animation, the vector layer's data is not update between frames, so I just get a series of identical output images.
Is there a way to force the animation exporter to read the updated attribute values before each render?


texts: Copying some entries from one column to another with QGIS field calculator?

texts: I have some rows in the AttributeTable where I want to copy some entries from one column to another. I thought that the best way to do this was in the field calculator and use a CASE WHEN expression followed by IN (and the references I want to update) but I'm not sure where to go next? I thought I'd use THEN to select column I want the information from but it doesn't seem to work?


texts: Preserving geometry of Union output with Virtual Layer in QGIS 2

texts:
There is a polygon layer "Union23" that was created with 'Union' tool.

Layer's attribute table includes the following values.

My effort is simply to aggregate those polygons by "FLAECHEID" field with the sum and count of a field "In_Value1" via Virtual Layer. So I am using the expression

SELECT *, SUM(In_Value1), COUNT(In_Value1)
FROM Union23
GROUP BY "FLAECHEID"

However, I am getting really weird output and I do not understand why I am losing a part of geometry on the way? Any suggestions?

When I use the solution offered by @JGH I am somehow getting better output but still with small issues. It is most like related to the original geometry. Improving geometry is necessary as was mentioned by @Stu Smith.

References:

Select all columns except geometry using virtual layers
Unable to show geometry in virtual layer qgis


texts: Digitizing two layers of same table with different projection Fails in QGIS?

texts:
When I try to digitize two layers of different projection coming from the same table fails to digitize.  I am using a PostGIS table as a data source. The table has two geometry columns with different projection, which I want to digitize. What happens is the following:

When I add layer 1, it works in digitizing.
Then when I add the second one while the first one is there, the second one work.
But when I get back to layer 1, it draws an invalid geometry.

Unfortunately, after that, key things do not work, especially featureAdded and beforeCommitChanges signal and layer.rollBack().
I am using a custom form by supressing the default one using curr_layer.editFormConfig().setSuppress(1) this stops working. When I save data , I use beforeCommitChanges to do the saving directly using SQLAlchemy. Then I excute the undo stack so that QGis thinks there is nothing to save. However, it ignores beforeCommitChanges also try to save the same feature.

It there any known issue in this regard and a solution?


texts: I want to create a workflow in ModelBuilder. The workflow I work on is this: first I create a fishnet on my land use raster's extent. Then I assign dominant land use classes on every grid as table with majority statistic type in "Zonal Statistics as table". Then I want to join the output table coming from zonal statistics into the empty fishnet (Figures Below). But I could not join together these two different features with "Add Join" command.


texts:
I have received a drawing with a type of projection system which I do not recognize.

The drawing says that it's LO 19/WGS 84.

I would like to convert this to decimal degrees to be able to plot the drawing in Arcmap.