texts: Leaflet Zoom on WMS-Layers - Error: Attempted to load an infinite number of tiles GridLayer.js

texts:
I use Leaflet 1.5.1 for showing WMS layers on a GeoServer 2.15.1. I can only zoom up to zoom level 14. If I continue to zoom into the map, I see the following error message.

Error: Attempted to load an infinite number of tiles GridLayer.js:658:44

I found this: https://github.com/Leaflet/Leaflet/issues/6416
But I am sure, that is not the problem for me. I set zoom correcly.

Do I have to pay attention to zooming WMS layers?

Edit:
I noticed the mistake myself.

I defined a CRS

    var crs = new L.Proj.CRS('EPSG:25833',
        '+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs',
        {
            resolutions: [
                8192, 4096,
                2048, 1024, 512, 256,
                128, 64, 32, 16,
                8, 4, 2, 1, 0.5,
                0.25, 0.125, 0.0625, 0.03125,
                0.015625, 0.0078125
            ],
            origin: [0, 0]
        })

and I made a mistake. I used in the last two resolutions a comma as the decimal separator instead of the dot. The entry 0,015625, 0,0078125 instead of 0.015625, 0.0078125 leads to the error message.


texts: Extracting blank space between buffers regions

texts:
I have a set of parallel lines with standard distance of 1.4m and I have created a buffer of 0.70cm for each row of the lines. Now I need to identify the regions where the buffers do not touch, which is where the distance between the lines is greater than 1.4m.
Image for example,  I need to save the blanks between the buffers:

I created a convex hull of shp and then used the symmetric difference to get the blanks, the problem is that all spaces are united in one feature, I want to separate them.

I do not know if there is another way to do it, any suggestions are welcome! I use QGIS and I'm learning postgis.


texts: How to set up 'Extent' environment variable with arcpy when using CellStatistics?

texts:
I have a set of rasters of different extent and I have written a script to combine them all into one raster using CellStatistics:

import arcpy
from arcpy.sa import *

arcpy.env.workspace = r'C:\ESRI\ArcGIS_Pro_Projects\DATA\Optimity_Alf\Work.gdb'

#Path where rasters are located
path = r'C:/ESRI/ArcGIS_Pro_Projects/Optimity_Alf/DATA/Work.gdb'

#List of raster names from workspace
rasters = arcpy.ListRasters('*', 'All')

#Building list of raster complete paths
paths = []
for raster in rasters:
    if raster[0] != '_':
        complete_path = path + "/" + raster
        paths.append(complete_path)

#Combining rasters into one unique raster
path_combined_raster = r'C:\ESRI\ArcGIS_Pro_Projects\Optimity_Alf\DATA\Work.gdb\Combined_Raster'
combined_raster = CellStatistics (paths, 'MINIMUM', 'DATA')
combined_raster.save(path_combined_raster)

However, when I run this script it throws the following error:

    Traceback (most recent call last):
  File "C:\ESRI\ArcGIS_Pro_Projects\Optimity_Alf\Arcpy_Scripts\Combining_Viewsheds.py", line 21, in <module>
    combined_raster = CellStatistics (paths, 'MINIMUM', 'DATA')
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\sa\Functions.py", line 3173, in CellStatistics
    ignore_nodata)
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\sa\Functions.py", line 3169, in Wrapper
    [function] + Utils.flattenLists(in_rasters_or_constants))
RuntimeError:  ERROR 999998: Unexpected Error.
 Failed to execute (CombinedViewsheds).

I have no idea what all that error message means. However, I tried to do the same operation in ArcGIS Pro and got an error regarding the extent of the rasters. Therefore I went to Environments and set the extent to 'Union of Inputs'. That change made the thing work fine. My question is, how can I set the extent in arcpy? The instruction "CellStatistics (paths, 'MINIMUM', 'DATA')" does not accept more parameters.


texts: Can the field mapping be shortened in ArcToolbox tools?

texts:
I have converted a model to python in ArcGIS 10.3, and have found that the tool outputs a full and very very long field mapping.  Is there a way to shorten this somehow?  While I want all field names/types/lengths to stay the same, I don't necessarily want all fields to be output.

arcpy.TableToTable_conversion(SourceTable, TargetDB, "xxOutputTable", "", "RecordID \"RecordID\" true false false 8 Double 0 19 ,First,#,Database Connections\\\\GISDatabase.sde\\GISDatabase.dbo.v_InputTable,RecordID,-1,-1;Descript \"Descript\" true true false 65 Text 0 0 ,First,#,Database Connections\\\\GISDatabase.sde\\GISDatabase.dbo.v_InputTable,Descript,-1,-1;Status \"Status\" true true false 30 Text 0 0 ,First,#,Database Connections\\\\GISDatabase.sde\\GISDatabase.dbo.v_InputTable,Status,-1,-1;RecordName \"RecordName\" true true false 50 Text 0 0 ,First,#,Database Connections\\\\GISDatabase.sde\\GISDatabase.dbo.v_InputTable,RecordName,-1,-1;Owner \"Owner\" true true false 2000 Text 0 0 ,First,#,Database Connections\\\\GISDatabase.sde\\GISDatabase.dbo.v_InputTable,Owner,-1,-1;", "")

While the python works fine with it all there, I would like to be able to make it as simple as possible to modify if/when required, and I do find the long field mapping quite hard to read when scrolling through it across the screen.

Is there a way to simplify this expression?


texts: Why would a perfectly usable raster mosaic produce an unusable clipped raster?

texts:
I'm getting funky results from Extract By Mask (spatial analyst).  I'm clipping a Raster Mosaic of DEMs from the National Elevation Dataset with polygons.  I have a script to run through 15 separate polygons, however, after getting strange results I've resorted to trying to get a single mask to run.

I'm getting a warning that there is a Datum conflict, but I was careful to make sure to use the same Datum...

My inputs:

My output:


texts: Select and copy features in ArcMap using Python add-in tool

texts:
I'm trying to add a tool to an add-in toobar to select features from an existing feature class and copy them across to another feature class. I need two tools, one for a point and one for a rectangle. Add-in tools only handle rectangles, so I'm trying to generate the point coordinate off the onMouseDown event and simply use the X and Y feedback I get. The problem is that I can't seem to pass these coordinates to anything usefull to extract the features. I've tried using Environment extents and then doing a simple CopyFeatures, but that's not working as it copies the entire fc and doesn't seem to honour the environment settings I set.

Can I pass coordinates to the Select by Location tool somehow or is there another way of passing coordinates to something to extract by that extent?

This is what I have now:

import arcpy
import pythonaddins
arcpy.overWriteOutput = True

# Replace this with the SDE layer once it's available.
global ELAtemplate, fc
ELAtemplate = r'C:\Data\nsw_map_units.shp'
fc = ""

class DefineUnitsbyPoint(object):
    """Implementation for DefineUnitsbyPoint.tool (Tool)"""
    def __init__(self):
        self.enabled = True
        self.shape = "Rectangle" # Use onMouseDown to get initial extent of the rectangle.
    def onMouseDownMap(self, x, y, button, shift):
        # fc can be assigned from a combo box selection in another class
        global ELAtemplate, fc
        if fc == "":
            pythonaddins.MessageBox('Choose a layer to select from.', 'Choose a Layer', 0)
        else:
            mxd = arcpy.mapping.MapDocument("CURRENT")
            pointGeom = arcpy.PointGeometry(arcpy.Point(x,y), mxd.activeDataFrame.spatialReference)
            arcpy.SelectLayerByLocation_management(ELAtemplate, "INTERSECT", pointGeom, "", "ADD_TO_SELECTION")
            arcpy.CopyFeatures(ELAtemplate, fc)


texts: PyQGIS problem to find the path of the Python script

texts:
I have two Python scripts. The main script import the second script
test_credenziali.py (main script)
#!/usr/bin/env python
# coding=utf-8
# Copyleft Roberto Marzocchi - Gter srl Innovazione in Geomatica Gnss e Gis
import os,sys,shutil,re,glob, getopt
dirname0 =os.path.dirname(sys.argv[0])

#dirname = os.path.dirname(os.path.realpath(__file__)) #NameError: name '__file__' is not defined

dirname2= os.getcwd()
print('The path of the python script is ', dirname0 )
#print('The path of the python script is ', dirname )
print('The path of the python script is ', dirname2 )

# here i like to use the dirname variable
sys.path.insert(0, r'C:\Users\assis\Desktop\script_test')
from credenziali import *
print(host)

credenziali.py
host='domain.it'

I save both script in the same folder (e.g on my desktop) C:\Users\assis\Desktop\script_test
and if I use the following line before the import of credenziali.py it works
sys.path.insert(0, r'C:\Users\assis\Desktop\script_test')
Nevertheless I like to use relative path in order to call the credenziali.py script. Using python normally it works or I can use the command
dirname = os.path.dirname(os.path.realpath(__file__))
to obtain the path of the main script.
Using the QGIS Python console it doesn't work.


texts: Continents geometry data / extracting contours of continents from country geometry

texts:
I have an ESRI Shapefile data from Natural Earth (Admin-0 Countries). I converted it to geoJSON format. There is no world or continents geometry there, only the county one. And I'd like to extract eg. only Europe contours, without any country borders. And there's another problem - what about Russia? I'd like to include Europe-Asia border in Ural mountains, and there's no such data in my files.

Now I'm thinking that maybe the whole idea of converting that NE data is wrong - I bet that somebody already done that and there's world & continent geometry available to download somewhere on the internet - but I can't find it...


texts:
Image a polygon layer with fields calculated from values found in a point layer. No witchcraft, there, as that's what we have all been doing for years.
The polygon layer needs to be populated with counts and aggregates coming from the point layer. This is an operation I will repeat many times during my work.
Here are some of the operations I will want to include. There are more like these.
Sum myNET
aggregate( layer:= 'MyPointLayer_aa66a666_a66a_666a_a6ab_6a6a6aa666', aggregate:= 'sum', expression:= myNET, filter:= contains(geometry(@parent), $geometry))

Sum myNET if myHZyear older than 2010, exclude WP in myHZtype
aggregate( layer:= 'MyPointLayer_aa66a666_a66a_666a_a6ab_6a6a6aa666', aggregate:= 'sum', expression:= myNET, filter:= contains(geometry(@parent), $geometry) AND ("myHZyear" < 2010 OR "myHZyear" is null) AND "myHZtype" != 7206)

Sum myNET if myHZyear older than 2010, exclude WP and FW in myHZtype
aggregate( layer:= 'MyPointLayer_aa66a666_a66a_666a_a6ab_6a6a6aa666', aggregate:= 'sum', expression:= myNET, filter:= contains(geometry(@parent), $geometry) AND ("myHZyear" < 2010 OR "myHZyear" is null) AND "myHZtype" != 7206 and "myHZtype" != 7208)

myNETperArea
myNET / $area / 10000

