texts: Computing global DEM topographic indices and compromising projections?

texts:
I would like to compute global topographic indices based on SRTM elevation data. These indices are computed along river/channel networks extracted from a DEM and are a functions of both local slope and upstream area obtained from a flow accumulation routine. For those that are familiar them, the indices I am interested in are a stream power and channel steepness index.

To compute these on a global grid (worldwide from ca. 60°N to 60°S) and to be able to compare slopes and areas in different regions of the world, I was thinking about using compromise projections such as a Dymaxion/Fuller map, which would limit both area and angle distortions. These would allow me to treat the earth as a whole and ease the processing. However I have a hard time getting a handle on the distortions these kind of projection induce and the associated uncertainty or bias in slope/area (I would be happy with less than a few percent).

Alternatives could be to split the map into small regions where UTM projections would provide good alternatives, but this would mean a more complicated workflow. Another alternative is to go for a conformal projection to get the slopes right and then run a flow accumulation routine that is weighted by a correction grid which would account for changes in area with changes in latitude.

Does anyone have comments or suggestions about these workflows?

As these global estimates are computationally intense I'd better think twice about the right projection to use.


texts: How to view ArcGIS Online-created web mapping application on my device?

texts:
I have an ESRI Organizational account and an iphone 4s (iOS 7.1.1).  On my phone using the ArcGIS app, I am able to see all maps created in My Groups, but I cannot see a single web mapping application.  All applications have been shared to the group that I am a member of just as the maps are.  I can see them all just fine from my desktop using ArcGIS Online.

Any ideas why they are not shown on my phone?  What am I missing? How do I view my applications on my device?


texts: How to set the icon for a Leaflet marker based on a variable

texts:
I'm just beginning to learn coding in d3 and leaflet, i have some code that i'm experimenting with and facing some problems with setting the icon.

I wanted to change the leaflet icon depends on the data value, but i'm stuck how to show it.

This is my code,

var ActiveIcon = L.Icon.extend({
    options: {
        iconUrl: 'images/marker-icon_green.png',
        iconSize: [17, 30]
    } });
var NotActiveIcon = L.Icon.extend({
     options: {
      iconUrl: 'images/marker-icon_blue.png',
      iconSize: [17, 30]   } });

d3.json("data_db.php", function(error, data) {
  // add map markers to map layer
   data.forEach( function(d,i) {
    // create a map marker if the lat lng is present
    if (d.latitude!=null && d.latitude!=undefined) {
      var title = d.site,
          mark = L.marker(([d.latitude, d.longitude]), {title: title});
          mark.bindPopup("<strong>" + d.site + "</strong>Status : " + d.status);
      markersLayer.addLayer(mark);
      clusterLayer.addLayer(mark);
    }
});

How can i change the icon from the data value i have, if my d.status is equal to "Active" change as green icon, and Inactive as blue icon?


texts: Creating a preview of a data set from the arcpy command line

texts:
I am creating an ArcPy tool that will be run from the commandline, and I would like the user to be able to easily create a map of results. This would essentially be equivalent to the "Preview" window in ArcCatalog.

Is there any simple way of creating such a plot while using ArcPy outside of the ArcGIS Pro GUI? I have searched through ESRI's online documentation but could not find a clear answer.


texts: GFS SIG/SFC file format of siganl.gdas/sfcanl.gdas analyzing/decoding problem

texts:
I am creating a database from datafiles. I encountered many fails and errors, which I have already fixed (also with the help of this Q&A web). Now I am stucked a bit in decoding of 6 files because I totally don't know in what format they are coded.

Those files are:

sigf06.gdas.2010123100
sigf00.gdas.2010123100
siganl.gdas.2010123100

sfcf06.gdas.2010123100
sfcf00.gdas.2010123100
sfcanl.gdas.2010123100

They are binary and when I open them in hexa (in console with od -c "filename"), I see the first 3 files have

  "\0  \0  \0       G   F   S       S   I   G      \0 003 017   ="

on the first line and the next 3 files have

  "\0  \0  \0       G   F   S       S   F   C      \0 003 017   ="

on the first line and after that all those 6 files have binary data.

All I know is that I downloaded them from directory
http://nomads.ncdc.noaa.gov/modeldata/cmd_HIC/2010/201012/20101231/ (any of the 4 .tar.gz files)
and the link to this directory I have from
http://nomads.ncdc.noaa.gov/data.php#CFSR-data -> High Resolution Initial Conditions HTTP link.
All the datafiles contain geo data in a grid, I think that those 6 files store land values.
First 3 files have 227MB, the other 3 files have 117MB.

I am already quite insane how non-documented the files are, but it's still the best source for the needed data to my database.
There is probably some documentation of the files, but I haven't found any. There was about 20 dataset files in those .tar.gzs, the majority I already decoded. If I knew just the format of these files, I would convert them easily, but link to the documentation of these files with structure and filetypes would be the best.


texts:
I am trying to create markers with text-labels inside in Tilemill. As I don't want close markers to overlap I have the

marker-allow-overlap:false;

Now, this does not affect the text labels of the markers, so some text is lying unmotivated in plain nowhere.

Labels for hiden markers are still visible

When setting text-allow-overlap to false the labels are now respecting the overlap but not visible on the markers anymore.

Is there a way to get markers with text-labels where the overlap-function is working on both simultaneously?

My full CartoCSS reference:

#tmp {
  marker-width:20;
  marker-fill:#f45;
  marker-line-color:#813;
  marker-allow-overlap:false;
  marker-ignore-placement:false;
  ::labels{
    text-name:"[values]";
    text-face-name:"Arial Bold";
    text-allow-overlap:true;
  }
}


texts: Geopandas Polygon to matplotlib patches Polygon conversion

texts:
Unfortunately geopandas plotting is extremely slow and takes a lot of resources, hence I would like to use instead matplotlib for plotting.

When I use pure Fiona to open and read the shapefile I have no trouble to extract the Polygons as matplotlib patches but now I would like to use as starting point the geopandas dataframe to get my matplotlib Polygons.

I am currently using something like:

with FI.open(df_map_elements, 'r') as layer:
    for element in layer:
        key = int(element['id'])
        if key not in dict_mapindex_mpl_polygon.keys():
            dict_mapindex_mpl_polygon[key]=[]
        for tp in element['geometry']['coordinates']:
            q = np.array(tp)
            polygon = Polygon(q) # matplotlib Polygon NOT Shapely

For plotting polygons with matplotlib:

from matplotlib import pyplot as plt
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection


texts: Geodatabase doesn't export all features into XML workspace

texts:
I try to export the database into XML workspace and when I import the XML into a new geodatabase I find that some features of latest update are missing (at least what I have identified).
I'm facing this issue only in this database. I'm using SQL Server Express + ArcGIS 10.5.
I have some basic database management skills, I'm not able to identify if it is a SQL Server issue or an ArcGIS issue.
Has anybody have had similar issue?


texts: ArcGIS Network Analyst OD Cost Matrix - how to change network?

texts:
I am using Network Analyst OD Cost Matrix to calculate the distance between a number of locations in dataA and a number of locations in dataB.

I wonder how to change the network used for these calculations. The default under Network Analysis -> Network Data Source is arcgis.com. I assume that the Input Network Analysis Layer is where I choose along which network the distances should be calculated. By default, it says OD Cost Matrix.

I have my own road map that I would like to calculate the distances along, in shp-format. I have added this to my map, so I can see the roads. But I cannot figure out how to use this as th network in OD Cost Matrix calculations.


texts: Detecting automatically where rivers flow, directly into the sea or into another river

texts:
I am using this dataset of rivers in France. I'd like to automatically distinguish between rivers flowing directly into the sea versus rivers that flow into another river and draw a map of only those rivers. How can I do this?
Each river in this dataset has a standard SANDRE code and can be looked up online (example), but they don't seem to record the information I need. I don't mind having to write code to do this, but I would like something reliable (i.e. no "guesswork"). For example, is there a way of knowing whether an endpoint of a river is on another river? If so, I could look for rivers in which neither endpoint is on a river. However, from what I can tell shapefiles don't store whether a vertex is "on" another line.
I'm open to any kind of approach. I'd like to use this dataset though since I think it's probably the most exhaustive (it's a French government dataset).


texts: Removing NaN value from legend in matplotlib map plotted from GeoPandas DataFrame

texts:
I have a GeoPandas GeoDataFrame which I want to plot against a specific column.
for day in days_list:
    df2 = csv_file[["state_id", "state_name", "tmx_corrl"]].loc[csv_file["days"] == day]
    new_df = df1.merge(df2, on = 'state_id', how="left")
    new_df.drop(columns=['state_name_y'], inplace = True)
    new_df.rename(columns={"state_name_x": "state_name"}, inplace = True)

    #new_df = pd.merge(df1, df2, on="state_id", how="left")
    #df_new = new_df.dropna(axis="index", how="any")
    df_new = new_df.fillna(0)
    #df_new = new_df.replace(np.nan,0)

    cmap = ListedColormap(["red", "orange", "blue", "green"])

    fig, ax = plt.subplots(figsize=(12,8))
    df_new.plot(ax=ax, linewidth=1, color='black')

    df_new.plot(
    ax=ax,
    column=new_df["tmx_corrl"],
    linewidth=1,
    cmap=cmap,
    scheme='User_Defined',
    edgecolor='black',
    alpha=0.9,
    classification_kwds={'bins':[0, 0.25, 0.50, 0.75, 1]},
    legend=True
    )

    ax.set_title("Day "+ str(day)+" Rainfall")

    plt.tight_layout()

    plt.savefig("day_{}.png".format(day), dpi=300)

This column have some NaN values, like at column no 7 in this dataframe

which I am removing using several methods

Dropping the column with NaN value

    df_new = new_df.dropna(axis="index", how="any")

Filling the NaN value to Zero

    df_new = new_df.fillna(0)

Replacing the NaN value to Zero

    df_new = new_df.replace(np.nan,0)

NaN value changed to zero

But it is still shows the NaN value in the map legend which I want to remove i.e. the interval " nan - 0.00 " in the legend, as it is altering the whole order of color coding.

Can somebody suggest, what am I missing here?


texts: Opening EOS netCDF4/HDF5 file with correct format using xarray?

texts:
I am using xarray (version 0.12.2) to work with the MODIS/Terra+Aqua MAIAC Land Aerosol Optical Depth dataset (MCD19A2.006). xarray successfully opens that dataset with the correct variables and dimensions but the dataset has projection metadata that is currently not being read. Has anyone else encountered this problem or have solutions to parsing the metadata to construct the correct coordinates?

Here is an example:

import xarray as xr

ds = xr.open_dataset('MCD19A2.A2000057.h09v07.006.2018013034454.hdf')

print(ds)

Which outputs the following dataset information:

<xarray.Dataset>
Dimensions:            (Orbits:grid1km: 1, Orbits:grid5km: 1, XDim:grid1km: 1200, XDim:grid5km: 240, YDim:grid1km: 1200, YDim:grid5km: 240)
Dimensions without coordinates: Orbits:grid1km, Orbits:grid5km, XDim:grid1km, XDim:grid5km, YDim:grid1km, YDim:grid5km
Data variables:
    Optical_Depth_047  (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    Optical_Depth_055  (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    AOD_Uncertainty    (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    FineModeFraction   (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    Column_WV          (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    AOD_QA             (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    AOD_MODEL          (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    Injection_Height   (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    cosSZA             (Orbits:grid5km, YDim:grid5km, XDim:grid5km) float32 ...
    cosVZA             (Orbits:grid5km, YDim:grid5km, XDim:grid5km) float32 ...
    RelAZ              (Orbits:grid5km, YDim:grid5km, XDim:grid5km) float32 ...
    Scattering_Angle   (Orbits:grid5km, YDim:grid5km, XDim:grid5km) float32 ...
    Glint_Angle        (Orbits:grid5km, YDim:grid5km, XDim:grid5km) float32 ...
Attributes:
    HDFEOSVersion:                     HDFEOS_V2.19
    StructMetadata.0:                  GROUP=SwathStructure\nEND_GROUP=SwathS...
    Orbit_amount:                      1
    Orbit_time_stamp:                  20000571645T
    CoreMetadata.0:                    \nGROUP                  = INVENTORYME...
    ArchiveMetadata.0:                 \nGROUP                  = ARCHIVEDMET...
    identifier_product_doi:            10.5067/MODIS/MCD19A2.006
    identifier_product_doi_authority:  http://dx.doi.org

You can see that there are no coordinates associated with dataset but in reality it does have geographic coordinates, it is just buried deep in metadata and needs to be extracted. My next step is to read the projection information from the  StructMetadata.0 attribute and construct the coordinates from there. But I wonder if I am missing something with xarray or if someone else has code to where it would extract the geographic information automatically.


texts:
how can I select only good data (Quality informations) from Modis product in Google Earth Engine.
I tried use this code:

var collection = ee.ImageCollection ('MODIS/MYD13Q1').select(['NDVI'],['SummaryQA'] == 0).filterDate(ee.Date('2013-01-01'),ee.Date('2015-12-31'))
Map.addLayer(collection)

But this does not work


texts: CartoDB named maps and events / infoWindows

texts:
I've created the named map, it's working well and I'm able to query it with setParameter() method and use some basic methods like setOpacity() and show/hide but I haven't succeeded to use any of the methods related to the events and infowindows. It appears that I can't derive a sublayer from the named map layer like this:

var sublayer = layer.getSubLayer(0);

This is what is done in all examples for creating infowindows. I tried layer.getSublayerCount() and it returns zero, so there's no way to get sublayer object out of the layer. Also none of the event methods work. Is there a workaround?


texts: Date information from point layer of fire occurrence from FIRMS VIIRS not ingesting properly in GEE

texts:
I just downloaded a Shapefile which contained point data of fire occurrence from the NASA FIRMS VIIRS database. I immediately ingest the files (.dbf, .prj, .shp and .shx) into Google Earth Engine as a table asset. After the asset has been ingested, I found that the column containing date information has been modified and is wrong. In the original file, for each point, date is stored in a QDate format (when inspected in QGIS) i.e. YYYY-MM-DD. However, when I view the Elements of this asset in GEE, I see very long integers instead in the ACQ_DATE (Long) column e.g. 1349247600000.
I am guessing this has to do with how date data can sometimes be represented as milliseconds since 1970-01-01 00:00:00. Is there a way for me to resolve this problem without manually recalculating the date information? This is because I intend to subset the data in GEE by querying specific time intervals, and thus would need this column attribute to recognized as a date and not a text string format. I have made this GEE asset public so that it can be accessed. The following is the script I am using.
//  Load VIIRS point data
var fire_raw = ee.FeatureCollection('users/earthorbgis/raw_viirs_fire_occurence_01_01_2012_to_31_12_2018');

//  Filter only to High and Nominal confidence
var fire = fire_raw.filter(ee.Filter.eq("CONFIDENCE",['h','n']));

//  Subset data for each year
var fire_2012_2013 = fire.filterDate('2012-01-01','2012-12-31');
var fire_2013_2014 = fire.filterDate('2013-01-01','2013-12-31');
var fire_2014_2015 = fire.filterDate('2014-01-01','2014-12-31');
var fire_2015_2016 = fire.filterDate('2015-01-01','2015-12-31');
var fire_2016_2017 = fire.filterDate('2016-01-01','2016-12-31');


texts:
I would like to know whether it´s possible to add or remove specified widgets on runtime.

form example i press any tool on my toolbar and add widget to sidebar .
press again to remove widget.


texts: Pass buffer as variable using arcpy.Buffer_analysis

texts:
I am trying to pass a buffer variable and seem to have run into some sort of limitation of the command.

The buffer will vary each time the script is run so I need to have this parameter variable

varBuffer = 5 # in this instance
varBufffer = '"'+ varQuan + ' Miles"'

I tried this but it does not work:

arcpy.Buffer_analysis(projectArea, MileBuff_shp, varBuffer, "FULL", "ROUND", "ALL", "")

Only this works.

arcpy.Buffer_analysis(projectArea, MileBuff_shp, "5 Miles", "FULL", "ROUND", "ALL", "")

Is there a work around for this?


texts: Looping through nested IGraphicsContainer

texts:
I'm looping through the IElement objects in an IGraphicsContainer and need to get the properties of another IElement by name.

I thought I could write a getElementByName function to loop through an IGraphicsContainer and return the IElement, but as soon as I call this function, the state of the IGraphicsContainer of the main loop changes.

How can I get an IElement by name while looping through an IGraphicsContainer?

Edit: While looping through one GraphicsContainer, I'm calling a getElementByName function which itself loops through a GraphicsContainer.  When the function calls GraphicsContainer.Reset(), the outer loop GraphicsContainer is also reset.  How do I stop this?

Edit: hnasr came up with a very good solution in VB.  Here it is in C#

using ESRI.ArcGIS.ArcMapUI;
using ESRI.ArcGIS.Carto;

namespace FooBar
{
    public class dotNetGraphicsContainer : Dictionary<string, IElement>
    {
        public void Refresh(IMxDocument pMXDoc)
        {
            base.Clear();
            IGraphicsContainer pGC = pMXDoc.ActiveView.GraphicsContainer;
            pGC.Reset();
            IElement pElement = pGC.Next();

            while (pElement != null)
            {
                IElementProperties pElementProp = (IElementProperties)pElement;

                if (!this.ContainsKey(pElementProp.Name))
                    this.Add(pElementProp.Name, pElement);

                pElement = pGC.Next();
            }
        }
    }
}


texts: Using ModelBuilder to Automatically Import Excel, Display XY and Append to Master Shapefile?

texts:
I'm trying to automate a daily workflow using ModelBuilder. I have several tasks that I do daily in this order:

Import an Excel file that contains 16 fields including two fields with x and y coordinates (xcoordinate and ycoordinate). The Excel is always stored in the same folder with date as the file name.
Map the information using the Display XY function
Append the resulting shapefile to a master shapefile

I've attempted to build a model that reflects this workflow (see attached).

The model appears to run, the resulting shapefile is blank. Upon re-runs to test, I get this message:

Traceback (most recent call last):
File "c:\program files (x86)\arcgis\desktop10.3\ArcToolbox\Scripts\ExcelToTable.py", line 254, in <module>
    arcpy.GetParameterAsText(2))
File "c:\program files (x86)\arcgis\desktop10.3\ArcToolbox\Scripts\ExcelToTable.py", line 225, in excel_to_table
    template=tmp_table)
File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\management.py", line 15242, in CreateTable
    raise e
ExecuteError: ERROR 999999: Error executing function.
The maximum record length has been exceeded.
Failed to execute (CreateTable).

Failed to execute (Excel To Table).

Outside the model, the excel imports into ArcMap without any issues. I've tried switching to the "1904" option in excel as well but it still does not run.

Any ideas on where the error is coming from?

I'm running ArcMap 10.3.


texts: Adding label to whole group of drainage features using QGIS?

texts:
I am using QGIS 3.2 Bonn and have a lot stuff imported from Google Earth.

Am I able to create a group of all drainages and label them with one click, or do I need to choose/select a new one each time?


texts: How do I store a reference to a workspace in a feature?

texts:
I have a feature class where each feature needs a reference to a raster in another workspace.

Is there some way I can get a string that uniquely identifies any workspace (ArcSDE, personal or file)?

Since I am writing from a geoprocessing tool IDEWorkspace2.ConnectionString looked promising. However I found that it returned an empty string when the workspace was a .gdb.

Another approach I thought may work was to call IWorkSpace.ConnectionProperties and then XML serialize. However this also seems to return an empty string when the workspace is a .gdb. Not sure if I did this correctly though.

Also if there is a way to get such an unique string, how do I get a IWorkspace object from this string?


texts: How to determine stream channel Reach locations?

texts:
I have a stream network and want to determine stream reach locations

In attached figure, Reach 4 is at the junction of sub2 and sub5, and Reach 3 is at the junction of sub1 and sub4. Any suggestions to determine these junctions based on confluence points using ArcPy or within ArcMap will be very much helpful.


texts: ArcGIS Pro 2.5.0 - Is there any way to snap a polygon / feature class to pixels within a Raster?

texts:
I need to clip a raster to a specific area of interest, and a specific number of pixels.
Whilst I can manually draw my AOI feature class and get the right number of pixels, it generally contains a number of partial pixels, which causes issues with later analysis (image attached).
Is there any way to draw a polygon that snaps to the pixels of a raster image, so that no partial pixels are present within the clipped raster?
I'm using a 10m resolution Sentinel-2 image, and trying to clip to a 16 x 16km square (1600 x 1600 pixels). The image has been reprojected from WGS 84 to British National Grid.
I have tried using maintain_clipping_extent and whilst it solves the partial pixels, unfortunately the pixels are resampled (which affects the spectral values I need for later analysis) and gives me 1599 rows and columns rather than the 1600 I need.


texts: Easy Labels do not show when opening project on another machine

texts:
Love the easy label plugin, once I got the hang of it. I opened the project file on my desktop after setting up all the labels on my laptop and they do not show. I assume they will once I open the Project on my laptop again.

Thoughts?


texts: How do I limit the number of rows that are loaded in a QGIS attribute table?

texts:
Every time I accidentally open an attribute table for a large table in QGIS it takes ages to load. The data is in Postgres, but that doesn't help when the table is trying to load millions of features. Is there an option to only load say the first 1000 rows, then load the rest as required, like pagination? Similar to ArcGIS.

I know that I can query first, or use psql or pgadmin, but sometimes you just want to see the attributes in QGIS. Sometimes it's an accident.


texts: Deleting lines crossing particular features from another layer using QGIS

texts:
I have two QGIS layers. Layer 1 has a floor plan with an attribute field holding feature types, e.g., wall, corridor etc. Layer 2 consists only the lines drawn between door centroids. I want to query the lines in Layer 1, which cross a wall or corridor from Layer 2, then delete these lines.
How would this be done using QGIS?
This is what I have:

and this is what I want to have:


texts:
Using ArcObjects SDK for Java, I have the full path of 2 shapefiles: String in1, String in2, and I want to merge them using the geoprocessing tool Merge.

I tried to several methods:

IGpValueTableObject vtObject = new GpValueTableObject();
vtObject.setColumns(2);
vtObject.setValue(0, 0, in1);
vtObject.setValue(0, 1, 1);
vtObject.setValue(1, 0, in2);
vtObject.setValue(1, 1, 2);

or

IGpValueTableObject vtObject = new GpValueTableObject();
vtObject.setColumns(2);
vtObject.addRow("'" + in1 + "' 1");
vtObject.addRow("'" + in2 + "' 2");

or

IGpValueTableObject vtObject = new GpValueTableObject();
vtObject.setColumns(2);
vtObject.loadFromString(in1 + ";" + in2);

Before executing the the gp tool:

Merge merge = new Merge();
merge.setInputs(vtObject);
merge.setOutput(out);
IGeoProcessorResult result = gp.execute(merge, null);

In each case, I get this error: AutomationException: 0x80004005 - Unspecified error

I guess I am wrong in my way to set up the Merge inputs but I didn't find any samples explaining how to do it, except in C# or VB.


texts: Edit KML lines in Illustrator or getting smooth curves

texts:
I have a problem, when I'm making my traces in an application like Google Earth, I can't make smooth curves, only straight line.

So I thought maybe I should export the trace to Illustrator then modify the curves in it and then export them.

Is this the right way to do it ? How to do this ?


texts: subtract polygons from other polygons with SQL

texts:
I have a MSSQL table that houses a bunch of polygons. I want to select a few ( with a WHERE clause ) and subtract another bunch from these first ones.

Compare it with having all plot-outlines and building-footprints in a table.

How can i get the difference between these two sets? ( in the example, return all gardens (where a garden = plot - house)

Something along these lines: (this is obviously not working - syntax is wrong too)

SET plots = SELECT geom FROM tableA WHERE type= 'plot'
SET houses = SELECT geom FROM tableA WHERE type= 'house'

SELECT plots.STDifference(houses) AS gardens FROM tableA

I can't change the data structure, so proposing to move plots and houses into a separate table is not helpful.


texts: how to publish 15000 TIFF files on GeoServer

texts:
have ortophotos of an area that I need to serve as a baselayer of my map in openlayers, i'am using geoserver and php5.

i have 15000 files (1000MB each) in raw TIFF format with corresponding TFW word file - 15Tb of data. Coordinate system is EPSG:4326 (wgs84).

on geoserver there is "imagemosaic" option wich take the whole folder of tiffs to publish them - but it does not work for me. and i can create mosaic using QGIS, just the probleme it will take me forever to make one mosaic of these files


texts: Query 200-300 viewshed rasters from a vector layer with z-values with QGIS

texts:
What I got:

200-300 invisibility depth viewshed rasters with information about how tall buildings can be without being visible from the viewshed points
A vector layer with a planned building as a polygon with z-value (building height)
QGIS 2.18.12

What I want to do:

Query all my viewshed rasters with the planned building vector layer to find out if the planned building intersects with any of the viewshed rasters

Could anyone point out an appropriate way to solve this?


texts: How to fill in the parameters for kernel density estimation?

texts:
I am desperately trying to calculate animals home range with kernel. Unfortunately, the result is obviously incorrect and I cannot identify my mistake (I would like to post a screenshot but its not allowed to me yet).

Is there an instruction how to fill in the parameters for kernel density estimation?

Is there a concept how the attribute table has to look like? (I have three columns: ID, longitude, latitude in decimal degrees)

Using QGIS 1.8


texts: Merge shapefiles without creating extra fields

texts:
I have a road shapefile with fields like line, name, line id, line type (like nh,sh) etc.

I also have another line shapefile which is the updated one of the older one. Here one line is changed from nh to sh and vice-versa. I want to know how can I merge the new shapefile with the older one without increasing the the number of fields. I have been doing it manually.

Is there  any tool to help me out or any short-cut procedure?

Rectifying 100s of lines is taking too much time.


texts: Leaflet dialog box that contains info on 1 to n layers that were clicked on

texts:
I have a Leaflet application where the user can click on a feature and the info for that feature will be updated in dialog box.

var info = L.control();

info.onAdd = function (map) {
    this._div = L.DomUtil.create('div', 'info');
    this.update();
    return this._div;
};

info.update = function (props) {
    this._div.innerHTML = '
Property Info
' +  (props ? 'Address: '+ props.address +'
Property Type: '+props.buff+'
Buffer Distance: '+props.buff_dist: 'Click on a Property Pal');
};
info.addTo(map);

var theMarker={}
var main_layer;
map.on('click', function(e) {
    map.eachLayer(function(layer) {
        if (layer instanceof L.GeoJSON) {
            layer.eachLayer(function(l){layer.resetStyle(l);});
        };
    });
    var str = "Latitude: " + e.latlng.lat.toFixed(5) + " Longitude: " + e.latlng.lng.toFixed(5) + " Zoom Level: " + map.getZoom();
    var lat = e.latlng.lat;
    var lng = e.latlng.lng;
    ProcessClick(lat,lng,str)
});

function ProcessClick(lat,lng,str){
    if (theMarker != undefined) {
      map.removeLayer(theMarker);
    };
    theMarker=L.marker([lat,lng],{icon: pineapple_icon}).addTo(map);
    map.eachLayer(function(layer) {
        if (layer instanceof L.GeoJSON) {
            layer.eachLayer(function(layer2) {
                var click_intersection = turf.booleanPointInPolygon(theMarker.toGeoJSON(), layer2.toGeoJSON());
                if (click_intersection) {
                    layer2.setStyle({fillColor: 'lightblue', opacity: 0.8});
                    var i=layer2.feature.properties;
                    console.log(layer2.feature.properties);
                    info.update(i);
                };
            });
        }
    });
    $("#map_coords").html(str);
};

The code is setup so the user can click on multiple overlapping features but only the last feature that is clicked gets its info displayed in the dialog box on the top right.

Picture below explanation: the childcare and religious buffers layers are displayed on the map. in this example I click on an overlapping childcare and religious buffer. When I click on these features they get highlighted and the info from these two features updates the top right dialog box (would this be called a dialog box?).

What I want: A dialog box that stores all the info from the features that are clicked. A way to scroll through the results and when info from one feature is displayed in the dialog box, only that feature will be highlighted.


texts: Spatial coordinates - How to convert DMS coordinates to WGS84/UTM coordinates?

texts: Work colleagues handling on-site surveys unusually collected survey points with DMS coordinates such as : N15°58.796'. Usually they sending us points with WGS84/UTM coordinates. I would ask how I convert DMS coordinates to WGS84/UTM coordinates? And specifically to WGS84/UTM zone 20N (EPSG 32620).


texts: How to create polygons from unclosed polylines?

texts:
I have a layer of driveway that is in polylines unclosed . They are just two parallel lines outlining the width of the driveway. Then again, there is a gap (sidewalk layer) and driveway continues until it touches the road.
So my question is how to extend the line? and how to close all those lines to form a polygon?

I also have a lot boundary that has a unique ID. So I can summarize the lines that fall within the same boundary ID. But I am encountering the problem to connect the lines and make it a polygon layer.


texts:
I'm trying to measure how far along a river polyline a feature is using the measurement tool in ArcGIS. Is it possible to trace the polyline (in a similar method to the trace tool in an edit session) to get the total actual distance rather than straight line /estimated distance?

(I'm aware that you can have segments with the measurement tool, but over large distances this quickly becomes inaccurate compared to the feature length)

Currently my work around is to create a copy of layer, edit the layer, split the line and re-calculate the length which seems excessive for a quick query for a colleague.


texts: Clicking to create multiple icons on map using OpenLayers 4?

texts:
I am trying to build python-based (flask) application where multiple icons can be
placed on an openstreetmap by clicking on the various spots where I want an icon. I am trying to
achieve that using OpenLayers 4. Here is what I have tried:

var iconStyle = new ol.style.Style({
   image: new ol.style.Icon({
      size: [21, 25],
      offset: [0, 0],
    //  src: 'C:\flask_apps\click_coord\images\mark.png'
      src:'{{url_for('static', filename='mark.png')}}'
   })
});

var map = new ol.Map({

   layers: [
      new ol.layer.Tile({
        source: new ol.source.OSM()
      })
    ],
   target: 'map',
   view: new ol.View({
     center: [0, 0],
     zoom: 2
   })
 });

var vectorSource = new ol.source.Vector();

var source_name = 'iconSource';
vectorSource.set('name',source_name);

var vectorLayer = new ol.layer.Vector({
   source: vectorSource,
   style: iconStyle
});

var layer_name = 'iconLayer';

vectorLayer.set('name',layer_name);

map.addLayer(vectorLayer);

map.on('click', function(evt) {
   var position = evt.coordinate;

   var mlayers = map.getLayers();
   var length = mlayers.getLength(), l;

   for (var i = 0; i < length; i++) {
      l = mlayers.item(i);

     if(layer_name ===l.get('name')) {

       var iconFeature = new ol.Feature({
         geometry: new ol.geom.Point(ol.proj.transform([position], 'EPSG:4326',
         'EPSG:3857'))

       });

       iconFeature.setStyle(iconStyle);

       this.vectorLayer.getSource().addFeature(iconFeature);
       this.vectorLayer.clear();
       //   this.vectorLayer.redraw();
       //   this.vectorLayer.changed();
       //   this.vectorLayer.refresh({force:true});

     }
   }

 });

I have checked using view page source that the framework is resolving the path to the icon image
properly. I had previously hard-coded the location of the icon image. Irrespective of what I try, nothing
happens when I click the map. I have tried refreshing the map (I may not be doing this properly), that
doesn't seem to work either.

Could someone point out what I am doing wrong?

I am using OpenLayers 4.2.0


texts: Updating Layer tree after changing symbology with PyQGIS

texts:
When I change symbology for a given QgsVectorLayer (QGIS 2.14.16 LTR) by

r=iface.activeLayer().rendererV2()
r.symbol().symbolLayers()[0].setColor(0,0,255,255)

as shown here:

and refreshing the map content by

iface.mapCanvas().refreshAllLayers()

I get the expected result, but the Layer Tree does not update as shown here:

What am I missing, and how can I update the Layer tree from PyQGIS programmatically?


texts: How to create trigger to check geometry validity?

texts:
I want to create trigger to check geometry validity in spatial table. I have a table:

CREATE TABLE vec_geometry
(
  id bigserial
  is_valid boolean,
  invalid_reason character(128),
  geom geography
)

When I try to create trigger function I get the error:

CREATE FUNCTION geom_valid_check() RETURNS trigger AS
$BODY$
WITH test as (SELECT * FROM ST_IsValidDetail(NEW.geom))
UPDATE vec_geometry
SET is_valid = (SELECT valid FROM test),
    invalid_reason = (SELECT reason FROM test)
WHERE id = NEW.id;
RETURN NULL;
$BODY$
LANGUAGE plpgsql VOLATILE NOT LEAKPROOF;

An error has occurred:

ERROR:  syntax error at or near "WITH"
LINE 3: WITH test as (SELECT * FROM ST_IsValidDetail(NEW.geom))

How to rewrite function to make it work?


texts:
Continuing what is written in this post, I show gdalinfo output

Size is 1200, 1400
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",42],
    PARAMETER["central_meridian",12.5],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-598884.437500000000000,650790.250000000000000)
Pixel Size = (996.479919433593750,-999.677795410156250)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2018:04:01 00:10:06
  TIFFTAG_DOCUMENTNAME=/datamet/data/tmp/export/RADAR/36/RRN/VMI/VMI.geotif
  TIFFTAG_IMAGEDESCRIPTION=IDL TIFF file
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=IDL 8.3, Exelis Visual Information Solutions, Inc.
  TIFFTAG_XRESOLUTION=100
  TIFFTAG_YRESOLUTION=100
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -598884.438,  650790.250) (  4d32'13.20"E, 47d34'44.14"N)
Lower Left  ( -598884.438, -748758.664) (  5d56'18.16"E, 35d 4'36.70"N)
Upper Right (  596891.466,  650790.250) ( 20d26'12.34"E, 47d34'50.73"N)
Lower Right (  596891.466, -748758.664) ( 19d 2'23.69"E, 35d 4'40.94"N)
Center      (    -996.486,  -48984.207) ( 12d29'17.00"E, 41d33'32.31"N)
Band 1 Block=1200x1 Type=Float32, ColorInterp=Gray

When I open it using GeoTools I get the following CRS:

PROJCS["unnamed",
  GEOGCS["WGS 84",
    DATUM["World Geodetic System 1984",
      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
    UNIT["degree", 0.017453292519943295],
    AXIS["Geodetic longitude", EAST],
    AXIS["Geodetic latitude", NORTH],
    AUTHORITY["EPSG","4326"]],
  PROJECTION["Transverse_Mercator"],
  PARAMETER["central_meridian", 12.5],
  PARAMETER["latitude_of_origin", 42.0],
  PARAMETER["scale_factor", �],
  PARAMETER["false_easting", 0.0],
  PARAMETER["false_northing", 0.0],
  UNIT["m", 1.0],
  AXIS["Easting", EAST],
  AXIS["Northing", NORTH]]

I made an attempt and I get it

gdalwarp -s_srs "+proj=tmerc +lat_0=42 +lon_0=12.5 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -t_srs "EPSG:4326" source.tif target.tif

ERROR 1: Too many points (441 out of 441) failed to transform, unable to compute output bounds.
Warning 1: Unable to compute source region for output window 0,0,1200,1400, skipping.
0...10...20...30...40...50...60...70...80...90...100 - done.


texts: QGIS3 Unable To Get Local Issuer Certificate

texts:
I just launched QGIS3 (2.99) after an sudo apt-get upgrade in Ubuntu 16.04. But I get an error on startup that says:

* Unable To Get Local Issuer Certificate: The issuer certificate of a locally looked up certificate could not be found

Screenshot below:

Is it OK to ignore it or is there something on my Lubuntu that needs to be fixed?


texts: When the mouse moves on the grid, so does the corresponding point on the map

texts:
I am trying to create a specific functionality in mapping application.

The application currently lets users draw a profile line on a map and click run. The map then takes a DEM file and creates a grid below the image that shows "the lay of the land;" in other words it shows a 2 dimensional display of the elevation on an x/y axis that traces the rise and fall along the profile line that the user initially created.

What I want to do is have a ball appear on the profile line in the map that traces the mouse movements when the user's mouse cursor is on the grid (I already have crosshair, highlight ball, and elevation tooltip that uses the Flot plothover functionality).

Does anyone know of an application that already has this functionality or know of a way to achieve this; even in simplistic terms?

I am currently able to pull out the coordinates that trace the initial profile line on the map (decimal degrees), but I am not sure how to relate them to the grid coordinates (not real coordinates, basically just elevation values and values relating to points on the x-axis).


texts: Convert node data to lat/ long > DBF to CSV conversion

texts:
I have a DBF table with nodal data locating cycle injuries in Portland, OR. I'd like to convert the table to a .CSV. Transferring the attribute data isn't a problem, but I'd like to get lat/long values for the nodal values (circled in enclosed image).

What exactly does the node data describe, the roads as a network, and these points nodes in the network?

Is there a way to do this in bulk? If the answer involves using QGIS or ArcGIS, I'd prefer to use QGIS.


texts:
I have a question about Landsat 8 I am going to work on...

I am thinking about stacking all the bands however there are 15 m, 30m, and 100m

Do you stack all together ? This is to give the end-user to change their band ratio in ArcGIS.

or should I just stack first 3 band ratio for each such as

According to this website : http://blogs.esri.com/esri/arcgis/2013/07/24/band-combinations-for-landsat-8/

4,3,2

7,6,4

5,4,3

and so on

I'd like to hear your input


texts: Classifying image using image segmentation in google earth engine

texts:
I am trying to classify a Landsat 8 image in google earth engine. I have to classify them object based, since I am building on an earlier research where they have used object based classification.

Create objects has been done in the same way as in the segmentation session of the user summit https://code.earthengine.google.com/189e4388e81f01c967b72c0bd78d5d45. I also took ndvi and evi in to account.

I have ended with a bunch of objects, which I made into vectors. Is there a way that I can assign these polygons to a class?

Can I code which particular polygon would represent a particular class, so that my generated objects will become my training areas?


texts: Building a web-app/service compatible application with QGIS

texts:
I have been tasked with building an application that uses QGIS compatible geoprocessing libraries to do automated GIS tasks. The application should follow the criteria below:

Coded in python
Possible to run on a local machine
Coded in such a way that future integration as a web app or web service should be smooth
Complete with GUI features
As efficient as possible. Meaning that the functionality should be called from python in as low level a way as possible. This is why I am a bit skeptical about the iface API.

I really have no idea where to start. Should I even be thinking in terms of QGIS or should I be looking at individual geoprocessinglibraries such as GDAL or GRASS? Isn't QGIS just a platform with a mix of different geoprocessing tools? Should I be using iface?


texts: How to draw a boundary feature of an area?

texts:
Which shape file is more appropriate for drawing a boundary of an area. I mean would I take polyline feature or polygon feature for drawing a boundary?

I want to digitize a map of metropolitan area. First I have made a manual map and now I am trying to digitize.In this map there are moujas, thanas, rivers, roads, different important location such a market place, police station. So in this case, there is overlapping of borders. Because a thana border can be a district border at the same time. So In this case which feature I chose.


texts: Legend colour not matching as assigned in the bins in matplotlib map plotted from GeoPandas DataFrame

texts:
I have a GeoPandas GeoDataFrame which I am plotting against a specific column and the colour of the output/map is not matching with the colour assigned to different bins/intervals. It should ideally be unique for distinct intervals but the legends are taking same colour for different intervals
##colour

    cmap = ListedColormap(["red", "orange", "blue", "green"])

## bins/intervals

    classification_kwds={'bins':[0.25, 0.50, 0.75, 1.00]}

Same colour legend for 2 different intervals is coming in the output. Why it is happening so?

The code is as follows

for day in days_list:
    attribute = "far"
    df2 = csv_file[["state_id", "state_name", attribute]].loc[csv_file["days"] == day]
    new_df = df1.merge(df2, on = 'state_id', how="left")
    new_df.drop(columns=['state_name_y'], inplace = True)
    new_df.rename(columns={"state_name_x": "state_name"}, inplace = True)
    df_new = new_df.fillna(0)

    cmap = ListedColormap(["red", "orange", "blue", "green"])

    fig, ax = plt.subplots(figsize=(12,8))

    df_new.plot(
        ax=ax,
        column=attribute,
        linewidth=1,
        cmap=cmap,
        scheme='User_Defined',
        edgecolor='black',
        alpha=0.9,
        classification_kwds={'bins':[0.25, 0.50, 0.75, 1.00]},
        legend=True
    )

    ax.set_title("Day "+ str(day)+ " " + str(attribute))

    plt.tight_layout()

    plt.savefig("day_{}.png".format(day), dpi=300)

This is the Geopandas GeoDataFrame


texts: Allocating (or adjusting) points to closest polygon QGIS

texts:
Hi I've recently posed a question on how to create buffer areas for coastal borders of country polygons in order to match all my points (of which many are just outside country polygon borders) to a country. However, expanding country polygons through buffer zones does not provide a reliable solution.
Using higher quality global shapefiles helped a little but did not work for all my points.

Is there a way to adjust points to the closest polygon (e.g. change their spatial location to sit within the closest polygon)? Or is there another way to reliable count all points either within or next to the closest polygon? All my points have to be associated with one country, even if they sit 1 or 2km outside the polygon boundaries.

I've played around with NNJoin a little, but couldn't find a way to get it to work. Which doesn't say much as i'm a beginner.


texts: Using aggregate function in QGIS graphical modeler?

texts:
I'm looking to create five fields of grouped means within 33 shapefiles. I'm using the Graphical Modeler to save me using the Field Calculator 5 x 33 times.

The new fields are the means of the traffic volumes ("AADT") of five modes, grouped by road types ("DESC_TERM") in the shapefiles. These functions work when calculating 'manually':

mean( "AADTallc_1" ,  "DESC_TERM" , "AADTallc_1" != 0.0)
mean( "AADTalllgv" ,  "DESC_TERM" , "AADTalllgv" != 0.0)
mean( "AADTTaxi" ,  "DESC_TERM" , "AADTTaxi" != 0.0)
mean( "AADTMotorc" ,  "DESC_TERM" , "AADTMotorc" != 0.0)
mean( "AADTLtBus" ,  "DESC_TERM" , "AADTLtBus" != 0.0)

I've created this flow on Graphical Modeler:

Where for LGVs (where the model fails);

With the failure message "An error occurred while evaluating the calculation string: Cannot use aggregate function in this context":

What could be going wrong?

I tried changing the parameters, but nothing seems to work.


texts: Writing ArcPy UpdateCursor based on row test?

texts:
I'm working on building a series of updates to street feature class based on fields inside the class. I was attempting to do this with update cursor to get better performance.

When the Streets_GC_FullName is null, the STREETS_NW_SLIPRD = 2 and STREETS_NW_RAMP=1, then change the Streetname field to "Exit Ramp"

The Cursor runs for a little while but then none of the rows are actually updated.
I had tried adding the print statment to see if it was stopping on those rows that matched, but none printed.

So I'm assuming that my row test is failing?

fc = r'C:\GIS\Sept2018_MergeLA.gdb\NW_GC'
fields = ['Streets_GC_FULLNAME', 'Streets_NW_FEATTYP', 'Streets_NW_SLIPRD', 'Streets_NW_RAMP', 'StreetName','Streets_NW_ID']

with arcpy.da.UpdateCursor(fc, fields) as cursor:
    for row in cursor:
        if (row[0] is None and row[3] == 2 and row[4] == 1):
            row[5] == "Exit Ramp"
            # print("Updated"+row[6])
            cursor.updateRow(row)
print ("Processing complete")


texts: Symbology sample box in QGIS is not showing required symbol?

texts:
The sample box that is used in the legend doesn't look like the actual map when I use point pattern fill in QGIS. The box is mostly empty with a partial point showing in the corner.

I have tried different distances/displacements but none of them move the point in the sample box.


texts: ArcGIS Online/ArcGIS Pro: Turn off symbol by attributes for each attribute

texts:
When I am working in QGIS I have the option to catogorize the symbology by attribute and that I can turn off the symbols in the layerspanel on and off easily like in the image below.

However, when I use ArcGIS Pro or ArcGIS online I want to do this too. But in the contents pane there is no option for this. Is there an easy fix for this?


texts:
wet<-readOGR(dsn=".",layer="shp") #wetlands
r.raster <- raster()
extent(r.raster) <- extent(wet)
res(r.raster) <- 100
wet.r <- rasterize(wet,r.raster) #rasterize
plot(wet.r)
b500<-buffer(wet.r,width=500) #buffer of the wetlands 500m
b1000<-buffer(wet.r,width=1000) #buffer of the wetlands 1000m
sum<-b1000 | b500 #aggregate them with "OR"
plot(sum)  #plot them

I have to change the values of the cells of these rasters, b500 and b1000. How to do that? Furthermore, how do i check the current values of a raster?


texts: Missing the results from grass tools when used in python console

texts:
I am running a grass tool from the python console to obtain a buffer from a vector polygon layer. But the result shapefile is not being created.

Running the same tool manually from the Processing tools as a GUI works just fine.

I am using Windows-QGIS 2.14.11-GRASS 7.2.0.

The vector layer is in variable part1 and the result should be going to a temporary folder (output is specified as None to the script). The call to the tool is:

result = processing.runalg('grass:v.buffer.distance', part1,'10','0.01',True,True, extent, -1.0, 0.0001, 0, None)

The command runs without problems and the variable result contains the name of the processing, the name to the result file:

{'result': u'C:\\Users\\username\\AppData\\Local\\Temp\\processing1dbc8a694a98406aa1f7656068b6f239\\985c82bb1ac3468f9a8d49b49ad53def\\output.shp'}

But there are no files there.

The same problem happes at least with grass:v.clean.

Anyone else has this problem?

Since the tool works when run manually from Processing as a GUI tool, is there a bug somewhere where in the python scripts?


texts: Java or .net platform is good for creating interactive web page

texts:
I have just now completed my degree in Spatial Information Technology, and learned how to use the ARCGIS platform and QGIS. Now I wanted to create a server/client webpage using mapserver for a project work. In the future, I am planing to go more into the GIS developer side. Please advise on whether you think java or .net platform will be better for creating interactive web pages.

I would like to use opensource tools like mapserver rather than ArcGIS Server due to cost.

Thank you


texts: Labels with halos become faded and transparent in CMV when using an aerial image

texts:
I am using CMV to create a web map but when I use esri's world imagery or another image service as my base map, my labels (only with halos - Aerial 10pt with 2pt halo) become faded and transparent -  see image below.

I have change the font style, size, colors, halos (note that the labels are fine without halos).


texts: Fast, simple assignment of flow directions to flats in GRASS

texts:
I have some enormously huge DEMs which I'd like to assign flow directions to. These DEMs are guaranteed to have no sinks, just flat areas.

It seems as though the only ways of assign flow directions to flats in GRASS involve doing depression filling (e.g. r.terraflow and r.watershed and r.fill.dir). In some cases, the fill is implicit, as in r.terraflow.

I'm looking for a simple, O(n) assignment of flow directions, similar to that presented by Barnes et al. (2014). Does GRASS include such a functionality?


texts: Compare and modify changes amongst GIS and MIS data

texts:
import arcpy
import sqlite3

# the column names
columnNames = ["ASSETNUM", "DESCRIPTIO", "INSERVICE", "ENGINEID"]

data = []
with arcpy.da.SearchCursor(r'E:\mypathtoshp\Fire_Hydrants.shp', columnNames) as sCursor:
    for row in sCursor:
        data.append(row)

# create a connection to your sqlite database
conn = sqlite3.connect(r'E:\mypathtodb\MyDB.db')
cursor = conn.cursor()

# create a table within the database
cursor.execute(''' CREATE TABLE Hydrants_MIS(ASSETNUM text, DESCRIPTIO text, INSERVICE text, ENGINEID text)''')

# insert your new records
c = ['?' for c in range(len(columnNames))]
c = ','.join(c)
sql = 'INSERT INTO {} VALUES ({})'.format('Hydrants_MIS', c)
cursor.executemany(sql, data)
conn.commit()

I have the above code to create a table in SQLite and populate the data from shp file. The SQLite table serves as a MIS repository and the shp serves as GIS repo. My client will induce changes in GIS data.

I want to run a python code every 1 hour to compare the 'ASSETNUM' values from the MIS data(SQLITE) and GIS data(shp), checking whether any data has been modified/added/deleted. If there is any modifications, perform the following operations:

If there are any new rows in GIS data, pick the assetnum and update the fields in MIS data. Eg: There is a new asset added by the client in GIS data. I need the script to update this new data in the MIS table.
If some rows in GIS data were deleted, update the fields in MIS data. Eg: Some assets were not in use and were removed by the client in GIS data. I need the script to update this in the MIS table.
If there are modifications in the data, the updation should happen at MIS table as well.

So far, I have fetched the 'ASSETNUM' field values from both the sources to compare. I am stuck at structuring the conditional statements.

    import arcpy
     import sqlite3

    conn = sqlite3.connect(r'E:\mydbpath\MyDB.db')
    cursor = conn.cursor()
    MIS_result = cur.execute("SELECT ASSETNUM FROM Hydrants_MIS")
    rows = MIS_result.fetchone("ASSETNUM")

    with arcpy.da.SearchCursor(r'E:\myshppath\Fire_Hydrants.shp', ["ASSETNUM"]) as sCursor:
    for row in sCursor:
        if <there are new assets >
            add the same records in sqlite table
        elif <some assets are deleted>
            delete the same records from sqlite table
         elif <change in field data of some assets>
            update the same changes in sqlite table
        else
            no change


texts: What is wrong with the time format used in TimeManager

texts:
I have used the TimeManager plugin in Qgis. Trying to set the time layer to animate features

However I have something of the following

I believed I have fulfilled format YYYY-MM-DD HH:MM:SS.ssssss with 2020-05-10 03:56:02.859499931.  So anybody can enlighten me what is wrong? Hope is not because of wrong typing

The time (originally in sec) is converted by python into datetime format and saved as a csv(delimited text file) and then loaded into qgis

A sample is reproduced here

Thanks


texts: WFS Filtered : The filters that worked for me in 1.8 - gives no results in 2.0

texts:
Using Quantum GIS 1.8.0 and 2.0.1 desktop on windows XP 32bit

Server: http://vmgeoserver.vd.dk/geoserver/VD/ows?

Setting filter for VD:vma_admdata:

"BESTYRER"  = '230' AND  "KODE_ADMIN_VINTERVEJSERVI"  = '230_V1' AND  "KODE_ADMIN_VEDLIGEHOLD_NI" = '230_N'

In 1.8 this returns the results I'm looking for.
Metadata for the resulting layer states Source for this layer:

http://vmgeoserver.vd.dk/geoserver/VD/ows?version=1.1.0&SERVICE=WFS&VERSION=1.0.0&REQUEST=GetFeature&TYPENAME=VD:vma_admdata&SRSNAME=EPSG:25832&FILTER= BESTYRER 230 KODE_ADMIN_VINTERVEJSERVI 230_V1 KODE_ADMIN_VEDLIGEHOLD_NI 230_N

Doing the same in Quantum GIS 2.0.1 gives no error messages but the resulting layer is empty. Checking metadata the Source for this layer is:

http://vmgeoserver.vd.dk/geoserver/VD/ows?version=1.1.0&SERVICE=WFS&VERSION=1.0.0&REQUEST=GetFeature&TYPENAME=VD:vma_admdata&SRSNAME=EPSG:25832&FILTER= BESTYRER 230 KODE_ADMIN_VINTERVEJSERVI 230_V1 KODE_ADMIN_VEDLIGEHOLD_NI 230_N

i.e.. identical.

What am I doing wrong?


texts:
I am trying to do some image analysis in Python but am unable to install GDAL due to limited rights. I have attempted copying some of the files and folders into my python folder but am constantly encountering the same error

ImportError: No module named _gdal

The searching that I have done always seems to point to some sort of installer for GDAL, but I am only able to copy the needed .py folders into my python folder.

Is there a way to do this on Windows and if not is there an alternative image reader that can properly open rasters?

sys.path results:

['C:\\Users\\ûser\\Documents\\Python Image Extractor\\Image Interpreter\\src', 'C:\\Python26\\ArcGIS10.0\\lib\\site-packages\\xlutils-1.5.2-py2.6.egg', 'C:\\Python26\\ArcGIS10.0\\lib\\site-packages\\xlrd-0.8.0-py2.6.egg', 'C:\\Python26\\ArcGIS10.0\\DLLs', 'C:\\Python26\\ArcGIS10.0\\lib', 'C:\\Python26\\ArcGIS10.0\\lib\\plat-win', 'C:\\Python26\\ArcGIS10.0\\lib\\lib-tk', 'C:\\Python26\\ArcGIS10.0', 'C:\\Python26\\ArcGIS10.0\\lib\\site-packages', 'C:\\Python26\\ArcGIS10.0\\lib\\site-packages\\PIL', 'C:\\Python26\\ArcGIS10.0\\lib\\site-packages\\win32', 'C:\\Python26\\ArcGIS10.0\\lib\\site-packages\\win32\\lib', 'C:\\Python26\\ArcGIS10.0\\lib\\site-packages\\Pythonwin', 'C:\\windows\\system32\\python26.zip', 'C:\\Program Files (x86)\\ArcGIS\\Desktop10.0\\bin', 'C:\\Program Files (x86)\\ArcGIS\\Desktop10.0\\arcpy', 'C:\\Program Files (x86)\\ArcGIS\\Desktop10.0\\ArcToolbox\\Scripts']


texts: Looking for best practices for determining turn penalties

texts:
Shortest path routing will often return routes with a lot of turns (zigzag shape).

If you want to get routes with fewer turns, you can use turn penalties to make turning more "expensive".

What are best practices for determining such turn penalties?

The use case:
I would specifically like to get fuel-efficient routes for vehicles. Slowing down, turning and accelerating costs more fuel than driving straight ahead. (Traffic signs - such as stop signs - don't have to be considered.)


texts: Issue editing nodes on a PostGIS table via QGIS

texts:
I launch a PostGIS table via QGIS DB Manager tool [PostgreSQL 9.4.17 - PostGIS 2.1.8 - QGIS 2.18.19]. The table geometry is MultiPolygon. I can toggle both Toggle editing and Node tool buttons. When I click on map canvas in order to move nodes, I get the following message:

And the Log Messages Panel shows:

I can update attributes values and move features.

The issue seems to be on editing nodes only.

--

After having read several posts about editing a PostGIS table from QGIS, I've tried the followings without success:

Connection parameters to DB using postgres default user or other user
PostGIS table has both a primary key and an index on the geometry column
Function ST_IsValid() returns only TRUE values (no invalid geometries)
I can't connect directly via the "Add PostGIS layers" button's dialog for the DB is too big (over 50 schemas) and displaying the content is very long... Thus I use DB Manager instead.
Settings > Options > Digitizing: tried different tolerance values both in map units and pixels. Also tried to set "To vertex" or "To vertex and segments".
Settings > Snapping options: same tests.

What could I be missing?


texts: Filter edges in a geometric network trace?

texts:
I need to trace a geometric network of sewer pipes, but I'm only interested in storm pipes. More specifically storm pipes that connect directly to a starting flag, and not those downstream/upstream of another pipe-type. Each starting flag might have several pipes connected to it of different types (e.g. sanitary, storm, CSO, etc.). Unfortunately, I can't query out the storm pipes because traces only accept geometric networks and not feature layers.

I've been reading about "in_edge_along_digitized_weight_filter" (in trace geometric network tool), but I can't figure out how to use it. The drop-down box is empty for this option. My geometric network was not built with weights, so I'm guessing that's probably why.

I would like to avoid rebuilding the geometric network if possible.


texts: Snap point features based on attributes of another feature class - ArcPy

texts:
I am trying to snap a point feature class to another point feature class based on matching attributes.
I have a script tool built and it works for a sample area where all records match, however when trying to run the entire dataset it will fail when it hits a row that cannot be found/matched.

I'd like this script to iterate through each row in pts_Survey and see if there is a match to the attributes in pts_MDP.

If there is a match, update the geometry of pts_Survey to snap to the geometry of pts_MDP.

How can I change the script to account for this?
Can I add an IF statement or should I consider re-working the code entirely?

import arcpy
from arcpy import env
from arcpy.sa import *

# Check out the spatial extension
arcpy.CheckOutExtension("Spatial")

pts_MDP = arcpy.GetParameterAsText(0)
id_MDP = arcpy.GetParameterAsText(1)
pts_Survey = arcpy.GetParameterAsText(2)
id_Survey = arcpy.GetParameterAsText(3)

pts_Defined = {str(row[0]):row[1] for row in arcpy.da.SearchCursor(pts_MDP,[id_MDP,"SHAPE@"])}
with arcpy.da.UpdateCursor(pts_Survey,[id_Survey,"SHAPE@"]) as cursor:
    for row in cursor: # loop through points
        row[1] = pts_Defined[str(row[0])] # move point geometry to match point geometry from dictionary

    cursor.updateRow(row) # update the geometry


texts: gdal_calc - FutureWarning: comparison to `None` will result in an elementwise object comparison in the future

texts:
I use gdal_calc in the loop:

 for inp, out in zip(input, output):

    input_str = ''.join(inp)
    print input_str
    output_str = ''.join(out)
    print output_str

    gdal_calc = 'python gdal_calc.py --format GTiff -A '+input_str+' --outfile='+output_str+' --calc="10*log10(A) " --NoDataValue=0 --type=Int16'
    os.system(gdal_calc)

Warning appears

FutureWarning: comparison toNonewill result in an elementwise object comparison in the future.
  myNDVs=1*numpy.logical_or(myNDVs==1, myval==myNDV[i])

Does anyone know what the problem is?


texts: Extracting country's administrative hierarchy from OpenStreetMap

texts:
I would like to extract all the cities along with their administrative hierarchy within a country.

I was thinking to use a mix of tag and geom information since tags are often inconsistent/incomplete:

Maybe start searching with the boundary=administrative tag which should exist for major subdivision, then list the node[place=city] within the area.

That would take quite a while though, especially when getting to smaller admin divisions.

Is there a straight way to do this?


texts:
Which unit in result of calculating the polygon's geometry programmatically using .Area property of IArea object?

For example, there are 2 polygons. Let's calculate geometry manually (in some projected coordinate system) - two areas are:

3.918650 [sq km]  and 8.984022 [sq km].

Next, I calculate the polygons area programmatically this way

        IFeature ftr1 = fc.GetFeature(0);
        IGeometry g1 = ftr1.Shape as IGeometry;
        IArea ar = ftr1.Shape as IArea;
        double strar = ar.Area;

The result is:

0.000543 and 0.001243.

The areas are not proportional in this 2 ways (manual and programmatical):

8.984022 / 3.918650 = 2.292632

0.001243 / 0.000543 = 2.289008

It's very close, but not equal.

using C#, arcgis 10.2.2


texts: Accounting for the heading of my position when calculating a new lat/long

texts:
Suppose that I am in a position where I know my directional heading, latitude and longitude. I throw a ball, and know that it lands X meters in front of me and Y meters to my left.

From here, I want to calculate the ball's landing site in latitude/longitude. I know I can calculate a vector of where the ball lands relative to me, but I'm unsure of how to also account for the directional heading of my position.

I'm not a GIS person by any means, and am trying to write a Python script that does something similar.


texts: Trigger on SQL Server 2008 ArcSDE 10 and ArcGIS Server 10

texts:
I need to build a sql server db trigger to update a date field with the current date when a record is inserted or updated, I think.  I will be using the date field as a way to recognize updates/inserts programmatically.  I went to sql server mngt studio and navigated to the database in the tree view >> opened up the programmability folder >> opened up the Trigger folder and right clicked >> Create new trigger >> I then copied the below sql into the window and hit F5.  It messaged back to me that the command was successful but when I tried to open the Triggers folder I didn't see it there.  When I tested in arcmap inserting a record I did not get the date updated to my date field.

What am I doing wrong and are there special considerations for creating Insert/Update/Delete triggers against SDE?  Thanks...the sql is below.

>CREATE TRIGGER t_I_TestTable
>ON dbo.TestTable
>FOR INSERT
>AS
>
>UPDATE TestTable SET TestTable.LastUpdated = Getdate()


texts: Spatial join based on by location and by attribute

texts:
I have tried to solve my problem using previous similar issues but I didn't find anything.

I want to get direction for a shp (shapefile) with routes (polylines) and then to join this information in shp with points (for each route correspond two points and each route has a unique id same with the points that lay on this route).

I use the shp with routes in order to get direction (steps: first I split the routes in each segment and them I use the COGO tool in order to get the direction by ArcGIS Desktop).

How can I join the information from the direction to each points?

The problem is that while I have the unique id from each route when I split them then I have the same unique id in many segments so it makes difficult to join each point with a certain segment from the routes.

I cannot use direct the spatial join tool since it can join e.g. a segment with two points (when the points from two or more routes touch each other)

E.g. in the picture below you can see two routes (lines) and three points. the blue lines is associate with the two blue points and the yellow line with the point. I want to join the information from the blue routes to blue points.


texts: Copy and paste certain cell characteristics of +500 cells in Excel at same time to create polygon shapefile?


texts: Show attached images inside popup with ArcGIS Online?

texts:
In ArcGIS Online popups how does one display attached images inline as an image instead of just a link? In case of multiple images I'm okay with just first one showing inline.

Under "Configure Popup" checking " Show feature attachments as links" results in text-only links.  If the "Pop-up Media" control is used to add an image, it shows inline, however one has to know the URL ahead of time. The drop down selector does not show attachments as a possibility.

Desired result

Popup showing 1st image inline, hidden beneath fold are rest of attachments as links.


texts: Creating multiple maps from a single layer - ArcGIS 10.1

texts: I am looking for an add in for ArcGIS 10.1 to create multiple maps from a single layer that consists of a number of polygons. I have come across such an add in for MapInfo. Is there such an add in available for ArcGIS, or is there a way do this within ArcGIS?


texts:
Can anybody tell me what wrong in this code. Here I am passing point geometry to query task :

JS code :

var districtName = $('#ddlDistrictd').val();
function showresults() {
    $.ajax({
        type: "POST",
        url: "Agriculture.asmx/GetDistrictpoints",
        data: "{'districtname':'" + districtname + "'}",
        contentType: "application/json; charset=utf-8",
        dataType: "json",
        success: OnSuccess,
        //error: OnError

    });

    function OnSuccess(data, status) {
        var myObject = JSON.parse(data.d);
        var querytext = JSON.stringify(myObject);
        var querytext2 = "shape  in (" + querytext + ")";

        map.graphics.clear();

        var graphic;
        var queryTask = new esri.tasks.QueryTask('http://[my_server]/arcgisserver/rest/services/CRD/CRD2/MapServer/0');
        var query = new esri.tasks.Query();
        symbol = new esri.symbol.SimpleMarkerSymbol();
        symbol.setStyle(esri.symbol.SimpleMarkerSymbol.STYLE_SQUARE);
        symbol.setSize(10);
        symbol.setColor(new dojo.Color([255, 255, 0, 0.5]));
        query.returnGeometry = true;
        query.outFields = ["*"];
        // query.InSpatialReference = { "wkid": 102100 };
        query.where = querytext2;
        queryTask.execute(query, function(featureSet) {

            if (featureSet.features.length == 0) {
                alert("No Results Found");
                // esri.hide(dojo.byId("loadingImg"));
            } else {
                var extent = esri.graphicsExtent(features);
                if (!extent && features.length == 1) {
                    // esri.getExtent returns null for a single point, so we'll build the extent by hand
                    var point = features[0];
                    extent = new esri.geometry.Extent(point.x - 1, point.y - 1, point.x + 1, point.y + 1, point.SpatialReference);
                }

                if (extent) {
                    // assumes the esri map object is stored in the globally-scoped variable 'map'
                    map.setExtent(extent)
                }

            }
        });

    }
}

I am getting the below error :

and response error is


texts: ArcGIS layers not scaling when attribute table is open

texts:
I've recently discovered an issue in Chrome and I'm wondering if anyone else has had this issue or if anyone can tell me how to fix it.  Basically, when I zoom in and out on my map when the AttributeTable widget is closed, everything scales correctly as shown in the video below.

Correctly Scaling

However, when I open the AttributeTable and try to do the exact same thing, nothing scales correctly and begins to offset.

Incorrectly Scaling

Has anyone else encountered this issue?  If so, what did you do to fix it?


texts:
I can't figure out why the EasyButton does't work. The following is the result of many modifications but none of them work as expected.

    Simple Leaflet Map

<link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/leaflet-easybutton@2/src/easy-button.css">


texts:
I have successfully setup a tile server as described here: https://switch2osm.org/serving-tiles/manually-building-a-tile-server-14-04/

What's troubling is the lack of clarify surrounding OSM Bright and default styles.

I've followed the documentation as described, I've managed to install TileMill and MapBox Studio, but I cannot seem to find a precompiled stylesheet that somewhat resembles the default look of openstreetmaps.org.

Following the directions for OSM Bright, I get a map of continents and oceans.  Does a precompiled XML file that works with mapnik exist which handles more than continents and oceans?

I had to delete many  xml nodes as they were throwing errors when attempting to run renderd and I'm assuming these XML nodes described the features I'd like to see styled.

Can anyone elaborate on how these components fit together and what I am doing wrong/not seeing?  Please point me toward some resouce(s) that further explains how I can compile a default streetmap with basic styles for a North American dataset.


texts: Why does QGIS v.clean toolset create result many times more than the original layers?

texts:
I am trying to correct overlap topology error type of a polygon layer. The red line marks the places where there were overlapping errors on this particular layer.

I tried to use the v,clean toolset from the QGIS toolbox to correct the errors from the layer. The tool corrects the errors from the layer but on opening the attribute table i noticed that the result was almost 3 times the number of the original number of polygons. Pls can anyone advice on how to use the tool or any other approach in correcting overlap topology error of polygon layers in QGIS.


texts: How much sidelap is adequate for aerial imagery? My AOI has few buildings taller than 3 stories. I want to minimize building lean. I've seen imagery with 40% sidelap and the building lean is pretty bad along the edges. Is 60% side and endlap overkill?


texts: '.Simplify' causes 'null geometry' errors (Python, OGR)

texts:
Question relates to an earlier one: Can't retrieve geometry (Python, OGR, Rtree)

Im trying to simplify a dataset with 180k+ polygons and place the simplified geometries in a separate dataset. With the use of Check Geometry tool in ArcMap the original dataset contains one error:

at feat_id = -1 PROBLEM: could not find spatial index

But I can't select this feature and, for example, remove it. When I continue and apply .Simplify() to each feature iteratively in my script, it introduces 'null geometry' errors. Is my use of .Simplify() as shown below correct or not? I tried both .Simplify(10) and .Simplify(5). The former results in more errors (71) compared to latter (8).

if.os.path.exists('class.shp'):
    driver.DeleteDataSource('class.shp')
outFile = driver.CreateDataSource('class.shp')
outLayer = outFile.CreateLayer('class.shp', srs, geom_type = inLayer.GetLayerDefn().GetGeomType())

selectLayerDef = selectData.GetLayerDefn()
for i in range(0, selectLayerDef.GetFieldCount()):
outLayer.CreateField(selectLayerDef.GetFieldDefn(i))

outLayerDef = outLayer.GetLayerDefn()
for i in range(0, selectData.GetFeatureCount()):
    inFeature = selectData.GetNextFeature()
    outFeature = ogr.Feature(outLayerDef)
    for i in range(0, outLayerDef.GetFieldCount()):
        outFeature.SetField(outLayerDef.GetFieldDefn(i).GetNameRef(),
                            inFeature.GetField(i)
    geom = inFeature.GetGeometryRef()
    geomSimple = geom.Simplify(10)
    outFeature.SetGeometry(geomSimple.Clone())
    outFeature.SetField('FEAT_TYPE', 'LAV')
    outLayer.CreateFeature(outFeature)
    inFeature.Destroy()
    outFeature.Destroy()

UPDATE:

I changed the above code to add a check on if the simplified geometry has issues. When this is the case it would clone the original geometry instead of the simplified geometry. It's certainly not pretty but it seems to work.

xmin, xmax, ymin, ymax = 0
outLayerDef ´outLayer.GetLayerDefn()
for i in range (0, selectData.GetFeatureCount()):
    inFeature = selectData.GetNextFeature()
    outFeature = ogr.Feature(outLayerDef)
    for i in range(0, outLayerDef.GetFieldCount()):
        outFeature.SetField(outLayerDef.GetFieldDefn(i).GetNameRef(),
                            inFeature.GetField(i))
    testFeature = ogr.Feature(outLayerDef)
    for i in range(0, outLayerDef.GetFieldCount()):
        outFeature.SetField(outLayerDef.GetFieldDefn(i).GetNameRef(),
                            inFeature.GetField(i))
    geom = inFeature.GetGeometryRef()
    geomSimple = geom.Simplify(10)
    testFeature.SetGeometry(geomSimple.Clone())
    testGeom = testFeature.GetGeometryRef()
        if testGeom.GetEnvelope() == (xmin, xmax, ymin, ymax):
            outFeature.SetGeometry(inFeature.GetGeometryRef().Clone())
            outFeature.SetField('FEAT_TYPE', 'LAVBEBYG')
            outLayer.CreateFeature(outFeature)
            testFeature.Destroy()
            inFeature.Destroy()
            outFeature.Destroy()
        else:
            outFeature.SetGeometry(geomSimple.Clone())
            outFeature.SetField('FEAT_TYPE', 'LAVBEBYG')
            outLayer.CreateFeature(outFeature)
            testFeature.Destroy()
            inFeature.Destroy()
            outFeature.Destroy()


texts: Converting from Boolean to point features using ArcPy?

texts:
I used geometry objects and da.searchCursor to check if randomly created points fall within polygon. Now I get True and False of each points that falls within or outside of the polygon. I made a list of them. Now I am facing problem in filtering  only the true ones and updating the point features class. I finally need a point point featureclass in geodatabse(I have already created personal geodatabse "mygdb"); that contains all the points that fall within the polygon. Here is the code where, I create a "point" feature class in gdb.  I need to update this featureclass with all the true coordinate pair (or points). I have very little knowledge to work with geometry objects.

fc = "C:/Users/pp/Desktop/proj/poly.shp"
arcpy.FeatureClassToGeodatabase_conversion(fc,mygdb)
fc1 = arcpy.CreateFeatureclass_management(mygdb, "point", "POINT")
polyGeom = arcpy.CopyFeatures_management("poly",arcpy.Geometry())

plist = []
with arcpy.da.SearchCursor(fc1, "SHAPE@") as cursor:
    for row in cursor:
        #print row[0].within(polyGeom[0])
        pp =  row[0].within(polyGeom[0])
        plist.append(pp)
    print (plist)enter code here


texts:
I have two raster layers. They come with the same cell size (100x100m), the same extent and are projected in the same CRS (+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs).
My problem is, that they are not exactly congruent, they are off by ca. 50m.

How can I make them congruent using ArcGIS 10.3 on Windows 10?


texts: Using Multiprocessing inside a plugin in Qgis

texts:
I spent quite a bit of time looking on how to use the multiprocessing package, but couldn't find anything on how to use it inside a plugin in QGIS.
I am developing a plugin that does some optimization for several elements. I would like to parallelize it.

I found a useful link on multi-threading inside a python plugin (http://snorf.net/blog/2013/12/07/multithreading-in-qgis-python-plugins/), but nothing on using the multiprocessing module, which might be easier?

I have been trying with a very basic example. I am only showing the run function from the plugin here:

def run(self):
    """Run method that performs all the real work"""
    # show the dialog
    self.dlg.show()

    # Run the dialog event loop
    result = self.dlg.exec_()

    # See if OK was pressed and run code
    if result:
        #Get number of cores
        nProcs = mp.cpu_count()

        #Start a Process
        p = mp.Pool(nProcs)

        #Define function
            def cube(x):
                return x**3

        #Run parallel
        results = p.map(cube, range(1,7))

When I run this code from the plugin in QGIS, it opens several QGIS windows, which then return errors (can't load layers, etc.).

What am I missing?

Do I need to start a worker first on another thread and then use multiprocessing there?

Or would we use another function from multiprocessing?

I am working under windows 7, using QGIS 2.10.


texts: Calculating distance between points in two tables and including name of second point in infowindow in Carto?

texts:
I am plotting community gardens in CartoDB and want to show the nearest subway entrance in the infowindow. I have two tables with the location of the gardens and the subway stations. I am able to get the closest subway station to each community garden by using the following query:

SELECT   nyc_greenthumb_community_gardens_fusion.cartodb_id,
MIN(ST_Distance(nyc_subway_entrance_export.the_geom, nyc_greenthumb_community_gardens_fusion.the_geom)) as MinDistance
FROM nyc_subway_entrance_export, nyc_greenthumb_community_gardens_fusion
GROUP BY nyc_greenthumb_community_gardens_fusion.cartodb_id

How do I also return the name of the subway station?

It is in the nyc_subway_entrance_export table as nyc_subway_entrance_export.name.

I am assuming that once I have the name of the subway entrance, I can do a join to add the name and distance as new columns in the Community Garden table.


texts: Error when clipping raster in QGIS: "'ascii' codec can't encode character"?

texts:
I have done it many times (clipping a raster with a polygon layer or by extension) but today I got an error for the first time and I don't know how to solve it. I am always using the same type of raster layers. When running the clipping tool I got the following message:

gdal_translate -of GTiff -ot Byte -projwin 410282.185971 4694950.41162 410783.882475 4694524.96815 -co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 -co TFW=YES "C:/Users/FERRAN ALA/Desktop/Justificació Fontanals/QGIS/Topo 1_5000/bt5mv20sd0f283077st1r031.sid" "C:/Users/FERRAN ALA/Desktop/Justificació Fontanals/QGIS/background/image.tif"

GDAL command output:

  'ascii' codec can't encode character u'\xf3' in position 201: ordinal not in range(128) See log for more details

Does anybody know what that is?


texts: GeoServer YSLD - Point Stacking with multiple features on map

texts:
So, there isn't much information out there: https://docs.geoserver.org/stable/en/user/styling/ysld/reference/transforms.html#point-stacker

And I have successfully used this doc to create a simply style that stacks soil sample points on my map (Geoserver for Boundless Suite).

What I have been given now is dataset with multiple types of features where I need to use the
filter: ${entity_feature_here = 'Example'} multiple tikes to have each of the (10+) features render on the map. When I do this, it obviously looks cluttered. What I can't figure out is how to use Point Stacker to cluster each of these points until they are zoomed in.  When I do try, all my data seems to vanish.

name: YSLD_Point_Staclk
feature-styles:
- name: PointStacker
  transform:
    name: vec:PointStacker
    params:
      cellSize: 64
  rules:
  - name: 'Example 1'
    filter: ${entity_feature_here = 'Example1'}
    title: Example 1
    scale: [min, max]
    symbolizers:
    - point:
        size: 10
        symbols:
        - mark:
            shape: circle
            fill-color: #some color
  - name: 'Example 2'
    filter: ${entity_feature_here = 'Example2'}
    title: Example 2
    scale: [min, max]
    symbolizers:
    - point:
        size: 10
        symbols:
        - mark:
            shape: circle
            fill-color: #some color
  - name: 'Example 3'
    filter: ${entity_feature_here = 'Example3}
    title: No Plant Area
    scale: [min, max]
    symbolizers:
    - point:
        size: 10
        symbols:
        - mark:
            shape: square
            fill-color: #some color

I have tried it with and with setting the [min,max] and I have tried getting rid of the scale all together as is shown in the example in the Geoserver doc. Neither has worked, my data just disappears on the map.


texts:
I have BGS 50k vector data which includes a field referring to the LEX code on the BGS website. Obviously the web address will be unique for every geological formation so I need to insert the field into the map tips window.

I have laid out the information in table format which is why I have
 and  tags. I have tried typing the following with no success.

BGS Link:

or

BGS Link:

and

BGS Link:	<a href="[% "LEX_WEB" %]">

or

BGS Link:	<a href=[% "LEX_WEB" %]>

The last two produce the hyperlink at the very top of the map tip but when clicking it, my web browser does not follow it. Strange thing, I can't even reproduce the hyperlink outside the table format!

If you are keen to test it out I have left the my table styling below. I propose setting up a shapefile with a 'Link' field and entering something like 'www.google.com'.


texts: Using Update Cursor within Search Cursor with ArcGIS Pro?

texts:
I am new to using ArcPy/Python and I'm running into problems with a script I'm trying to write.

I'm using IDLE Python 3.6.8
ArcGIS Pro 2.4.0

Here is what I'm trying to do - I have 2 feature classes: 1. Water Mains, 2. Multi Ring buffer.

Within the Water Mains attribute table there is a "CriticalityScore" field. This is an integer value that is assigned based on the diameter(inches) of the pipe. Small diameter pipes = low value, larger diameter pipes receive higher values.

I am tasked with writing a script that will calculate this score for newly added pipes. I created a coded domain for a field called "Participation" to indicated "Yes/No". New pipes that are added into the attribute table by default receive "No."

What I have done so far is to assign if diameter = x then criticality score will = y.

Here's how the scores might look

2" = (score of) 5
4" = 9
10" = 15 etc..

What I'm struggling with now is:

Only updating "CriticalityScore" if "Participation" = "No"

Here's what I have so far

#set your workspace, where is this pulling data from? (gdb file path)
  ws = env.workspace = r"H:\Python Scripts & Tools\Script Tester - Crit Scoring.gdb"

  #set the name of the shapefile path with which you're working (water main .shp)
  fclass = r"H:\Python Scripts & Tools\Script Tester - Crit Scoring.gdb\Water_Mains.shp"

  buffers = r"L:\Engineering\Brice\CriticalityScoringwPipeTestBrice.gdb\BufferParcelEmergancyResponceFac"

  #The third item in "fields" (index position [2])should be the Yes/No field in the feature class
  fields = ["Diameter", "CriticalityScore", "Participation"]

  with arcpy.da.UpdateCursor(fclass, fields) as cursor:
      while row[2] == No:
          for row in cursor:
            if row[0] == 2:
                row[1] = "5"
            elif row[0] == 3:
                row[1] = "7"
            elif row[0] == 4:
                row[1] = "9"
            elif row[0] == 6:
                row[1] = "11"
            elif row[0] == 8:
                row[1] = "13"
            elif row[0] == 10:
                row[1] = "15"
            elif row[0] == 12:
                row[1] = "17"
            elif row[0] == 16:
                row[1] = "20"
            else:
                row[1] = "22"
            cursor.updateRow(row)
          print("The script has finished.")

The script was working before I added the "While row[2] statement. As of right now it will run through the script and give me the ending print statement, but the .aprx/attribute table is never actually updated.


texts: ArcPy SearchCursor getValue leads to NoneType object has no attribute lstrip?

texts:
I am having a little bit of a dilemma. I am having a problem with using getValue() and striping the value that is returned. For instance, take the following code.

import arcpy as ap
mxd = ap.mapping.MapDocument("CURRENT")
mylayer = ap.mapping.ListLayers(mxd, "Structures")[0].name

    for row in ap.SearchCursor(mylayer):

        str =  row.getValue("Struct_ID")
        newstr = str1.lstrip('0')

    if str == "662":
        print(str.lstrip('0'))

The error I get no matter what is:

Traceback (most recent call last):
  File "", line 9, in
AttributeError: 'NoneType' object has no attribute 'lstrip'

The code works fine if you take out the newstr = str1.lstrip('0') line. The   print(str.lstrip('0')) statement works fine. I just cannot use str.lstrip anywhere else in the code other then inside of print.

Anyone know what I am doing wrong or is it something to do with Esri?


texts: How to disable addition of database name to layer name when adding a vector layer in QGIS?

texts:
I am using QGIS version 2.18.14.

When I add a vector layer from a file geodatabase, I choose the OpenFileGdb driver and then pick the layers I would like to have on the map.

However, the layers are added with the name following the format <%DB_NAME% %LAYER_NAME%>. This makes it impossible to use the DBManager SQL editor directly for the added layers as I need to rename the layers first. This is because the database table cannot contain a white space. So, I need to rename "NYC streets" into "streets" or "NYC_streets".

Is there any setting in QGIS that can let me disable prefixing the layers I add with the database name from which they come?


texts: Exhibit "Attribute Table" Field values when mouse hovers on top of a feature in Shapefile in QGIS2.8?

texts: I am drawing a line to represent a Fibre Cable. I would like to see the Attribute Table of this line pop up the field values assigned to this line representing a Fibre Cable (Aerial or Underground).


texts: How I can load tileLayer in Leaflet framework using local tiles?

texts:
I'm using file:///C:/path/tiles/{z}/{x}/{y}.png like url of tileLayer but it doesn't work and no error is shown.

var layer = new L.tileLayer("file:///C:/OSM/tiles/{z}/{x}/{y}.png", {
      minZoom: 1,
      maxZoom: 15
}).addTo(map);

Any idea?


texts: ArcGIS10.2 Field calculator problem with strings

texts:
I am trying to do a simple concatenation using the ArcGIS 10.2 field calculator with Python code. I have a field of census tracts numbers (tract). I want to concatenate this with the numeric codes for the state and county so that I can join this table to the spatial geometry. The new field is defined as 'text'.

My code is "06085" + "!tract!" This yields the correct concatenation PLUS .0 on each record. So 0608512345.0

I have also tried str(06085)+str(!tract!) which yields the same results. What am I missing?


texts: Selecting points within polygon using ArcGIS Desktop?

texts:
I am using ArcGIS Desktop for the first time in many years.

I have a postal sector layer under a series of merged drive time boundaries.

How do I select just those postal sectors which fit within the drive time boundaries?


texts: Uploading maps to QGIS cloud and control scale?

texts:
I have managed to upload my maps on QGIS but when viewing them on the web or mobile platforms I am not able to zoom out as much as I could on QGIS desktop

What could be the problem? what settings should I change so that they can be view-able and allow zooming in and out.


texts: Is it possible to create CAD file from shapefile and also create CAD blocks from the database fields

texts:
I am trying to find out if there is any software or plugin/code that I can bundle that will take a shapefile and export it to a CAD file with CAD blocks generated so our CAD mappers can just move the blocks and finalize the drawing.

For example, I have the following data in the shapefile:

lat
long
detail1
detail2
detail3
...

I want to have it output:

detail1
detail2
datail3


texts: ArcPy raster output spills over environment extent by one row and one column

texts:
I'm resampling a raster, and I would like it to line up and match the dimensions of another template raster. For whatever reason, the output of the Resample_management function has an added row and column.

import arcpy
from arcpy import env

# Set file paths
in_raster = "inputRas.tif"
templateRaster = "templateRas.tif"
out_raster = "outRas.tif"

# Set 'Snap raster' and 'Extent' environment variables
arcpy.env.extent = arcpy.sa.Raster(templateRaster).extent
arcpy.env.snapRaster = templateRaster

# Get cell size from template raster (which is 1 m)
cellsizeX = arcpy.Describe(templateRaster).children[0].meanCellWidth
cellsizeY = arcpy.Describe(templateRaster).children[0].meanCellHeight
cellsize = str(cellsizeX) + " " + str(cellsizeY)

# Resample input raster
arcpy.Resample_management(in_raster, out_raster, cellsize)

# Print environment extent
print arcpy.env.extent

# Print extent of the output raster
print arcpy.sa.Raster(out_raster).extent

The resulting extents:

582568 5816053 583282 5816790 NaN NaN NaN NaN
582567 5816053 583282 5816791 NaN NaN NaN NaN

The function has essentially added one row of cells to the top of the output and one column of cells to the left. Why?


texts:
I want to use 'Geometry by expression' in a model I'm working on but I couldn't figure out how to pass the layer name parameter. Basically I have two  polygon layers (A and B) and I want to do a feature-by-feature difference using a shared attribute. The expression I have below works outside the modeller.

Inside modeller, since the two layers are intermediate outputs, and not an input, I do not know what to use for the layer name parameter ('layer B').

difference( $geometry, geometry( get_feature( 'layerB', 'attributeB', attribute($currentfeature, 'attributeA'))))


texts: Find percentage overlap (by category) between census tract polygons and circular polygons in QGIS?

texts:
Context
This is my second attempt to post this question since the answers to my first question were not as useful as I had hoped. The link to my first question is here: first question attempt
I am trying to solve a problem at a national level but am starting with Rhode Island to solidify my approach. Thus, everything in this question will be focused on the state of Rhode Island.
Here is the map I am starting out with:

The polygons are census tracts and the points are company office locations.
I’d like to figure out the following: For each census tract in Rhode Island, what percent of land area is covered by a certain company (total of 3 companies). Making sure that water area is not included in the calculations.
And just as an FYI, I am extremely new to GIS in general and have been going at this problem for a few days now. I am using QGIS 2.18.9.
Inputs
I have 3 input files:

The 2010 Census Tract shapefile (RI_Tracts) for Rhode Island with edited fields for simplicity:
A CSV file (RI_Tract_Tiers) that classifies each census tract into tiers 1-5:
A CSV file (RI_Offices) for office locations with the following fields:

Desired Output
A table that looks something like this:
What I've Tried

Join the RI_Tracts table with the RI_Tract_Tiers table to classify the Tracts into tiers
Color code the census tracts by Tract_Tier
Import and buffer the company locations in the RI_Offices file

Buffer distance is 10,000ft (3048 meters), so I enter “0.030480” into the distance field of the Fixed Distance Buffer tool

Dissolve the company locations by Office_Category to get rid of any overlaps b/w the same company
Color code the company locations by Office_Category

Here is the map I end up with after performing the steps above:
This is where I get stuck. I have tried using the following tools to analyze the map and get the output I want (with no success):

Union
Intersect
Statistics by Categories


texts:
I am using ol-cesium and I want to receive an event when the map is rotated. So my approach is to use the change:rotation event of the ol.Map, but it is never called.

this.map = new ol.Map({
        layers: [
            tileLayer
        ],
        controls: [],
        target: 'map',
        view: new ol.View({
            center: ol.proj.fromLonLat([12.372542, 47.426919]),
            zoom: 12,
            minZoom: 12
        }),
    });

    this.map.getView().on('change:rotation', () => {
        console.log("fired");
    });

I am also using typescript but this shouldn't be the problem. My development browser is Google Chrome.


texts: HeatMap to identify hotspots on cached service from ArcGIS Server?

texts:
ArcGIS Server lets you create cached map services.
These can be consumed by various applications (ArcGIS Desktop, web apps, etc)
I would like to get an idea on which areas of my cache are getting hit the most.
I would like to visualise these results using a heat map overlaid onto the same cached map.
I was thinking of utilising something like the HeatMapsAPI and one of the APIs provided by ESRI (Javascript API should be sufficient)

I know the cached maps might be slightly changing at v10 of ArcGIS Server, but im trying to think of the best way of determining which tiles get requested the most without putting a significant load on either ArcGIS Server or the web server (IIS).

I was originally thinking looking at the actual cache directory under the arcgiscache folder, and perhaps hooking into the Date Accessed attribute (but have little idea on how to actually implement this).

I would like to hear from anyone who has experience with ArcGIS Server, on if they have any suggestions on how to do this?


texts: Route From Network Analyst Not Fully Displayed

texts:
I am working on a project for finding the optimized route solid waste pickup within a defined neighborhood.
These are my setting for route properties. The curb approach must always be on the right side of the vehicle.

The system gave me a route once I hit the solve button and if I offset the route and add arrows to it to show the flow of the direction. I got the following

How come it didn't show the missing route (I drawn in red)?


texts: GRID ERROR in raster calculator of ArcGIS Desktop?

texts:
I wrote and ran this formula in raster calculator in ARCGIS. But the following error gets me into trouble.

Does anyone have a suggestion to overcome this?

Here is the formula I wrote and the error that appears.

What is wrong with this formula?


texts: How to create a rule-based gradient using an expression

texts:
I am mapping points with values from 0-100 and want to display these using a gradient color scheme. However I also want to specify the breaks as regular intervals along this range. For a gradient QGIS automatically chooses breaks but doesn't allow me to specify the start, stop, and interval size. Alternatively I can use rule-based styling but it becomes tedious to update the point sizes, colors, etc.

Is there are smarter way to get QGIS to create a gradient style but retain the flexibility of the rule-based cutoffs as well? The points share the same attributes except for color, which is determined by gradient.


texts: Getting polygon country boundaries in decimal degrees?

texts:
Using Google Maps, I'd like to create draggable polygons (on the fly) of country boundaries (they don't need to be particularly precise).  Something similar to this, but for countries.  So I want the boundaries in decimal degrees, but I can only find the data in .shp or kml/kmz files via sites like:

shp:

http://www.naturalearthdata.com/
https://hiu.state.gov/data/data.aspx
http://www.gadm.org/

KML:

https://www.google.com/fusiontables/DataSource?docid=1N2LBk4JHwWpOY4d9fobIn27lfnZ5MDy-NoqqRpk

... but I'm unsure about how I go about extracting decimal coordinates from those data sources.

How do I solve this problem?


texts: How do I get the layer from a featureSource within a QGIS plugin?

texts:
I'm writing a QGIS processing toolbox plugin 3.14 and need to modify its features:
So am doing something like this:
source = self.parameterAsSource(parameters, 'LAYER', context)
features = source.getFeatures()
for current, feat in enumerate(features):
  # do something to each feature

But how do I get the layer from the source? To make it work, I'm doing something like this:
canvas = iface.mapCanvas()
layer = canvas.currentLayer()
layer.startEditing()
# Then inside the loop above
  layer.changeAttributeValue(feat.id(), 2, value)
# And then finally:
layer.commitChanges()
layer.reload()

All works fine but am confused about not being able to get to the layer from the source? Is this the only way?
Thanks for any help!