Instead of doing a lot of clicking each time I have a point layer and a polygon layer and I need to do these calculations betwixt them, I would much prefer to store the calculations in one single operation. I tried in graphical modeler, but that produces one output layer per edited field and also a bunch of new layers.
Is it possible to chain field calculator output so that the fields of the input layer are calculated in one go?
I have started to make a model with the graphical modeler with also some steps before what is described here and have also exported that model to a python script. If needed, I can of course dump the whole script here, but wouldn't want to do so right away if the question is clear enough as is.
Some questions I have considered, but nothing much came out of it:
Counting number of features within another feature using aggregating function in field calculator of QGIS?
Calculating point layer values within polygon features in QGIS 2


texts: Setting options in "PDF Export Options" to export an atlas to PDF using PyQGIS

texts:
While saving an atlas as PDF in QGIS, I can manually change the settings shown the image using "PDF Export Settings" window.

How can I set those options using PyQGIS?

Minimal code:

layout_manager = QgsProject.instance().layoutManager()
layout = layout_manager.layoutByName('LAYOUT_NAME')

layout_exporter = QgsLayoutExporter(layout)
layout_exporter.exportToPdf(path, QgsLayoutExporter.PdfExportSettings())


texts:
I've been working on making mountain biking and hiking maps using Illustrator and I am now trying to move to QGIS so there's real GIS data powering the maps, making updating/scaling them easier than it is Illustrator.

One issue I've had is that when trails converge the trail creators/maintainers usually give the trails two colors while they are joined. Here's an example from the Illustrator version:

I'm using the following solution to style the trails with two colors, but I'm curious if there are other ways to accomplish this:

The trail features have two color attributes: color and color2
I enter a string value for the "color" of the trail: red, blue, etc.
If the trail segment has a secondary color I put it in "color2"
I style the layer using the "color" attribute
I duplicate the trail layer and change the styling of the new layer to use "color2" instead of "color".
I change the duplicate layer to use a dashed style of line instead of solid.

Now both layers are pointing to the same Shapefile but styling them differently.

Is this the best way to accomplish this type of styling?


texts: Using Google Earth based dynamic web generated KMZ file (FCCInfo tower data) in QGIS 3.10.6

texts:
I'm trying to figure out how to use the downloadable kmz file in QGIS. The file is usable and searchable in Google Earth. When I try to import this to QGIS it says that it's an invalid format. The data is essentially all towers within the US registered with the FCC.
I assume this is because the content is fetched in a web database then generated on google earth. I assume if you can get a similar thing with google satellite imagery or terrain data, there should be a way to direct an XYZ tile to a similar url. I don't fully understand how the file works in Google Earth so I'm not sure exactly what to be searching for.
I've tried unzipping the kmz and importing the kml which doesn't work either. I've searched around for similar datasets, this FCC obstacle data is the closest, but does not have the full data, and it would also slog down the system so bad to have the full data. So, having it dynamically generated like a XYZ tile would be much better.


texts: setting label expressions causing duplicate labels

texts:
kind of new to VB.NET.

Trying to set my label expression for a layer. Essentially a user picks a layer from a combo box, hits a button, expression is applied on the layer and then the single shp file is converted to annotation. my issue is it seems to work for 75%. I am labeling creeks and want to label the [TYPE] vs the [NAME]. It seems to label the Type for most but for whatever reason there are a few names being labeled as well. Did I screw up some where?

it seems to label the DEFAULT field as well as my expression when I just want it to label my expression.

Shared Sub ConvertLabelsToAnnotationSingleLayerMapAnno(ByVal pMap As IMap, ByVal layerIndex As Integer)
    Dim pConvertLabelsToAnnotation As IConvertLabelsToAnnotation = New ConvertLabelsToAnnotationClass()
    Dim pTrackCancel As ITrackCancel = New CancelTrackerClass()
    'Change global level options for the conversion by sending in different parameters to the next line.
    pConvertLabelsToAnnotation.Initialize(pMap, esriAnnotationStorageType.esriMapAnnotation, _
                                          esriLabelWhichFeatures.esriVisibleFeatures, True, _
                                          pTrackCancel, Nothing)
    Dim pLayer As ILayer = pMap.Layer(layerIndex)
    Dim pGeoFeatureLayer As IGeoFeatureLayer = pLayer

    If Not pGeoFeatureLayer Is Nothing Then
        Dim pAnnotateLPC As IAnnotateLayerPropertiesCollection
        pAnnotateLPC = pGeoFeatureLayer.AnnotationProperties
        Dim pAnnotateLP As IAnnotateLayerProperties
        pAnnotateLP = New LabelEngineLayerProperties
        Dim pLabelELP1 As ILabelEngineLayerProperties2
        pLabelELP1 = pAnnotateLP
        pLabelELP1.Expression = String.Format("[Type]")
        pAnnotateLPC.Add(pAnnotateLP)
        pGeoFeatureLayer.DisplayAnnotation = True

        Dim pFeatureClass As IFeatureClass = pGeoFeatureLayer.FeatureClass
        'Add the layer information to the converter object. Specify the parameters of the output annotation feature class here as well.
        pConvertLabelsToAnnotation.AddFeatureLayer(pGeoFeatureLayer, pGeoFeatureLayer.Name + "_Anno", _
                                                   Nothing, Nothing, False, False, False, False, False, "")
        'Do the conversion.
        pConvertLabelsToAnnotation.ConvertLabels()
        'Turn off labeling for the layer converted.
        pGeoFeatureLayer.DisplayAnnotation = False
        'Refresh the map to update the display.
        Dim pActiveView As IActiveView = pMap
        pActiveView.Refresh()
    End If
End SubationSingleLayerMapAnno(pMap, cmbAOA.SelectedIndex)


texts: Linking Excel Sheet to Point Data reference table using ArcMap?

texts:
I am trying to join an excel spreadsheet to the attribute table of project locations in ArcMap. The excel CSV file is in my layers and I can open and view it as an attribute table.
However, it is when I try and join the two attribute tables by FID that I encounter some problems.
After selecting that I want to join by the FID field, I am unable to select the corresponding field in the table. (See pictures)

Essentially I a trying to create a reference map of all my projects. I dropped points at all the project locations on ArcMap and now I would like to link the data table with all the project info to those points.

If there is an easier/better way to do this, I am open to suggestions. I have a Microsoft Access database with the project info in it, but I don't really know how to access and link it in ArcMap.


texts: Creating polygons representing specific range of values within raster using QGIS?

texts:
Using QGIS, I started with LIDAR points (every 5 m) on a 90’000 ha area and then converted it in a raster of elevation (interpolating to fill some holes in LIDAR data). Now, I want to get a polygon representing the surface for which the cell value is equal or below a given threshold value.

I tried two methods:

Creating contours at the threshold value from the raster and converting these contours to polygons.
Polygonizing the raster reclassified as zeros (cell > x) and ones (cell =< x) and then simplifying it to a more contour-like shape.

Using method 1, the result is really close to what I want. Polygons are created but polygons that should be holes within a polygon, or islands within a hole, are created as separate polygons. Is there a way to convert contour lines into polygons while keeping holes and islands? If not, is there a way to correct this afterwards?

Method 2 seems a good approach in theory although very time consuming (in processing time). I found instructions proposed previously by « Lars » using GRASS (GRASS v.generalize method=douglas not working as expected). I’m satisfied with the result I got using a sampled polygon, but it doesn’t seem to work when using complete dataset : the processing of « v.generalize » seems to either crash or stop, and it happens when topology is building. I wonder if it’s only a matter of file size and computing limitations or if it’s a little more complicated than that.

v.generalize -c input=Vectorised_RASTER@username output= generalized_RASTER method=douglas threshold=5


texts:
I have a number of rasters that I need to calculate histograms for. These are continuous rasters from 0 to 294.208 - many decimal places.  I need to be able to classify the histogram down to at least the first decimal.

But! I cannot get the calculate statistics tool to run.  It says that the file cannot be dropped, yet doesn't say why.

I have tried:
-Moving the tif raster out of the geodatabase.
-Using Int() to make it an integer raster
-Converting to a GRID raster.

Nothing works! I have no idea what the tool wants, and the help file doesn't say. Can you help?

Running ArcGIS 10.2 for Desktop.


texts: Can ArcMap factor in legends and data frames when placing labels?

texts:
Is there a way to factor data frames and the legend when having placing labels using the Maplex Label Engine?

I am talking specifically about the labels for roads and contours. I really don't want to have to convert the labels to annotation and manually move them, and these are the labels that always seem to find a way to hide partially under an inset map, even when there is no apparent need to have the label in that one exact spot.

My current way around this is to create a polygon with no fill or border that fits under the inset and legend and give the polygon a suitable weighting. This works, and sometimes is quite handy at making the map look decent to show what is underneath the legend/inset map, but it seems like the placement of the labels should already factor in the map elements, and the process is becoming tedious enough to ask the question here.

Am I missing something in the labeling toolbar that allows this?


texts: Creating polygon GeoJSON from PostGIS and rendering in Leaflet

texts:
Drawing the polygons on QGIS and then import the layer directly on the PostgreSQL / PostGIS db. Then I create a JSON file via PHP and read it to the Leaflet. Everything works but zooming I see a lot of difference in rendering of the polygons which are much more simplified in the web-gis.
I set "smoothFactor": 0 but it's the same.

    function strutture_pol () {
    global $conn;
    $sql = 'SELECT strutture_pol.*, public.ST_AsGeoJSON(public.ST_Transform((geom),4326),6) AS geojson, datazione.id_data_gen as data_gen, datazione.id_data_spe as data_spec, datazione.id_pref1 as prefisso1, datazione.id_sec_in as datain,  datazione.id_pref2 as prefisso2, datazione.id_sec_fin as datafin   FROM strutture_pol ';
    $sql .= 'INNER JOIN datazione ON strutture_pol.id_data = datazione.id ';
    $sql .= 'WHERE strutture_pol.id_data = datazione.id ';
    $stmt = $conn->prepare($sql);
    $stmt->execute();
 # Build GeoJSON feature collection array
    $geojson = array(
   'type'      => 'FeatureCollection',
   'features'  => array()
    );
# Loop through rows to build feature arrays
    while ($row = $stmt->fetch(PDO::FETCH_ASSOC)) {
       $properties = $row;
        # Remove geojson and geometry fields from properties
        unset($properties['geojson']);
        unset($properties['geom']);
          $feature = array(
            'type' => 'Feature',
             'geometry' => json_decode($row['geojson'], true),
             'properties' => $properties
             );
          # Add feature arrays to feature collection array
           array_push($geojson['features'], $feature);
      }
      header('Content-type: application/json');
       header('Access-Control-Allow-Origin: *');
      echo json_encode($geojson, JSON_NUMERIC_CHECK);
    }


texts:
This python script for qgis should reproject all the raster files in a folder and resample it to 10 m. Raster file type should be Float32. The script doesn't produce any result.

import os
import processing
from qgis.core import *

inputDir = "R:\\LAI_Biom_Maps\\kriging_Gladys\\in"

outputDir = "R:\\LAI_Biom_Maps\\kriging_Gladys\\out"

for i in os.listdir(inputDir):

    layercount = os.path.join(inputDir, i)
    raster_layer = QgsRasterLayer(layercount, "raster")

    fileOut = outputDir + "\\" + i

    processing.runalg("gdalogr:warpreproject", \
    raster_layer, \
    "EPSG:4326", \
    "EPSG:27700", \
    10, \
    0, \
    None, \
    5, \
    fileOut)


texts:
I am using GRASS 7 for windows, and I'm having trouble running v.in.ascii, although I'm fairly sure I've run it successfully before. When I run v.in.ascii, I get an error: Windows can't open this file.

This happens through the command line and the gui. Other v.*.exe commands seem to run properly.

I thought maybe there was a path clash somewhere, so I generally cleaned up my path variable (e.g., removing OSGeo4W/bin). I also reinstalled GRASS 7. I searched GRASS bugs, but didn't find something similar. Other than that, I haven't been able to find anything on google to help explain this, and I've run out of ideas.

Recently, IT installed hard drive encryption software. This and a couple problems popped up about the same time, so I have a sneaking suspicion that it's related.

Edit: some additional information that @GeoStoneMarten helped me tease out:

From the GRASS shell, if I cd to C:\grass70\bin, then v.in.ascii works
v.in.ascii is the only v.* command that gives me this error
C:\grass70\bin is being added to the path variable correctly upon calling grass70.bat


texts:
I'm trying to open Sentinel-2 bands directly from the .zip L2A files with GDAL, but this only returns a None. Unzipping the folder prior to reading the file yield better results, but still randomly returns some None from time to time. I'm using gdal 2.2.2.

import numpy as np
import osgeo.gdal as gdal
import zipfile, os

from libs.generic import filemanager as fm
os.environ['CPL_ZIP_ENCODING'] = 'UTF-8'

class s2img:

def __init__(self):
    return None

def readL2A(self, zippath):

    try:
        zipf = zipfile.ZipFile(zippath, 'r')
    except:
        raise IOError("Unable to open ZIP-file.")

    flist = zipf.namelist()
    zipf.close()

    bandnames = []
    bandnames.append( ['B02','B03','B04','B08'] )
    bandnames.append( ['B05','B06','B07','B8A','B11','B12','SCL'] )
    bandnames.append( ['B01','B09'] )

    bandres = ['10m','20m','60m']

    fnames = [f for f in flist if 'IMG_DATA' in f]
    fnames = [f for f in fnames if f.endswith('.jp2')]
    for r in range(len(bandnames)-1):
        res = bandres[r]
        resbnames = bandnames[r]
        fresnames = [f for f in fnames if res in f]

        for b in range(len(resbnames)):
            fname= [f for f in fresnames if resbnames[b] in f ][0]
            chn_fn = "/vsizip/%s/%s" % (zippath,fname)
            gobj = gdal.Open(chn_fn, gdal.GA_ReadOnly )

            # DO STUFF
    return None


texts:
I have a bunch of raster maps from a general plan of a city.

Examples:

Areas reserved for federal, regional and municipal buildings
Special land use zones
Zoning
Risks of occurrence of emergencies
Power supply scheme
Industry, energy sector and tranportation zones
Transportation infrastructure scheme
Water supply and sewage scheme

I need to vectorize them such that I can store all elements (lines, polygons, points) of the map listed in the legend in a database like PostGIS such that I can

find all elements of a master plan, which affect a certain region on the map (e. g. if I have a house at location X and the government wants to build a road in the same place, I want to notice it) and
compare two versions of a master plan by running database queries.

The ultimate purpose is to bring the data contained in these maps into a format that can be understood by non-technical people - inhabitants of the city so that they can express agreement or disagreement with the general plan.

One obvious approach is to simply lay raster map over a vector map of the region (there is a vector map of that city in OSM) and manually transfer the lines from raster to vector form.

What tools (preferably open source) can I use to vectorize the maps manually and save the data in PostGIS?
Are there any alternatives, which are faster? Can this process be automated (fully or partially) using open source or inexpensive commercial software?

Update 18.01.2013: FYI: I asked the copyright holder (city administration) of the maps to provide them in vector format. I'll wait for their answer and then decide on further actions.

Update 27.01.2013:

Today I made a first experiment in georeferencing for these maps.

Inputs:

OpenStreetMap map
Raster map

Process: You can see how I do georeferencing in this video.

Results:

GCP points for both maps
Modified raster map

I have following questions:

How many points do I need to enter in order to have good enough quality of the mapping? Are there any rules of thumb for that?
What can I do to improve the quality of the mapping (remove mapping defects shown below) ?

Sample mapping defect #1

Sample mapping defect #2


texts: Error when trying to publish on QGIS Cloud

texts:
I get this error when trying to publish on QGIS Cloud:

Name contains invalid characters
Traceback (most recent call last):
  File "/Users/rubenkalmbach/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/qgiscloud/qgiscloudapi/qgiscloudapi.py", line 657, in request
    response = urllib.request.urlopen(req).read()
  File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/urllib/request.py", line 223, in urlopen
    return opener.open(url, data, timeout)
  File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/urllib/request.py", line 532, in open
    response = meth(req, response)
  File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/urllib/request.py", line 642, in http_response
    'http', request, response, code, msg, hdrs)
  File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/urllib/request.py", line 570, in error
    return self._call_chain(*args)
  File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/urllib/request.py", line 504, in _call_chain
    result = func(*args)
  File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/urllib/request.py", line 650, in http_error_default
    raise HTTPError(req.full_url, code, msg, hdrs, fp)
urllib.error.HTTPError: HTTP Error 422: Unprocessable Entity

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/rubenkalmbach/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/qgiscloud/qgiscloudplugindialog.py", line 659, in publish_map
    map = self.api.create_map( self.map(), fname, config)['map']
  File "/Users/rubenkalmbach/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/qgiscloud/qgiscloudapi/qgiscloudapi.py", line 304, in create_map
    content = request.post(resource, data)
  File "/Users/rubenkalmbach/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/qgiscloud/qgiscloudapi/qgiscloudapi.py", line 569, in post
    return self.request(resource, method='POST', data=data)
  File "/Users/rubenkalmbach/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/qgiscloud/qgiscloudapi/qgiscloudapi.py", line 683, in request
    raise BadRequestError(msg)
qgiscloud.qgiscloudapi.qgiscloudapi.BadRequestError: Name contains invalid characters


texts: Specify a field in attribute table with raster calculator

texts:
How can I specify the attribute field to use in the raster calculator of ArcMap 10.2?
In 9.x, I used the dot notation.

i.e. To use the "New_Value" field, I could do the following:

[MyRaster].New_Value + 1 = new raster

This doesn't appear to work in 10.2?


texts: Matching extent (samples & lines) to get multi-temporal stack

texts:
I downloaded several Landsat scenes (4x3) to create some multi-temporal stacks to see the phenological difference in April/March - mid summer - autumn. I actually use The EnMap Box which has been developed by my university.

I alreadly stacked the 7 bands of each scene without problems. Now the programm tells me that the 3 scenes don't have the same spatial size. In the .hdr data the samples and lines are not exactly the same what should cause the problem.

How can I solve that issue?


texts: How to translate a QGIS project using a *.qm file? / my QGIS project is not translated?

texts: Usually when you create a QGIS project, you will name the attributes, layer etc. in your own language. You might want to deliver the project to another user in another language.


texts:
I am trying to write a processing script that takes a template file as input and writes out the current canvas content as a PDF. The processing script works fine. But when I use the script in a model, the map output is blank.

The model starts with a blank canvas. A layer is generated by the model and loaded in the canvas. Following that step the algorithm containing the script is called to generate the pdf. I have also set the previous step of the model as the parent, but no luck.

Below is the processing script.

from PyQt5.QtCore import QCoreApplication
from qgis.core import QgsProject, QgsLayout, QgsLayoutExporter, QgsReadWriteContext, QgsMapSettings, QgsProcessingAlgorithm, QgsProcessingParameterFileDestination, QgsProcessingParameterFile
from PyQt5.QtXml import QDomDocument

class ExportLayoutAlgorithm(QgsProcessingAlgorithm):
    """Saves the attributes of a vector layer to a CSV file."""
    OUTPUT = 'OUTPUT'
    INPUT = 'INPUT'

    def flags(self):
          return super().flags() | QgsProcessingAlgorithm.FlagNoThreading

    def initAlgorithm(self, config=None):
        self.addParameter(
            QgsProcessingParameterFile(
                self.INPUT,
                self.tr('Template File'),
            )
        )

        # We add a file output of type CSV.
        self.addParameter(
            QgsProcessingParameterFileDestination(
                self.OUTPUT,
                self.tr('Output File'),
                'PDF File (*.pdf)',
            )
        )

    def processAlgorithm(self, parameters, context, feedback):
        source = self.parameterAsFile(parameters, 'INPUT', context)
        pdf = self.parameterAsFileOutput(parameters, 'OUTPUT', context)
        # Load template from file
        p = context.project()
        l = QgsLayout(p)
        l.initializeDefaults()

        with open(source) as f:
            template_content = f.read()
        doc = QDomDocument()
        doc.setContent(template_content)

        # adding to existing items
        items, ok = l.loadFromTemplate(doc, QgsReadWriteContext(), False)
        title = l.itemById('Title')
        title.setText('My Map')
        map = l.itemById('Map')

        exporter = QgsLayoutExporter(l)
        exporter.exportToPdf(
             pdf, QgsLayoutExporter.PdfExportSettings())

        return {'OUTPUT': pdf}

    def name(self):
        return 'export_layout'

    def displayName(self):
        return self.tr('Export Layout As PDF')

    def group(self):
        return self.tr(self.groupId())

    def groupId(self):
        return ''

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

    def createInstance(self):
        return ExportLayoutAlgorithm()


texts: How to get the nodes on the absolute boundary of a TIN with ArcObjects

texts:
I am trying to get only the nodes on the absolute boundary of a Tin. But GetNeighbor or LeftTriangle methods return edges whose nodes have NaN Z and very weird X-Y coordinates, so the row is inserted in the feature class without a feature.

Is there any other method or am I doing anything wrong? Thanks, Mar

Dim i As Integer
    For i = 1 To pTinAdvC1.TriangleCount

Set pTinTriangleC1 = pTinAdvC1.GetTriangle(i)

    Dim j As Integer
    For j = 0 To 2

        If pTinTriangleC1.Edge(j).GetNeighbor Is Nothing Then

        Dim pTinEdge As ITinEdge
        Set pTinEdge = pTinTriangleC1.Edge(j)

        Dim pBorderP As IPoint
        Set pBorderP = New Point
        Set pZAwarePoint = pBorderP
        pZAwarePoint.ZAware = True

        pTinEdge.FromNode.QueryAsPoint pBorderP

        Set pPointLayer = pScene.Layer(2)
        Set pPointFC = pPointLayer.FeatureClass
        Set pFeat = pPointFC.CreateFeature
        On Error Resume Next
        Set pFeat.Shape = pBorderP
        pFeat.Store

        End If

    Next

Next


texts:
I can observe a strange behavior with my PostGIS query. I am trying to accomplish the following:

I have one polygon layer
I have 2 different line layers

layer with horizontal lines
layer with 45° lines

the goal is to cut the lines at the polygons to get the line parts of each line per polygon

Here is an example what that looks like:

I am using a query which looks like this:

CREATE TABLE split AS
    ( SELECT row_number() over() AS id,
             a.id AS old_id,
             st_intersection(a.geom, b.geom) AS geom
     FROM lines a,
          polygons b
     WHERE st_intersects(a.geom, b.geom) ) ;

The query takes ~3 minutes for horizontal and vertical lines but 5+ hours (cancelled there) for 45° lines.

Can anyone think of a possible cause? All layers have spatial indices, pk and so on. I am using a Postgres 9.5 with PostGIS 2.2.

I am really clueless about that behavior.


texts: Distance between centroid and farthest point of polygon

texts:
I have a village polygon layer which has more than 6,00,000 records.
I have calculated centroid of each village.
I want to find the distance between centroid and the farthest node of each polygon.
Check the image below for reference. Black lines are polygon boundaries.


texts: Load Maperitive's tiles.json into HTML for getting bounding parameters

texts:
I've created a local tileset of a city in Maperitive, and a basic leaflet.js webpage to browse it offline. There's a tiles.json also created by Maperitives that has all the essential bounding and zoom limits that I want.. at present am having to copy-paste them manually into the html file. Please show how to import the tiles.json and load the limits etc directly from there. This way a generic html is available that will display whatever content the tiles folder with it is carrying. (note: this isn't a typical GeoJSON.. doesn't have placemarks or polygons.)
HTML code:

<link rel="stylesheet" href="lib/leaflet.css" />

tiles.json contents: (generated by Maperitive when generating tiles)

// Generated by Maperitive v2.3.32
// For more information about TileJSON format, visit https://github.com/mapbox/tilejson#readme
// TODO: Update the 'tiles' to reflect your actual path on the Web server
// TODO: Update the default map location (longitude, latitude, zoom)
{
    "tilejson":"1.0.0",
    "name":"My Map",
    "description":"Made with Maperitive",
    "attribution":"Map data © OpenStreetMap contributors",
    "tiles":
    [
        "tiles/{z}/{x}/{y}.png"
    ],
    "minzoom":8,
    "maxzoom":15,
    "bounds":
    [
        73.741131698096311,
        18.443790655688225,
        73.955365096533811,
        18.672919555595893
    ],
    "center":
    [
        73.848248397315061,
        18.558355105642057,
        8
    ]
}

Requirement : HTML should load the tiles.json file, and get the center, zoom, minZoom, maxZoom, maxBounds data from there rather than author having to manually type them in.

Edit: Since the tileJSON format has been created by Mapbox, I put this same use-case in another question, targeted at Mapbox, and some solutions seem to be coming together there. So please check out : Load external tileJSON file in Mapbox


texts: Is it possible to merge raster base maps with a DTM?

texts:
In ArcMap 10.1

I have two raster layer 1) a raster base map and 2) a DEM raster.

Is it possible to merge the two to create an attractive visual showing the background map and terrain?  i.e. in Global Mapper there is an option to multiply layers together.

I have tried putting the base map on top making it slightly transparent so the underlying DTM is visible underneath.  But the result is quite faint thereby missing off base map detail.


texts: How to use GeoPandas buffer function to get buffer zones in kilometers?

texts:
I am having trouble with using GeoDataFrame.buffer function. That is probably because of my lack of understanding crs. Here is the case:
I have events dataset (points longs & lats) and locations dataset (points longs & lats). I want to check which of the events occurred within 2km radius circle from locations. I know how to perform spatial join in geopandas, but the problem is in buffer(distance, resolution) function. When creating both GeoDataFrames I used crs: 'epsg:4326' and in distance parameter use 0.2 (without any good reason). So my question is:
How to adjust crs and distance parameter to make it in kilometers?
Here is my sample code. The intersection should contain only events with id 1 & 2. However it does contain event 3 which is about 10km far from any of locations. Also it does not contain event with id 4, which is good, because that one is hundreds kilometers from any of locations.
events = pd.DataFrame([
    {'id': 1,
     'longitude': 54.606402,
     'latitude': 18.347709},
    {'id': 2,
     'longitude': 54.604681,
     'latitude': 18.346534},
    {'id': 3,
     'longitude': 54.544322,
     'latitude': 18.293896},
    {'id': 4,
     'longitude': 50.843233,
     'latitude': 11.130688}
])
locations = pd.DataFrame([
    {'id': 'a',
     'longitude': 54.604972,
     'latitude': 18.346815},
    {'id': 'b',
     'longitude': 54.605917,
     'latitude': 18.347249}
])

locations_gpd = gpd.GeoDataFrame(locations,
                                 geometry=gpd.points_from_xy(locations.longitude, locations.latitude).buffer(0.2, 16),
                                 crs={'init': 'epsg:4326'})

events_gpd = gpd.GeoDataFrame(events,
                              geometry=gpd.points_from_xy(events.longitude, events.latitude),
                              crs={'init': 'epsg:4326'})

intersection = gpd.sjoin(locations_gpd, events_gpd, how='left')


texts: Export csv-coordinates to JSON - square brackets

texts:
I´ve got multiple point-layers that I want to join in a single JSON-file used for webmapping. The layers is to be listed in column by column like this:

[55.778659,12.151137  ],    [55.77272358,12.14963029  ],

Is this the right setup when it comes to separators (commas)?

It's for Leaflet.


texts: Is there a way to draw a straight line on a map which also represents a straight line on the earth surface?

texts:
I'm trying to understand more on projections.

In this example, why must one calculate the curvature of the earth (great circle) if the map is already projected? Is there a way to draw a visibly straight line on that map which also represents a straight line on the surface of the earth? Why or Why Not?


texts: How to get the URL path to a layer in Geoserver?

texts:
From the Geoserver Layers Preview page, is there a way to get the URL of the services in the 'All Formats' drop-down?

Alternately, is there a way to change which formats are listed in the 'Common Formats' column?

I've been sharing the Layer Preview page to external users as a way to give a quick way to share our data in various formats, but would like to, for example, give the URL for a CSV or SHP rather than send instructions to pull the drop-down menu.


texts:
How can I do layerstack in QGIS 2.16.2?

I do not find the option under the Raster>Miscellaneous>Merge options.


texts:
I've been given Omega/Phi/Kappa data from drone surveillance photos and I need to convert it to Euler angles (or a quaternion, I'm not fussy).

I've been given the opk data in X right Y up format, and I need to convert this to an Euler angle in a X right Z up coordinate system. This is taken from the Autodesk Forge Viewer which is matching the Revit coordinate system.

I'm trying to match the camera angle in the Forge viewer to the drone's camera angle.

Am I right in assuming the opk angle is basically an Euler angle already?

In which case, should I use the default XYZ Euler order?

Finally, how do I convert X right Y up to X right Z up?

I assumed just switching the Y/Z axes was the correct transform here but that's not giving me the right result.

Edit: here's an image of the coordinate system I need to match:


texts: Access control for SaaS-based GIS, which is the right way to go?

texts:
I'm currently researching on some access control models for SaaS-based GIS, and in the future would like to implement one. The common use cases would be like: some of the services, resources, layers, features, areas on map, etc. will be accessible only by authorized users.

I have searched for a while and find the OGC GeoXACML implementation specification meet our needs quite well and looks promising. And also there's quite a few paper discussing about Geoxacml that I can refer to. But then I find another article which propose to use oracle VPD (we are using oracle spatial), and it points out that the GeoXACML has some inefficient consequences. There is another tool that might solve the problem is the GeoShield (https://sites.google.com/site/geoshieldproject/), but it seems to be in the early work and don't have to much documentation. So which option would you suggest? Is there other options available? If geoxacml is the right way to go, how long do you estimate it will take to implement it? Any response would be appreciated.


texts: Understanding Label syntax in ArcGIS for Desktop?

texts:
In the context of this question

Label specific lines in attribute table with syntax- ArcGIS?

i build this python phrase in ArcGIS 10.3:

def FindLabel ([OBJECTID]):
    if long([OBJECTID])==125 and long([OBJECTID]) ==190:
        label = [SHEM_KID]
    else:
        label = None
    return  label

in order to label only OBJECTID number 125 and 190.
My problem is that nothing happen, and all label disappear from the map. I try also to remove "long" from the code and it also  didn't worked.

Is it possible to do it with python code?

I tried this code also but with no success:

def FindLabel ([OBJECTID], [SHEM_KID]):
    if ([OBJECTID])=='125' AND ([OBJECTID]) =='190':
        label = [SHEM_KID]
    return  label


texts:
I'm using the maptools package in R to merge some data to a shapefile.

When I open the original input shapefile in QGIS, QGIS automatically picks up on the CRS.

When I open the output shapefile from R, I need to select the CRS.

This happens even if I specify the CRS in my R code. Is there a solution to this?

# Load libraries ----------
library("maptools")

# Read shapefile ----------
LSOA <- readShapeSpatial("c:\\Temp\\LSOA11.shp")

# Set CRS ----------
# http://spatialreference.org/ref/epsg/27700/
OSGB36.CRS <- CRS("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs")
proj4string(LSOA)<-OSGB36.CRS

# Check CRS ----------
proj4string(LSOA)

# Filter ----------
LSOA$Flag <- grepl("YORK",LSOA$LSOA11NM, ignore.case = TRUE)
LSOA <- LSOA[LSOA$Flag==TRUE,]

# Plot to check filtering ----------
plot(LSOA, col = 'light blue')

# Export shapefile ----------
writeSpatialShape(LSOA,fn="c:\\Temp\\LSOA_Fuel_Poverty")


texts:
I'm using the Python GEE API.
I'm able to crawl in my tasks through this :
 #search for the task in task_list
    filename = construct_filename(country_code, year)
    tasks_list = ee.batch.Task.list()
    current_task = ''
    for task in tasks_list:
        if task.config['description'] == filename:
            current_task = task
            break

Is it possible to do the same in my Assets and if yes how ?


texts: LiDAR data density, returns per square meter

texts:
How can I compute a raster showing the number of return points per square meter using SAGA GIS or other open source software?

My goal is to obtain a mask which would show me the places (cells) where there were no hits and would help me in interpreting LiDAR-derived data.


texts: Masking all object less than given connected pixel count

texts:
I am working on the nature of deforestation in the Philippines using the Global Forest Change (GFC) dataset of Hansen et al. (2013). I am trying to depict all connected pixels that has a minimum forest clearing of 5 hectares, which is approximately 56 pixels. First problem is that when I change the resolution of the map dashboard of the GEE, it seems that patches measuring below the 56 pixels slowly shows as I enlarge the image. One solution I found over the internet is to convert the projection into a nominal projection. However, this does not work well as I encountered a problem saying that the number of pixels requested from that variable exceeded the maximum allowed (2^31). I also attach the code that I made so that it would be convenient to comment down below.
var gfc2019: Image "Hansen Global Forest Change v 1.7 (2000-2019)"
var cc = ee.Number (30);
var pixels = ee.Number (6);
var lossPixels = ee.Number (56);
var forest = gfc2019.select (['treecover2000']);
var minCanopy = forest.gte(cc).selfMask ();
var contArea = minCanopy.connectedPixelCount ();
var minForest  = contArea.gte (pixels).selfMask();

var treeLoss = gfc2019.select (['lossyear']);
var treeLoss01 = treeLoss.eq(14).selfMask ();
var minLoss01 = minForest.and(treeLoss01).rename ('Forest loss in 01').selfMask ();
var contLoss = minLoss01.connectedPixelCount();
var largeLoss01 = contLoss.gte(lossPixels).selfMask();

var prj = gfc2019.projection();
var scale = prj.nominalScale();

Map.addLayer (minForest,{palette: ['green']}, 'Forest');
Map.addLayer (largeLoss01.reproject(prj.atScale(scale)),
{palette: ['red']}, 'Large Scale Loss');


texts: Performance problem to access QgsLineString vertex

texts:
I am writing a QGIS processing plugin.  I need to have accessible to Python in a list of xy tuple [(x0,y0),(x1,y1),(x2,y2),...], the x,y coordinate of each vertex of each QgsLineString.  The fastest way I found yet was with the following code:
qgs_pnts = qgs_line_string.points()
xy = [(qgs_pnt.x(), qgs_pnt.y()) for qgs_pnt in qgs_pnts]

Where qgs_line_string is a QgsLineString object.  I have to work with a lot of information. With the Python profiler I was able to estimate that the code above was taking 25% of the execution time just to load the QgsLineString into a list of xy tuple.  I understand why it's taking so much time the comprehension list is making two C calls for each vertex of the QgsLineString object.
My question: is there another way in pyqgs to load faster the vertex of a QgsLineString object  into a Python list of xy tuple?


texts: Splitting field with comma as separator into two fields in QGIS

texts:
How to split a "Field" into "Field 1" and "Field 2" when the comma is a separator in the field?

I am using Field Calculator in QGIS.


texts: Using M 'Measure' value in SQL Server 2008

texts:
MS Sql Server 2008 geometry type supports adding Z and M values to geometry objects, including linestrings. However, I cannot find any way to make use of these values, since methods such as      .STIntersection() return a 2D point.

My end goal is to find where a given point intersects with a linestring and then find the equivalent Z and M values.

Googling and scouring Technet has yielded no useful information other than to make me wish I was using PostgreSql instead, which seems to include the functionality I require (unfortunately I am constrained to Sql Server).

Does anyone know how to actually make use of the Z and M values within linestrings, or of any toolsets that provide this functionality?


texts: Minimize time for snow-covered image with `i.atcorr` module

texts:
Did anyone notice the time taken by i.atcorr module when running for a snow-covered image?

It is noticed that two to three hours for 30m resolution with the image of 3656*1000 size.

Both for VNIR(60bands) and SWIR(256 bands) it is taking a long time. What could be the reason?


texts: Intersect outputs shapefile with null values

texts: I have a model that iterates through a folder of field boundaries (polygons) and produces an intersect against a shapefile of points. The output is then a point shapefile for each boundary. Note: There are not points in every field boundary. My goal is to get an output of points ONLY if the points ALSO are in a field boundary. It works, except I get a set of points that have empty values for field boundaries that do not have point in them. How can I eliminate these empty point features?


texts: How to loop through the raster layers and create a PNG for each?

texts:
In QGIS 2.8.1 I have a number of raster layers covering the same extent, but for each year 1975-2014 (times three categories).

I am creating one print composer for the extent, but I need to print multiple PNG's - one for each raster layer.

How do I loop through the raster layers and send an output PNG for each?
Should I produce a python script (help needed) or somehow use the Atlas generator?


texts: Arcpy: Exporting PNGs with high quality using data driven pages

texts:
I am using data driven pages and want to automatically export my layouts/extents to PNG. Since you have to use arcpy to export to PNG with data driven pages, I used following code, as explained here: http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00sm00000008000000

mxd = arcpy.mapping.MapDocument("CURRENT")
for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):
    mxd.dataDrivenPages.currentPageID = pageNum
    arcpy.mapping.ExportToPNG(mxd, r"C:\Temp\ParcelAtlas_Page" +str(pageNum) + ".png")
del mxd

The problem is that the quality is very bad when I export like this. If I export to PNG manually, and set for example 300 dpi, it is way better than exporting automatically.

Would you know a solution to export in better quality using arcpy and data driven pages in ArcMap?


texts: Checking shape orientation using ArcGIS Desktop?

texts:
I have a shapefile with many minimum bounding box (envelope) polygons and I need to find out, if the boxes are rather landscape or portrait format. (See pic below, one of the selected boxes should be marked as landscape, the other as portrait)

I got the width and length of the boxes, however, the shorter line is always taken as the width by ESRI default so I cannot find out if the box is landscape or portrait using these values.

Anyone has an idea on how to solve this in an easy way?

All I can come up with is computing the corner coordinates, but not sure how to do it exactly.

I work with ArcGIS Desktop 10.1.


texts: Copy attribute values between same ID featutes QGIS

texts:
I'm trying to copy multiple attribute values between two features in a geojson with the same value of another field (like an id).

For example:
- id= 1 field = A

if in another record I repeat id =1 automatically set field = A

I think with autofields plugin i could make something like that but i don´t know how to build the expression i have to use.


texts: Quadtree based polygon indexing? (Python/Java)

texts:
I was wondering how to best represent the inside of any given polygon (GeoJSON coordinates) with some sort of indexing, like a grid or quadtree? See image below for an approximation of what I have in mind.

Is there any library in Python or Java that can help me with this?


texts: Calculating area by classified image in Earth Engine

texts:
I want to calculate the area for each classified color in earth engine.
Here are my codes.

var image = ee.Image(ee.ImageCollection('LANDSAT/LC8_L1T_TOA')
.filterBounds(roi)
.filterDate('2015-01-01', '2015-12-31')
.sort('CLOUD_COVER')
.first());
print(image);

Map.addLayer(image, {bands: ['B4', 'B3', 'B2'], max: 0.3}, 'image_L8');

var clip = image.clip(table);  //cliping image to shapefile

var newclassified =
vegetation.merge(water).merge(agri_bare).merge(builtarea);
print(newclassified);
var bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'B11'];

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

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

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

Map.centerObject(newclassified, 11);
Map.addLayer(clip,
{bands: ['B4', 'B3', 'B2'], max: 0.4},
'Landsat image');

Map.addLayer(classified,
{min: 1, max: 4, palette: ['green', 'blue', 'yellow', 'red']},
'classification');

Now I want to calculate the area covered by green, blue, yellow and red palette.
What code I required to do that??


texts: How can I render latitude-longitude coordinates on a map with D3?

texts:
I am trying to overlay dots of arbitrary sets of longitude-latitude coordinates on top of a United States map.

So far I have found the D3 cartography example, but when I try to place dots on X,Y pixel coordinates they appear way off the canvas.  I did find workshot lecture notes which include coordinate systems in d3, but I'm still unsure of how get lat/long coords to show up on the map.

This is what I have so far (pretty much only a map of the USA)

Advice on how to make this work would be greatly appreciated!


texts: How to use preconditions with row counts when some row counts = 0?

texts:
I have a model that iterates through a series of neighborhood polygons and clips several different datasets using the iterated neighborhood polygon. I am trying to determine the number of features that lie within the polygon and record these values to a table. So, my model clips the input layers (bike lanes, for example) to the shape of the neighborhood, counts the number of features using "get count", then uses "calculate field" to write these values to the table.

I obviously need the "get count" tool to execute before the calculate field tool, because I am using the row count as an inline variable for the field calculation. However, when I run the model, it does the field calculation BEFORE the get count tool runs. This makes my data useless.

I tried using a precondition to set the order of the model, but if there are no bike lanes in the neighborhood, the precondition kills the rest of the calculation because of the zero value instead of just writing the 0 value to the table and moving on (The process did not execute because the precondition is false). I cannot figure out how to control the order of the model without using a precondition. Can somebody tell me what I'm missing here?

Here is the section of my model that is breaking. In the actual model there are a lot of different clips that go on, and the model goes through all of the clips in the model, THEN to my Calculate Field, THEN to my Get Count.


texts: First call to transform() fails with inf, all subsequent calls are OK - what could be the reason?

texts:
I'm trying to get to the bottom of the following (weird) problem. When applying pyproj.Transformer transformation to a certain polygon, it fails the first time with (inf,inf) throughout. Second and all subsequent calls to the same work just fine!
I saw similar issues with inf which had to do with inappropriate CRS (out of bounds), but this is not the case here.
The question was originally prompted by geopandas' to_crs() failing occasionally; similar question by someone else here.
Example
import shapely, shapely.geometry, shapely.speedups, shapely.ops, pyproj
print('Shapely version:',shapely.__version__)
print('Shapely speedups enabled:',shapely.speedups.enabled)
#shapely.speedups.disable() # makes no difference
print('Pyproj version:',pyproj.__version__)
print('PROJ version:',pyproj.proj_version_str)

xy = [(532303, 181814), (532190, 181854), (532192, 181861), (532199, 181863),
      (532195, 181874), (532183, 181870), (532163, 181864), (532162, 181868),
      (532309, 181915), (532294, 181828), (532306, 181820)] #, (532303, 181814)]
geom = shapely.geometry.Polygon(xy)
OSGB1936 = pyproj.CRS('EPSG:27700')
WGS84 = pyproj.CRS('EPSG:4326')

proj = pyproj.Transformer.from_crs(OSGB1936, WGS84, always_xy=True)
func = proj.transform

# First call fails with inf, all subsequent calls are OK:
for ii in range(3):
    tp = shapely.ops.transform(func, geom)
    print(tp.exterior.coords[0])

Output:

Shapely version: 1.7.1
Shapely speedups enabled: True
Pyproj version: 2.6.1.post1
PROJ version: 7.1.0
(inf, inf)
(-0.0929010345368997, 51.519248856167074)
(-0.0929010345368997, 51.519248856167074)


texts: How to capture the event whenever the items inside a listview are checked

texts:
I'm developing a plugin in QGIS 2.8.2 using Pyqt 4.
Well the idea is using the processing 'qgis:extractbylocation' from the QGIS algorithm provider to capture the studied area

# Récupérer tous les iris qui s'intersectent avec la zone d'étude
res = processing.runalg('qgis:extractbylocation', irisLayer, layerZone, u'intersects', 0, None)

# Créer une couche vector issue de l'extraction des IRIS dans la zone
ObjectifVilleExtract.irisLayer = QgsVectorLayer(res['OUTPUT'], "IRIS dans la zone d'etude", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(ObjectifVilleExtract.irisLayer)

# Définir un ecodage Windows compatible avec le français
ObjectifVilleExtract.irisLayer.setProviderEncoding(u'iso-8859-1')

# Modifier l'encoage
ObjectifVilleExtract.irisLayer.dataProvider().setEncoding(u'iso-8859-1')

# Vérifier que le couche est valide sinon afficher que la couche est erronée
if not ObjectifVilleExtract.irisLayer.isValid():
    print "Layer failed to load!"

Then showing the results on a listview:

ObjectifVilleExtract.features = ObjectifVilleExtract.irisLayer.getFeatures()
ObjectifVilleExtract.model = QStandardItemModel(ObjectifVilleExtract.features)
for feature in ObjectifVilleExtract.features:
    item = QStandardItem(feature['nom_iris'])
    item.setCheckable(True)
    ObjectifVilleExtract.model.appendRow(item)
    self.dlg.listView_3.setModel(ObjectifVilleExtract.model)
    self.dlg.listView_3.show()

Now I want to capture checked item on the listview, to apply on its the following processes :

iter = ObjectifVilleExtract.features.getFeatures()
for feature in iter :
    q1 = "SELECT nom_com FROM france_full_data.contours_iris WHERE code_iris='" + feature['code_iris'] + "' "

I tried to use the following code :

for feature in ObjectifVilleExtract.features:
    while ObjectifVilleExtract.model.item(feature['nom_iris']) :
        if ObjectifVilleExtract.model.item(feature['nom_iris']).ischecked():
            ObjectifVilleExtract.iter = ObjectifVilleExtract.irisLayer.getFeatures()
            irisTab=[]
            for feature in ObjectifVilleExtract.iter :
                q1 = "SELECT nom_com FROM france_full_data.contours_iris WHERE code_iris='" + feature['code_iris'] + "' "

the result is empty !


texts:
I need to normalize the Haversine distance.

What is the maximum Haversine distance? The Earth circumference I guess? How much is it?


texts: Show only labels for a place if they wouldn't hide the label for a bigger place?

texts:
I have a layer of places with their name and population data. I'd like to dynamically label the places in QGIS: draw a label if it wouldn't cover up or hide the label of another feature with a bigger population.

So in the countryside, I might want a tiny village to show up, but not a village of the same size right next to the city. If I zoom all the way out, I'd expect the smallest villages to disappear, then all villages, then all the way out I'd only have the largest city in my area. I tried to play a bit with using the population value for priority, but wasn't able to make it work.


texts: Test Statistic, p-value interpretation for Local Moran's I

texts:
I am using Local Moran's I in the spdep in R.

See this question: https://stackoverflow.com/questions/44104008/spatial-autocorrelation-between-time-series-in-r?noredirect=1#comment75251231_44104008

Does this look right?
And I'm reiterating the question in the end of that post.

How do I collapse all the Iis to come up with one Ii for the entire time span for that cluster?

Averaging?

Fisher's method?

What do you do with a test statistic value that has some relevance (like Ii = -0.357682) but a high p-value (like p-value = 0.59281) {arbitrary values, but same general idea}?

Assume a 95% confidence interval.


texts: Microsoft Visual C ++ Crash when Install qgis.2.8.1

texts: I have installed and reinstalled QGIS ver.2.8.1. and always display a message as shown below. I hope someone can help me fix these problems.


texts: Is it possible to install zoo-project in Ubuntu 18.04?

texts: I want to install zoo-project in Ubuntu 18.04 (bionic), but it seems there aren't any programs to do that, for example I have a problem with zoo-kernel install because there is not "ibmozjs185-dev" for bionic. Which version of Ubuntu (14, 16 or 18) is suitable for installing zoo-project?


texts: GeoServer cutting off symbols close to tile edges?

texts:
When I use specific icons for point display using GeoServer, several icons are cropped automatically. This is may be because of tile size or else.

How can I resolve this?


texts:
I have an excel file containing coordinates of taxi positions which I would like to turn into lines in QGIS. The example data are as follows
no   trj_id       lat        lng
1    10162    -6.32044    106.837
2    10162    -6.32043    106.837
3    10162    -6.32048    106.837
4    10162    -6.32054    106.8369
.      .          .           .
.      .          .           .
.      .          .           .
.      .          .           .
438  10165     -6.25062   106.9883
439  10165     -6.25068   106.9881
.      .          .           .
.      .          .           .
.      .          .           .

This goes on for another 377000 points with total of 471 trj_id
I have attached the points represented in QGIS as shown below

The idea is to create a line based on the same trj_id in QGIS, but I am not sure how to implement it simultaneously. I previously separate the data points based on trj_id but that would take a lot files. Is there anyway to do this?


texts: Checking to see if database is there and delete if so using ArcPy?

texts:
I am looking for my python code to check for the polycheck.gdb and if it is there delete and create a new one. Also, I would like to set my workspace to that database. Not sure if this is the right way to go about it.

# Import arcpy module
import arcpy
import os
import sys

from arcpy import env

arcpy.env.workspace = "C:/PolyCheck.gdb"

if arcpy.Exists("PolyCheck.gdb"):
    arcpy.Delete_management("PolyCheck.gdb")

arcpy.AddMessage(arcpy.GetMessages())

arcpy.CreateFileGDB_management("C:/output", "PolyCheck.gdb")

# Script arguments

I get this error:

  ERROR 000732: File GDB Location: Dataset C:/output does not exist or
  is not supported Failed to execute (CreateFileGDB). Failed to execute
  (Script9).


texts: Apply raster contraints on single table with multiple files

texts:
I have a raster table (pre-created using CREATE TABLE) containing multiple individual raster files that I append piece by piece to a final table. The table itself is in its own tablespace 'gisdir' and schema 'gisdata'.
This is done with the following command:

raster2pgsql -s 32615 -F -M fileX.tif -t 100x100 -T gisdir -X gisdir -a gisdata.map > query.sql

This query is run multiple times and with different SRIDs and filenames. So the table gisdata.map finally should contain multiple raster files of the same format which all have different extents, SRIDs and names. This is also the reason, why I failed to provide constraints in the query.

Now I want to be able to query for instance the raster with the name fileX.tif in a query (in QGIS), but it fails. Probably because it was not probably registered as raster.

How can I do this given that the AddRasterContraints command requires single SRID to be specified?
Or is this even impossible (storing many raster files with different SRIDs and extents in a single table)?


texts: Extracting point information from raster in R

texts:
I'd like to work with climate data available from here and extract values of certain locations (for which I have lat and lon).

This is what I've done so far:

library(SDMTools)
library(raster)

test <- raster(read.asc.gz("C:/.../grids_germany_annual_frost_days_1951_17.asc.gz"))

> test
class       : RasterLayer
dimensions  : 866, 654, 566364  (nrow, ncol, ncell)
resolution  : 1000, 1000  (x, y)
extent      : 3280415, 3934415, 5237501, 6103501  (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : in memory
names       : layer
values      : 13, 288  (min, max)

Now I figured I might be looking for the extract function but:

> extract(test, c(3634415,5603501))
     [,1]
[1,]   NA
[2,]   NA

Adding the buffer argument results in:

extract(test, c(3634415,5603501), buffer=1000)
Error in .cellValues(x, y, ...) : unused argument (buffer = 1000)

Could you explain to me what I am doing wrong and how I can extract values from the raster?


texts: How to Export Map Accurately as a Raster Image in QGS 3.4

texts:
I'm struggling to accurately export a map I have using QGIS. I am using the Import/Export > Export Map to Image function. My map is based on a 600x800 pixel raster on which I have overlayed other stuff. I want to export the map as a 600x800 pixel image to exactly cover what I have created and to preserve the original GeoTIFF resolution.

In the Save Map as Image dialog I select "Calculate from Layer" to set the correct extent. The number of pixels at this point is wrong. I lock the output aspect ratio and then try to alter the number of pixels to 600x800, but that starts to change the extent, which is not what I want. When I alter the number of pixels, I want QGIS to alter the Scale or Resolution accordingly and leave the extent to what I set it at. How do I do that? My GeoTIFF has a resolution of 50 metres per pixel, so that is what I want to be able to specify in QGIS.

If I try and change the Scale manually, that also changes the extent.

The only option that leaves the extent unchanged is altering the Resolution. But the issue here is that the increments are not fine grained enough: I get 597 or 607 pixels wide, not 600.

What is this dpi thingy? Dots per Inch? That sounds like some archaic measurement system that only the USA uses, the rest of the world has long since moved on to metric. Anyway, what does DPI mean in this context? I am not printing or creating an image that has a real physical size. My image co-ordinate system is in metres, so if I have to provide a scale & resolution, then I need to set them in metric units too.

Why is this seemingly simple operation so frustratingly complex and not providing an exact result?


texts: What is the style of mapping used by Google Traffic?

texts: I'm doing a piece of work that involves mapping traffic flow, and will be displaying it in a way that Google Traffic does. I know how to make the maps using QGIS, but is there a specific name for this type of mapping? I need some literature to refer to (i.e. similar papers that use this type of mapping), but don't know what it's called.


texts: Pansharpening Landsat 8 tiles in Grass 7

texts:
I am trying to get some nice looking RGB images for visualization purposes, mainly using as background maps, from the Landsat 8 archive. I was following the tutorial provided by Markus Neteler at
   http://courses.neteler.org/processing-landsat8-data-in-grass-gis-7/ and    http://courses.neteler.org/processing-landsat-8-data-in-grass-gis-7-rgb-composites-and-pan-sharpening/ using WinGRASS-7.0.0.1.

The steps I used so far

Import the tiles for scene LC81740382014220LGN00 using r.in.gdal

Set the computational region using

g.region rast=LC81740382014220LGN00_B1

Convert the digital number data (DN) to reflectance

i.landsat.toar input=LC81740382014220LGN00_B output=LC81740382014220LGN00_refl metfile=D:\Landsat\LC81740382014220LGN00_MTL.txt

Balance the colors

 i.colors.enhance red=LC81740382014220LGN00_refl4@LandsatImport green=LC81740382014220LGN00_refl3@LandsatImport blue=LC81740382014220LGN00_refl2@LandsatImport

Display

d.rgb red=LC81740382014220LGN00_refl4 green=red=LC81740382014220LGN00_refl3 blue=red=LC81740382014220LGN00_refl2

which gives already a nice looking composite. But further on

i.pansharpen red=LC81740382014220LGN00_refl4@LandsatImport green=LC81740382014220LGN00_refl3@LandsatImport blue=LC81740382014220LGN00_refl2@LandsatImport pan=LC81740382014220LGN00_refl8@LandsatImport output=LC81740382014220LGN00_copy_brovey method=brovey

and then again

d.rgb red=LC81740382014220LGN00_brovey_red green=LC81740382014220LGN00_brovey_green blue=LC81740382014220LGN00_brovey_blue

results in reddisch rgb composite

balancing the colors

i.colors.enhance red=LC81740382014220LGN00_brovey_red@LandsatImport green=LC81740382014220LGN00_brovey_green@LandsatImport blue=LC81740382014220LGN00_brovey_blue@LandsatImport

results in

Any ideas where something went wrong?

In addition: I also tried installing the i.fusion.hp extension but failed. My Grass installation is located in C:\grass70 and i extracted the archive to C:\i.fusion.hpf-master . After switching to the extension folder running

make MODULE_TOPDIR=$GISBASE
Makefile:7: ISBASE/include/Make/Script.make: No such file or directory
Makefile:8: ISBASE/include/Make/Python.make: No such file or directory
make: *** No rule to make target
`ISBASE/include/Make/Python.make'.  Stop.

Where does the G in GISBASE go?


texts: Error 'GDAL extension not installed' while installing cartaro

texts:
I am trying to install Cartaro on my Ubuntu 12.04 LTS. I have followed every step from this guide: http://cartaro.org/documentation/installation

But on the 'Verify Requirements' page I get this warning which says ''GDAL extension not installed'. Also when I click on proceed with installation, on Setup database page I see mysql and SQLite as options type not PostgreSQL as shown here

Here is what I did:

After using

sudo unzip geoserver-2.3.5-war.zip geoserver.war -d /var/lib/tomcat7/webapps/

I opened /etc/default/tomcat7 and added the following code in the end

JAVA_OPTS="${JAVA_OPTS} -DENABLE_JSONP=true"

Then I edited geoserver/WEB-INF/web.xml and added the following code just after web-app starting tag.

  <session-config>
    <session-timeout>33333</session-timeout>
  </session-config>

Please tell me what did I do wrong and what is the correct way to install the GDAL extension.

Thanks


texts:
I am reviewing the geodata quality standards within my organization. During my review I keep coming across the ISO 19100 series. I am wondering if anyone has used these standards and if they have any input on their usefulness. The purpose of my review is as follows:

identify current organizational standards
identify potential short comings of my organizations standards
implement more robust geospatial data quality standard
implement/develop robust geospatial data quality process(es)
implement more current metadata and documentation standards

The organization is responsible for maintaining authoritative geospatial data, specifically:

parcel fabric
topographic features
ortho imagery
stereo imagery

Would anyone recommend venturing into the ISO 19100 standard for such a dataset, or would other standards be more appropriate?


texts: I need a library/framework that i can use to store a custom map(as tiles) locally on android and show it as a map . I would like to get a real-time position of an object(Point data - it's within a building ) from the server ( a compatible server of course) and project it on my app . My question is - Is there a library - Server combination that can do this . I'm still reading and learning about different types of servers and libraries, so can you please be elaborate


texts: Unable to connect to database from ArcMap

texts:
I am trying to connect to SQL Server database from  ArcMap from the same instance. I am able to connect one database; to another database I am unable to connect.

I am getting this error while connecting:

  Failed to connect to specified server. Do you want to
  proceed? Operation failed.

What might be wrong?


texts: wireless sensor network integration with GIS

texts:
I am working on a project in which I need to represent data on a map after collecting it from wireless sensors.The wireless sensors are placed in the field and it collects information regarding soil ph,moisture,temperature.It is fitted with a gps which takes the co-ordinates.These attributes are then transferred to a computer in the form of an excel sheet.

I want to know if it is possible to take this data directly into Qgis via  a pen-drive,bluetooth etc and plot it on the corresponding map and show the co-ordinates from where the data has been collected


texts:
I am experimenting with table service, the first restriction I thought about is the URL query length, as it work by http get request.

I need to send 5000 location, considering location is 35 characters long as for example:

  35.921870470047,31.949649689353226;

thus, I will be sending 175000 character in the url, which is not possible I think for all browsers.

I will pass the result of durations matrix into the vehicle routing problem solver of Google Optimization Tools with the duarion matrix of a 5000 input location.

Here is an example of my query for small number of locations:

http://router.project-osrm.org/table/v1/driving/35.921870470047,31.949649689353226;35.88357925415039,31.974007590177635;35.92055082321167,31.948830365146534

Any advice?


texts:
I would like to scrape geodata off three websites, which prove quite difficult for me.

https://www.aldi-suisse.ch/filialen/ - Here I would like to get all stores locations from the map. I checked all the network traffic in Chrome, but there seems to be no location data sent.
https://map.geo.admin.ch/?lang=en&selectedNode=&topic=inspire&layers=ch.astra.unfaelle-personenschaeden_alle&bgLayer=ch.swisstopo.pixelkarte-farbe&X=160125.00&Y=554500.00&zoom=1&layers_timestamp= - This one seems by far tougher. From what I can see, they are rendering the location markers server side into pictures, then overlay the maps with those pictures (no location data is sent at all). API calls are based on individual click locations on the map, making it impossible to click all locations.

Any ideas for either of those websites?


texts: How to overlay a flat image onto google maps with two points of reference?

texts:
The Android API allows for images to be superimposed onto a map.  To do so, the developer specifies a set of lat/long coordinates for an "anchor point" within the image, the width (or height) of the image in meters, and "bearing" value by which the image will be rotated before it is superimposed.  More info here:

http://developer.android.com/reference/com/google/android/gms/maps/model/GroundOverlay.html

Instead of an anchor point + width + bearing what I have instead is an image and two points of reference inside it, with coordinates for each.  So, for example, I might have a 1000 x 1000 pixel image and know that pixel (100,100) corresponds to coordinates (X,Y) and that pixel (600,600) corresponds to (X',Y').

What I'm unclear on is how to correctly calculate the anchor point, bearing and "width" values to plug into the Android maps API such that (100,100) in my image appears at (X,Y) on the map and (600,600) in my image appears at (X',Y') on the map.  Found this page:

http://www.movable-type.co.uk/scripts/latlong.html

with info on how to calculate the distance between two points and both starting and ending bearing.  Using the distance formula I can calculate the actual physical distance between my two points of reference.  Using basic geometry I can calculate the "pixel distance" between the two points in my image.  Taken together, this should give me a conversion factor between pixels and meters, so I should be able to calculate the "width" of the image.  Is that right?

But any "flat" image is going to be a projection of some sort, though, so I'm worried that this method of computing "width" won't be entirely accurate.  If that inaccuracy is small then I can live with it.  The actual physical size of the images I'm planning to overlay shouldn't be more than 20 miles x 20 miles max.

Then there's the issue of "bearing".  I can compute the forward azimuth between (X,Y) and (X',Y') and I can compute the angle between (100,100) and (600,600) in my image.  If I choose (X,Y) as my anchor and then subtract the latter angle from the former angle then that should give me the amount by which my image needs to be rotated before being its superimposed onto the map.  Right?

Again, though, I'm worried about curvature issues.  Is that a legitimate concern given I'm only talking about (approximately) 20 x 20 mile overlay images?


texts: Insert Geometry and update attribute using arcpy.da.InsertCursor

texts:
I'm trying to insert geometries from one feature class into another and then update an attribute associated with the geometry. Only inserting the geometry works fine, but I can't work out how to insert the new value associated with the geometry. Below insertRow() gives 'TypeError: cannot read geometry sequence'.

import arcpy

aLayer = arcpy.GetParameterAsText(0)
aValue = arcpy.GetParameterAsText(1)
target = "C:\\local.gdb\\target"
desc = arcpy.Describe(aLayer)
type = desc.shapeType

if type == "Polygon":
    with arcpy.da.InsertCursor(target, ["SHAPE@","aField"]) as iCursor:
        with arcpy.da.SearchCursor(aLayer, ["SHAPE@"]) as sCursor:
            for row in sCursor:
                iCursor.insertRow([row,aValue])


texts: NetCDF to Shapefile with Python and GDAL?

texts:
I need to save my outputs, which are netcdf files or panda python tables, to shapefiles. They cannot be rasters because the software I am loading them into accepts shapefiles only.

I have googled and looked around and couldn't find some clear answers. I understood that I could call GDAL from python, and I am digging into GDAL.
I know I have to reclassify the netcdf files before saving it to shapefiles. This question is not about which best shapefile is representative of a netcdf (countour, lines etc) but what is the correct script/command to use in GDAL/python to do something like:

convert netcdf_data.nc to shapefile_data.shp

I could use some more specific input.


texts: How to delete little geometries without hole in python script

texts:
I'm writing a python standalone script from a old PyQGIS script.
Input data is a raster of which I would like to delete little areas.
In the PyQGIS script, I vectorized it and used GRASS v.clean to delete/resorb little areas/geometries, without hole.

So I'm searching how to do this in my standalone script. Either directly from the raster or from the vectorized one.

Is anybody has an idea?

input 1 band raster

output


texts: How to display pop up image when click point in OpenLayers web map

texts:
I have basic WMS web map showing parcels and points layers which representing buildings of parcel. Then I want to add building images for each and every point. Then when a user clicks on a point image should appear as a pop up window.

I am using  OpenLayers 2.10 Beginner's Guide book and can't find any method to do the above thing.

Is this possible with OpenLayers? If so please give me some advise or some links to help me.

All layer are in WMS form and I am using GeoServer and Tomcat.


texts: Using layers attribute within QgsPropertyOverrideButton

texts:
In addition to my previous question, I am struggling to implement the possibility to choose an attribute from a layer within QgsPropertyOverrideButton.

I have a combobox where one can choose a layer to use. Getting attributes and geometries from this selected layer works just fine, but it is not recognized by QgsPropertyOverrideButton when I use .setVectorLayer() as described in docs. For example self.dlg.Isochrones_WalkSpeed_Override.setVectorLayer(selectedLayer)

A reduced code snipped:

layers = QgsProject.instance().layerTreeRoot().children() # Fetch the currently loaded layers
selectedLayerIndex = self.dlg.Isochrones_SelectInputLayer.currentIndex() # Fetch the selected layer in combobox
selectedLayer = layers[selectedLayerIndex].layer() # Use the in combobox selected layer
fieldnames = [field.name() for field in selectedLayer.fields()] # Receive fieldnames from selected layer
features = selectedLayer.getFeatures() # <-- Works just fine

self.dlg.Isochrones_WalkSpeed_Override.setVectorLayer(selectedLayer) # <-- Seems to not work
self.dlg.Isochrones_WalkSpeed_Override.updateFieldLists()   # <-- Seems to not work

# Do some stuff here

for feature in features:
    # retrieve every feature with its geometry and attributes
    print("Feature ID: ", feature.id())
    # fetch geometry
    # show some information about the feature geometry
    geom = feature.geometry()
    pointgeom = geom.asPoint() #Read Point geometry
    x = pointgeom.x() #Read X-Value
    y = pointgeom.y() #Read Y-Value
    print("PointX: ", x, " | PointY: ", y)

    # fetch attributes
    attrs = feature.attributes()
    # attrs is a list. It contains all the attribute values of this feature
    print(attrs)

    #Check where to gather attributes from: GUI or Layer?
    #WalkSpeed
    if self.dlg.Isochrones_WalkSpeed_Override.isActive() == True:
        Isochrones_WalkSpeed, IrrelevantSuccessStorage = self.dlg.Isochrones_WalkSpeed_Override.toProperty().value(QgsExpressionContext()) #Receiving Value from GUI: DataDefinedOverride
    else:
        Isochrones_WalkSpeed = self.dlg.Isochrones_WalkSpeed.value() #Receiving Value from GUI: SpinBox

    # Do more stuff here

Why does .getFeatures() recognize the layer but .setVectorLayer() doesnt? And how would I need to write my code to make .setVectorLayer() recognizing the selected layer, so users are able to choose an attribute from the selected layer within QgsPropertyOverrideButton?


texts:
In what coordinate system and unit measures the coordinates returned by the Google Maps API are?
I ask this because I have a page that returns the geocode lat / long for me to find in my postgis basis. But the query does not return anything, I need to know which projection should I work with my data.


texts: Giving the same color to different shapefiles in QGIS

texts:
I have two shapefile consist of many dots.  First shapefile's dot values are between -19.060 and 9.680, second shapefile's dot values are between -13.153 and 11.758.

I want to apply the same color ramp to shapefiles. I want the values of the two shapefiles to appear on the same scale. For example, -10 should be the same color for both layers. I had many attempts at how to do this, but I was not successful.


texts: ArcGIS Copy Features tool extremely slow when exporting to ArcSDE?

texts:
I have about 10 points in a file gdb that I am attempting to copy over to an sde feature set. In arcpy and in model builder it takes about a minute to execute. Do you know of any tips to speed this up? In SSMS I can query the table almost instantaneously, so there must be a bottleneck somewhere that I am not aware of, since the network is close to zero latency. In the progress bar, it reads "Writing Features 1" for the majority of the time, then seems to complete the rest of the task almost immediately.

If it is really supposed to be this slow, I would at least like to have some idea of what it is actually doing during this time. The file is only a few kilobytes total, so I assume this cannot be normal.

Update: I was trying this operation over a site-to-site VPN. It is a fast connection, so there is still no real reason it should take a full minute to upload 10 points. I tried this on a local connection (local to the SQL Server), and the same operation took 2 seconds. It seems that CopyFeatures has some sort of bug when working through a VPN connection. Any ideas on how I could find a workaround (I need to be able to copyfeatures over the network with a VPN connection)?


texts: File Geodatabase API .NET Wrapper: How to determine if Table/Feature Class Is Locked?

texts:
I have been using the file gdb API .NET wrapper, and it has suited my needs quite nicely (vs. using ArcObjects).  One of the things I have yet to figure out is if there is a proper way to determine if a table or feature class (fc) is locked before attempting to write to it.

What I have been doing is wrapping a update/insert in a try/catch, and if I get to 'can't acquire lock exception...', I know the object is locked by something else.

Am I missing a better way to tell if a given table or fc can be written to?


texts: Permanently sourcing data to an .mxd that was previously sourced to a .mpk

texts: I had to do some work away from the server that houses our GIS data. I created a Map Package so that I could take in on the road with me. Everything worked great away from the office. When I returned I open the .mpk and saved as an .mxd to the server and resourced all the layers to the appropriate geodatabase data on the server. Everything seemed fine, so I saved the .mxd and closed Arcmap. Unfortunately, the resourced data does not stay and I have to resource it every time that I open the .mxd. I even tried deleting the layer and adding the .lyr files again. this doesn't help.


texts:
I am trying to make a Stack in QGIS from the NDVI of several Landsat 8 images. But the following error appears.

How can I solve it?


texts: Problem while cloud masking in MCD19A2.006 data

texts:
I need to create a time series of aerosol optical depth using MCD19A2.006 data. I used the 'AOD_QA' band for masking the cloud. However, the following code is not working. If someone could assist me with the issues, I would be grateful.
var geom = ee.Geometry.Point(85.14,27.7172).buffer(10000);
function bitwiseExtract(value, fromBit, toBit) {
  if (toBit === undefined) toBit = fromBit
  var maskSize = ee.Number(0).add(toBit).subtract(fromBit)
  var mask = ee.Number(0).leftShift(maskSize).subtract(0)
  return value.rightShift(fromBit).bitwiseAnd(mask)
}
var image = ee.ImageCollection("MODIS/006/MCD19A2_GRANULES").first()

var qa = image.select('AOD_QA')
var Cloud_mask = bitwiseExtract(qa, 0, 2)
var mask =  Cloud_mask.eq(1) // Clear
var maskedImage = image.updateMask(mask)

var collection = ee.ImageCollection('MODIS/006/MCD19A2_GRANULES')
.filterDate(ee.Date('2019-01-01'), ee.Date('2021-7-31'));

var collection1 = collection.filterBounds(geom)
    .select('Optical_Depth_055y')
    .map(maskedImage)
       .map(function(img) {
           return img.multiply(0.001)
          .set('system:time_start', img.get('system:time_start'));
    });


texts: How to generate polygons upto the edge of other polygons?

texts:
I am doing some work on a pipe network (blue lines) and I am trying to attribute the information about the street to the pipes.  I have 2 polygon layers one of the blocks (green hatch) and one of intersections (yellow hatch). I am trying to create a set of buffers/polygons around the street segments (dark gray lines) that fill the space not already covered by either of the 2 other layers(white space).  It would be nice to get the double wide streets to have their polygons meet at the mid point between the 2 segments.
I have tried creating a buffer and erase what falls underneath the other two layers and with an over sized buffer but I end up with spill over into the cross streets.   I need to do this for an entire city, so I can’t just make an undersized buffer and snaping to the other layers never gets far enough and then manually adjust the vertices takes too long.  Plus with both methods I end up with overlapping polygons when it buffers around a double street segment which causes problems assigning the street information to the pipe segments.
What would be the best way to accomplish this?
Let me know if you would like more clarification.


texts: Problem Installing QGIS to Samsung Galaxy Tab3

texts:
I try to install stable and nightly(18/09/2014 last nightly) versions of QGIS to Samsung Galaxy Tab3 (Android) but unfortunately I can't. When the installation is completed asked me to confirm the installation of the necessary libraries (Ministro). I tried to open the program showed the following error "The application Ministro cannot satisfy the requirements of the application".
Any help or ideas for what is the problem or what I can do.
It is the same problem like this bug report http://hub.qgis.org/issues/6497


texts: How to detect roads from satellite images?

texts: OpenStreetMap has announced a Bing API to detect roads. I want to know what is the theory behind road detection and if there is any open source tool that do that.


texts: How to periodically refresh QGIS display with updated shapefiles?

texts:
I'm creating updated shapefiles for a fixed set of layers in a project every 15 mins. If the user has that project open in QGIS on a Windows 7 PC, I can't overwrite those shapefiles, so I rename the parent directory and recreate it with the new files. But the QGIS app is still looking at the old shapefiles. A manual refresh works fine, but we want the display to update automatically without user intervention.

Tried writing a plug-in to refresh the display every few minutes, but the display doesn't update until the plug-in completes. Also tried converting the plug-in to a stand-alone app, but seemed to have the same behavior. Someone else tried it on a Mac and it seemed to work a little better, but still had issues.

Now considering a hack to call a refresh app from a Windows task - are we missing something?

Using QGIS 1.7.4 originally, also tried QGIS 2.4 but no improvement.


texts: Format .csv or .txt file so that the following stops happening?

texts:
I have a large dataset compiled of point data. Previously, this data was contained in 3 different files that all were successfully made into layers on QGIS, however I just combined them into one large spreadsheet for ease of operations. I am attempting to load it onto QGIS so I can link a buttload of attribute data and perform further queries. The problem is: QGIS has been grouping my 4th column heading (the very important Y coordinate!) with the first line of data. [What appears in the "create a layer from a delimited text file window]

Here is what my data looks like. Clearly these 2 values are in seperate columns.

I have tried :
-csv, txt and xlsx formats/ all of the different options offered in the window/ adding commas to every cell, adding that as a csv and txt file/ reformatting all cells (number, custom, text)/ adding this data as vector layer, adding as shapefile/ changing column headings.

I'm by no means an expert at GIS but this seems like a silly glitch to get caught up on when everything else had been going so smooth! What am I missing here?


texts: Running Powershell Script from QGIS Actions or Python Console

texts:
Is anyone able to explain how the Windows/Open Actions work?

And also the python console?

I have a powerscript .ps1 file I was wanting to launch from either of these options and it just doesn't seem to work.

If I try to just run the following as a Windows Action.

powershell

I get the error message "Unable to run command powershell".

However if I try to run cmd I don't get an error and I know it is doing something in the background because running the following Windows Action launches notepad.

cmd /c ("notepad")

In the python console if I run the following script then notepad happily pops up.

import subprocess
subprocess.call(['notepad'])

But if I replace 'notepad' with 'powershell' I get the following error.

Traceback (most recent call last):
File "C:\OSGEO4~1\apps\Python36\lib\code.py", line 91, in runcode
exec(code, self.locals)
File "<input>", line 1, in <module>
File "C:\OSGEO4~1\apps\Python36\lib\subprocess.py", line 267, in call
with Popen(*popenargs, **kwargs) as p:
File "C:\OSGEO4~1\apps\Python36\lib\subprocess.py", line 707, in __init__
restore_signals, start_new_session)
File "C:\OSGEO4~1\apps\Python36\lib\subprocess.py", line 990, in
_execute_child
startupinfo)
FileNotFoundError: [WinError 2] The system cannot find the file specified

I expect this has something to do with an OSGEO4 shell that is running in the background or that is called by these commands?

I haven't been able to find detailed information about or good examples of how these functions work outside of doing things within QGIS, where I want to use QGIS to run certain functions outside of QGIS (via powershell scripts), I assume it is possible?

Also if someone does run powershell scripts from within QGIS could you please explain how it is done? And if you feed it variables I'd be very interested in that too (though assuming I can run ps1 scripts I assume the easiest way would be to output the variable into a text file and just have the ps1 script read in the text in there, or output a ps1 file and call it from another ps1 file?


texts: Converting classified raster map back into vector classes using QGIS?

texts:
I have this map (have whole country, this is just a sample, original map also contains color symbology for the value classes) and need to attach the values from it to points.

I have no problem to georeference it in QGIS but there is problem with the values.

I see two uncomfortable possibilities:

attach the values manually using the map legend.
manually vectorize the map in QGIS and add the value classes

But it would be nice, if I could somehow (semi-)automatically vectorize the map using the color (RGB?) values to build the vector classes in shapefile format. Is this possible?

I have here SAGA GIS, QGIS and I have also access to one PC with Arcgis 10 with Spatial and 3D analyst. Is it possible or do I have to use the manual way?

It is almost impossible to obtain source data (will try, but I am pessimistic in this way)


texts: I have some maps in 2 files format .ozf2 + .map that I can display in QGIS. However, when I open the same QGIS project in QField, they are not displayed.


texts: Creating layer from in-memory shape data using ArcObjects?

texts:
In my current project, a feature layer is created based on shape data generated after some complicated surface analysis procedure.

The sequence right now is like
[Surface Analysis Algorithm] -> [array of xyz data] -> [shp file written to disk] -> [create feature class from shape] -> [create feature layer from feature class]

I'm on the process of optimizing the project, so I wondered if there was any way to display directly the in memory shape data in a layer, bypassing the read/write to hard disk?


texts: Is there a DWG geometry type that cannot be expressed using WKT?

texts: Is there a geometry type (both 2D and 3D) that is supported by DWG and cannot (lossless-ly) be expressed using Well-Known Text? If there is, please provide an example.


texts: Get points, along a polyline, at defined distance

texts:
I need to store into an array the points, along a polyline, at each defined distance.

Example:
A polyline length is 50m, I set the interval to 5m. I have to store in the array all the points along the polyline each 5 meter. (points on 0m-5m-10m-15m....).

I'm using C#.
Any function to get these points?


texts: What raster is needed for ArcMap pan-sharpening?

texts:
I have a "private" satellite image I would like to pansharpen in ArcMap. However, I cannot check the box in Symbology - Show - RGB Composite and in the Image Analysis window the symbol stay grey.
So I wondered if I might have the wrong tif?
What kind of TIFF is needed and how can I find out if I have the right one?
Or what else could be my problem?