texts: Creating buffer zones in meters by reprojecting CRS but in new CRS layers do not align using QGIS

texts:
I am trying to create buffer zones in meters/kilometers around a series of sites in Montreal, Canada. I used Google maps to obtain the coordinates for each site, and saved them in a .csv file. After importing my .csv file using the open data source manager, and adding an .shp layer I downloaded from the government of Canada, I get an initial map that looks like this:

First inconsistency: I notice that even though my site layer CRS is projected in web mercator, the default from Google maps (EPSG:4326 WGS83), but my project and layer CRS is EPSG:3347 Statistics Lambert Canada, however both layers align just fine after I import my points file

When I try to create my buffer zones, I find I can only use degrees, and QGIS suggests reprojecting my points file into a local CRS.

I then change the points file CRS to match that of the project, and my sites end up somewhere in the Pacific ocean, despite the matching CRS (EPSG 3347)!

I tried the solution from a similar post at Delimited text layer doesn't align with shapefile, but it didn't do the trick since my points were already aligned with the .shp file after importing them.
What am I missing ?
How can I use a local coordinate system to allow buffer zones in meters/kilometers, while keeping both layers aligned ?


texts: How to prevent label rendering in bleed space in QGIS composer?

texts:
I'm making maps for a book where maps run to the edge of the page. The printer has advised a 3mm bleed space around images, therefore if my page is 127x177 mm, images should be 133x183 mm. This bleed space is to allow for errors when the book pages are cut, and in the case of maps should not have any essential information. So I use a composer page size that includes these extra 3 mm on all sides.

I don't want labels for towns to be cut off at the edge of the map, therefore I have disabled the "Show partials labels" option. However, labels are still being shown within this 3 mm "no-man's-land" that may be cut by the printer.

The solution I have considered is to make two overlapping, concentric and perfectly aligned maps, the lower map runs to the edge of the bleed and would not show the towns, the map on top would be the exact size of the book page (i.e. 3mm less on each side than the underlying map) and would show the town layer, therefore creating a 3 mm border where there are no towns (or labels) shown.

However, I can't make these two maps align perfectly in composer. I have tried calculating the extents for these two maps so they have the same center and the same scale, but QGIS composer seems to have a fixed 3 decimal place precision so I can never get them to align perfectly.

Another way to approach this problem would be a new feature which would allow not rendering labels within the atlas "Margin around feature".

Is there any way I can prevent labels from being rendered within a bleed border?


texts:
On Sql Server i have a table:

CREATE TABLE [dbo].[test](
[id] [int] IDENTITY(1,1) NOT NULL,
[geom] [geometry] NULL,
 CONSTRAINT [PK_test] PRIMARY KEY CLUSTERED
(
[id] ASC
)WITH (PAD_INDEX = OFF, STATISTICS_NORECOMPUTE = OFF, IGNORE_DUP_KEY = OFF, ALLOW_ROW_LOCKS = ON, ALLOW_PAGE_LOCKS = ON) ON [PRIMARY]
) ON [PRIMARY] TEXTIMAGE_ON [PRIMARY]

GO

With data:

DECLARE @g geometry;
SET @g = geometry::STGeomFromText('
CURVEPOLYGON(
  CIRCULARSTRING(379821 1582021, 379816 1582026, 379821 1582031, 379826 1582026, 379821 1582021)
  )
 ',0);

INSERT INTO dbo.test(geom) select @g

But when i publish table on geoservet in Layer Preview u got error:

Unknown WKB type 10

Geoserver cannot work with Curves? Any advises?


texts: Trying to get easting and northing points to be visible in QGIS

texts:
I'm a complete beginner on QGIS.
I have a load of Ordnance Survey NGRs which I have converted to eastings and northings of 5 digits each, imported them as a csv file into QGIS and tried to create them as a layer so that the individual site locations are visible.  The only other layer I have is a simple outline of the UK (OS high water polyline). The layer with the eastings and northings is at the top of the layers panel, but they are not showing up.  They are set to show up as simple coloured markers.
The project is set to OSGB 1936 (EPSG: 27700).
Any ideas where the location markers might be and how to make them visible?


texts: Importing data with empty cells without GIS replacing the empty cell with 0

texts:
I am importing data from excel into ArcGIS 10. The data are columns of GPS data, associated date and time, and a unique point identifier. For some of the rows the GPS data is missing and so the Lat and Long columns are empty but the other columns retain dates and times etc. When I import this data into ArcGIS 10 a 0 is placed in these cells and so Arc plots the points as 0,0.

Previously I have just removed the empty rows from the imported data table to get around this problem. However, I now need to import the full data table including the empty rows. How can I do this without Arc entering the empty cells as 0?


texts: MapBasic issue in creating Thematic with Shade

texts:
I'm trying to write a bit of code that looks the columns of a table and if they match by text, then generate a thematic of predefined ranges.
So far:

Include "MapBasic.def"

Declare Sub Main

Sub Main()
    Dim ColumnCount as integer
    Dim ColumnName as string
    Dim TableName as string
    Dim aColAlias as Alias
    Dim WinID as integer

    WinID = FrontWindow()

    TableName = tableinfo(1,TAB_INFO_NAME)

    'Assume table 1 for this test
    for ColumnCount=1 to tableinfo(1,TAB_INFO_NCOLS)

        ColumnName =  ColumnInfo(1,"COL" & ColumnCount,COL_INFO_NAME)

        if instr(1,ColumnName,"RSRQ") >1  and instr(1,ColumnName,"SS") > 1  then
            print TableName & "." & ColumnName 'works as expected
            aColAlias = ColumnName

            'Create Thematic
            Shade Window WinID TableName with aColAlias
            ignore 0
            Ranges
            -1 : -4 Symbol(108,6750207,10,"Wingdings",0,0),
            -4 : -7 Symbol(108,8454016,10,"Wingdings",0,0),
            -7 : -10 Symbol(108, YELLOW,10,"Wingdings",0,0),
            -10 : -13 Symbol(108,16757529,10,"Wingdings",0,0),
            -13 : -60 Symbol(108,RED,10,"Wingdings",0,0)

            exit For
        End If
    Next

End Sub

The Shade command always fails with 'Invalid variable Type'
I've tried replacing the Column Alias with the string of the column name, but that doesn't work either. I'm sure its something simple, but I just can't see it.


texts: How do you generate WKT from a graphic drawn on a map?

texts:
All,

In order to do some quick analysis using SQL spatial, I'd like to generate spatial objects using the STGeomFromText / ST_GeomFromText function.

With point data, I can do this by geocoding an address and using desktop GIS to give me the X/Y, which I then turn into a point using a syntax like this:

declare @xyPoint geometry;
set @xyPoint = geometry::STGeomFromText('Point(3175288.989  1711307.592)', 2232);

(This will create a temporary point geometry at the X,Y coordinates above in the 2232 (CO State Plane Central coordinate system)... I can then use the temporary point to drive distance calculations, etc.)

My question: is there a way to do something similar with a polygon: maybe draw a graphic in a desktop GIS map (ArcMap, QGIS, etc.), get the WKT from that graphic, which I can then use in a query to create a polygon object using the STGeomFromText function?

Basically, this is the opposite of the Quick Draw plugin for QGIS


texts: Missing checkbox when selecting layers using QGIS?

texts:
I'm an absolute beginner to QGIS.

I'm trying to select BRUNSWICK_1106_TreeTops.dbf as a layer, but it does not have a checkbox in front of it as the DRUHIG_0096.tif.

How can I get the checkbox in front of the BRUNSWICK.dbf so that I can view it as a layer?


texts: Splitting line if one or more polygon layers intersects at same time using QGIS?

texts:
I have a road network in the form of a line, that I need split up into sections, if a polygon intercepts.

This is fairly straight forward, but I have several layers of polygons where I need the full extent of the line, if it is continuous across several polygon layers.

I've tried to illustrate this with a picture, what I need is basically all the lines marked with Red on the picture, to be a separate new line.

As the picture shows, sometimes only a green polygon intersects, but other times it's a green polygon followed by a purple, and then I need it to be a full line across the several polygons. What is the easiest way to achieve this in QGIS?


texts: Python, OGR and Mapnik: Possibility of passing existing ogr.Open() datasource into Mapnik without using mapnik.Ogr()?

texts:
I'm trying to modify create_utfgrids.py to use an OGR VRT file, rather than shapefiles, and while the OGR stuff in this code responds to the VRT, the Mapnik stuff is not. It's failing with..

ERROR 1: Failed to open datasource `PG:host=127.0.0.1 user=postgres dbname=gis password=passw0rd tables=parcels_20140829_pmerc'.

I believe the VRT is fine, because it not only works with ogrinfo directly, but as stated above, the OGR-related code in the Python file also works.

While I'm pretty confident the VRT isn't the issue, for those who might ask, the VRT I tried looks like this. I know there's some redundant stuff in there, but that's because I thought maybe Mapnik was looking for something specific, so I was adding any node I could think of which might be relevant:

<OGRVRTDataSource>
  <OGRVRTLayer name="parcels_20140829_pmerc">
    <SrcDataSource>PG:host=127.0.0.1 user=postgres dbname=gis password=passw0rd tables=parcels_20140829_pmerc</SrcDataSource>
    <SrcLayer>parcels_20140829_pmerc</SrcLayer>
    <LayerSRS>EPSG:3857</LayerSRS>
    <SrcSQL>SELECT ST_AsBinary(wkb_geometry) as geomm, tms, owner_name, the_wkt FROM parcels_20140829_pmerc</SrcSQL>
    <GeometryField encoding="WKB" field="geomm"></GeometryField>
    <Field name="tms" src="tms" type="String"/>
    <Field name="owner_name" src="owner_name" type="String"/>
    <Field name="the_wkt" src="the_wkt" type="String"/>
  </OGRVRTLayer>
</OGRVRTDataSource>

I eventually got past the problem by passing in an OGR PostgreSQL connection string, parsing over that, and then pushing the values into a mapnik.PostGIS() constructor. Sure, this solves my immediate problem, but it's not quite as generic/universal as I would like.

So I'm wondering, since I can't get Mapnik to like my VRT, if it's possible to reference the existing ogr.Open() datasource object from Mapnik and avoid having to use something like mapnik.PostGIS() or mapnik.Ogr() in the first place.

To get this idea, this is kind of how I imagine it..

import mapnik
import ogr

# create_utfgrids.py uses shapefiles, hence "shppath".
# But for this exercise it represents a path to a VRT.
ds = ogr.Open(shppath)

# A bunch of OGR stuff here.

# This is where I want the magic to happen,
# pushing the OGR dataset into Mapnik..
ds = mapnik.Ogr(ds)

# And then Mapnik stuff follows.


texts:
I got a tool for ArcGIS which has its own license manager (lmtools / flexnet). Everything is working fine.

My problem is:
When i use ArcGIS and start the tool the license is in use. Ok. When i stop the tool the license is still in use until i restart ArcGIS.

Is there any way (with a script or workaround) to tell the licensemanager:
"Hey i finished working, somone else can use it now"

Its annoying to restart ArcGIS everytime :/


texts:
Is there any analog in ArcGIS Desktop (10.6) of Photoshop's magic wand tool? I have a raster layer that has clearly defined polygons is there any tool that is able to click to select the areas and create a polygon from them?

Also if this is not possible is it possible to write and add a tool that is able to actively be used or can you only implement scripts via creating custom geoprocessing tools?


texts: How to deal with flawed line geometries in Geopandas?

texts:
I have a big dataset containing line-geometries. I would like to delete overlapping lines using the intersection statement in geopandas overlay in the following manner:
import geopandas as gpd
from shapely.geometry import LineString

df1 = gpd.GeoSeries([LineString([(1,1), (0,2), (0,5), (0,6)]),LineString([(9,9), (10,12)])])
df2 = gpd.GeoSeries([LineString([(0,4), (0,7), (0,9), (0,12)]),LineString([(14,20), (13,13)])])

df1a = gpd.GeoDataFrame({'geometry': df1, 'df1':[1,2], 'Class':1})
df2a = gpd.GeoDataFrame({'geometry': df2, 'df2':[1,2], 'Class':2})

overlay = gpd.overlay(df2a,df1a, how='difference')

But it seems some lines do not intersect, maybe they are lying on a very small scale not on each other but parallel, indicated by checking with:
line1[['geometry']].intersects(line2[['geometry']]) --> False

How can I fix the geometry, so that all lines in my dataset are pefectly aligned (Make parallel overlap), that I can run the Set-Operation ? Preferably I would like to have a solution in python (Geopandas), but in worst case also QGIS would be ok.


texts: Learning variety and certain combinations of land use changes using ArcGIS Desktop

texts:
I have 3 Corine land cover rasters, which each contains types of land use, for simple example - 1 - arable land, 2 forest, 3 - urban areas etc.
How do I get not only the variety (which is easily obtained by cell stats) but also which type of land use changed to which, for all 3 of the cases.
So for example between years 2006 - 2012 from 1 to 2 and between 2012 - 2018 from 2 to 3.
The desired output could be a .dbf. I was thinking of using the "diff" tool, however I didn't come up with a reasonable and good looking output.
Any suggestions?


texts: "Item not found in this collection" when trying to set the Shape on IFeatureBuffer

texts:
I'm working on code to open a 3rd party data file and convert it into a collection of features in a custom layer based on FeatureLayer. I have had no issue importing the data itself (nodal data for each 'cell' as a feature, such as elevation).

Because the data files are on average rather large (10k+ cells), I'm using a cursor implementation with an IFeatureBuffer that is filled with each cell one at a time. However, when I try to set the geometry of the cell via IFeatureBuffer like so,

IGeometry cellGeometry = getCellGeometry(cell);

featureBuffer.Shape = cellGeometry;

I subsequently get the error "Item not found in this collection" called on the second line. My first thought was that I was creating the cell geometry incorrectly in some way, so I tested it by creating a simple point as the geometry. This failed with the same error. I next thought it might be an issue with ZAware. I tried setting it to both ((IZAware)cellGeometry).ZAware = true; and ((IZAware)cellGeometry).ZAware = false; with the same issue both times.

I'm thinking it might have to do with how the fields in the FeatureClass were setup, but I'm quite stuck at the moment.

UPDATE:

I've tried updating how the fields are created to fix any issues that may exist therein, and the issue remains. The fields are generated like so,

    {
        if (this.ImportedField == null)
        {
            return null;
        }

        IObjectClassDescription fcDesc;
        if (this.DatasetType == esriDatasetType.esriDTTable)
            fcDesc = new ObjectClassDescriptionClass();
        else
            fcDesc = new FeatureClassDescriptionClass();

        IFields fields = fcDesc.RequiredFields;
        IFieldsEdit fieldsEdit = (IFieldsEdit)fields;

        //Setup shape field
        IField shapeField = fields.Field[fields.FindField("Shape")];
        IFieldEdit shapeFieldEdit = (IFieldEdit)shapeField;
        shapeFieldEdit.Name_2 = "SHAPE";
        shapeFieldEdit.AliasName_2 = "SHAPE";
        shapeFieldEdit.Type_2 = esriFieldType.esriFieldTypeGeometry;
        shapeFieldEdit.Length_2 = 0;
        shapeFieldEdit.Precision_2 = 0;
        shapeFieldEdit.Scale_2 = 0;
        shapeFieldEdit.IsNullable_2 = true;
        shapeFieldEdit.Required_2 = true;
        shapeFieldEdit.Editable_2 = true;
        shapeFieldEdit.Domain_2 = null;
        shapeFieldEdit.DomainFixed_2 = true;
        shapeFieldEdit.GeometryDef_2 = BuildGeometryDefinition(BuildGenericSpatialReference(true), esriGeometryType.esriGeometryPolygon, true, false);

        var fieldNames = this.ImportedField.CoordinateData.ScalarData.Select(nodeData => nodeData.Name).ToList();
        fieldNames.AddRange(this.ImportedField.CellData.ScalarData.Select(cellData => cellData.Name));

        //Add fields from ImportedField
        foreach (var fieldName in fieldNames)
        {
            IFieldEdit fieldEdit = new FieldClass();
            fieldEdit.Length_2 = 1;
            fieldEdit.Name_2 = fieldName;
            fieldEdit.Type_2 = esriFieldType.esriFieldTypeSingle;
            fieldsEdit.AddField(fieldEdit);
        }

        return fields as IFields;
    }


texts: Calculate Average Slope over Specific Distance

texts:
I'm working with a development ordinance in which building cannot occur on areas where the average slope over a 50 foot distance exceeds 12% and these areas need to be mapped.

Calculating and categorizing the >12% areas on the DEM or CAD surface is straightforward, however it's determining where the 50' rolling average slope tapers to less than 12% that I'm having trouble with.

For example, if a 43-foot run with a slope of 5% steepens into a 7-foot run with 58% slope, this entire area along this 50 feet is considered a no-build area because the average is 12.2%.

I have the ability to work in QGIS, GRASS, SAGA, ArcGIS 10.3, or Civil 3D but I just don't know the workflow.


texts: Creating map using ArcGIS Online or ArcGIS for Desktop?

texts:
With what kind of tool/programm would you create such a map:

important criteria:

label and points need to be connected with a line
labels should have a bold headline
background color of labels should be changed individually (different layers only if not possible otherwise
the Map should automatically update edits in the excel file (.xls oder .xlsx) (attribte table are the Positions "Pos. 1 ... n" in ligns, attributes such as Gitter, FW,Holzpfosten etc und lon,lat in columns. In the cells is the amount of each attribute for the specific Position
last: particular own georeferenced maps should be added

Currently, I'm working with QGIS, where points 2,3,4 are not possible and bullet point 1 not as I wish (an extra layer is needed). Apart from that, it seems to be very userfriendly to me and is my major choice by now.
map created with QGIS

Personally, I'd prefer ArcGIS, because this is what I'm used to. But it is to complicated for the company I work for thus they only need it to label some points in the map. And they should be able to work with the map individually after I will have left the company (doing an internship here).
So I was thinking about alternatives such as ArcGIS online, ArcGIS for private user (but we use Macs here). I even tried google maps, but this doesn't work either.

What would you recommend?


texts:
I have two layers: lines and points. Points are snapped to lines.
The task is to split lines on the intersections with points.

In general, tool called Split Lines at Point works fine, as you can see on a picture (fig 1 and 2), but sometimes it makes mistakes like on fig.3 i.e. line isn't split.
I tried to specify the search radius but it didn't help.
What needs to be done in this case?


texts:
I'm using a leaflet plugin to prform on the flr geocoding of address fields. When I used a callback function to get the results, the fuction doesn't execute. Below is a sample code,

var googleGeocodeProvider = new L.GeoSearch.Provider.Google(),
  addressText = 'Amsterdam';

googleGeocodeProvider.GetLocations( addressText, function (data) {
console.log(data)

});

Getlocations is not getting executed and I'm not getting anything in console.log. Are there any solutions and can anyone help on this?


texts: Delete table rows in non-spatial database using arcpy

texts:
Working on a python script that does a bit of spatial analysis on a few layers, a spatial join of those layers, and then pulls data from a couple of other tables and holds them all in an in_memory table to be written out at the end. The idea is that this script would be run on a schedule every morning.

After all the processing is done, the in_memory table needs to be written to an existing table in a non-spatial database to be consumed by a third-party application.  I also need to be able to clear out that table before the new data is written.

I have found that UpdateCursor, Truncate, DeleteRows etc. don't work for me in this other database as the table is not registered with Geodatabase.  I am unable to do this as it's not stored in a geodatabase and I can't turn it into one.

How can I delete rows in this table in the non-spatial third-party database using arcpy?  Once the table is empty I can use arcpy.da.InsertCursor to insert new rows without any issue.

ArcGIS Desktop 10.3
Both the non-spatial and spatial database servers are running SQL Server 2012

# Delete/Truncate/Delete Rows of third-party DB table
nonGeoDB = "Database Connections/NonGeoDB.sde/NonGeoDB.dbo.ALoad"

# Delete Rows doesn't work
arcpy.DeleteRows_management(nonGeoDB)

# Update Cursor with deleteRows doesn't work
with arcpy.da.UpdateCursor(nonGeoDB, 'SEQ') as delcursor:
    for row in delcursor:
        delcursor.deleteRow()

  001323 : Input rows must be stored as a registered object in a geodatabase.

EDIT The main reason I'm trying to use arcpy over pyodbc is that I can then just use read-made database connection files that are stored securely on the server so that updates can be made to those connection files without having to modify connection settings within the .py directly.  If pyodbc can read the connection settings from one of those files, then that possibly would be the best option.


texts: I have a shapefile containing paths that I want to add to the line file prior to creating my OSM network using ArcGIS OSM Editor. How should I do this?


texts: Transform lat long to state plane coordinates

texts:
I am trying to convert a set of lat long to state plane coordinates using python and pyproj. The location is Tennessee and this gave me information about the central meridian, lat origin and the proj type which are 35.2500,34.333 and lambert conformal conic respectively. When i use National Geodetic survey's tool to do the conversion for a sample long lat(-86.911,36.316), i get this result:

NORTH(Y)    EAST(X)     AREA

220358.963  518192.337  TN

When i try doing this by code, the results are different:
Here is the piece of code:

from pyproj import Proj

p = Proj(proj='lcc',zone=4100,lon_0='35.35w',lat_1='34.33n')
x,y = p(-86.911,36.316)
print x
print y

Result x,y is:
-4367838.06379
5640715.01601

I cannot understand where i am going wrong. Please advice.


texts:
I have an output file data.json

 [{"5403": 2, "9734": 1, "2565": 5, "5087": 2, "9061": 1, "1761": 4, "4811": 7 "1111": 5, "3891": 3, "1962": 3}]

and a shapefile converted to geojson https://www.dropbox.com/s/uou0jxhpklgt8r6/pc4.json?dl=0

I use this code to generate the map in ipython notebook:

import folium
import numpy as np
from IPython.display import IFrame
from IPython.core.display import HTML

def inline_map(m, width=650, height=500):
    """Takes a folium instance and embed HTML."""
    m._build_map()
    srcdoc = m.HTML.replace('"', '&quot;')
    embed = HTML('<iframe srcdoc="{}" '
                 'style="width: {}px; height: {}px; '
                 'border: none"></iframe>'.format(srcdoc, width, height))
    return embed

    mapa = folium.Map(width=900, height=600, zoom_start=7, location=[52.3, 5], tiles='Mapbox Bright')

mapa.geo_json(geo_path='pc4.json',
              data_out='data.json',
              data=pcs,
              columns=["PC4","Count"],
              key_on="feature.properties.PC4",
              fill_color='PuRd',
              fill_opacity=1,
              line_opacity=1,
              legend_name='Demand',
              reset="True")

mapa.create_map('test.html')
IFrame('test.html', width=1000, height=700)

It creates the map fine, just as the output file, but it doesn't overlay any heatmap. It does show the legend though, but nothing else.

I know this question is similar to Geopandas/folium map not displaying but his solution doesn't do anything for me.


texts:
I have a layer with polygons - rectangles to be specific. For every rectangle there is an image specified. That's basically ImageMosaic format but ImageMosaic does not have feature info on those rectangles and I need it to be able to filter out specific rectangles.

I'd like to render a specific image for every rectangle feature and want this image to be stretched to the rectangle boundaries. This gives me repeated images in the rectangle

<GraphicFill>
  <Graphic>
    <ExternalGraphic>
      <OnlineResource
        xlink:type="simple"
        xlink:href="http://image.url" />
      <Format>image/png</Format>
    </ExternalGraphic>
  </Graphic>
</GraphicFill>

while I need one stretched image like:

<VendorOption name="graphic-resize">stretch</VendorOption>

No idea how to model this, maybe somehow with labels?


texts: ArcPy for loops increase duration with each iteration

texts:
I am running an ArcPy script with IDLE 2.7

For each line in a geometric network, it traces the path to a fixed location from that input line.
For each input line, it does this trace with all other lines as a barrier one at a time.

So 400 nodes will each perform 399 traces. Or 159,600 traces

I have found that each loop takes longer, on average, than the previous.  When I started the testing, I saw that they were taking ~2.5 min to do the 399 traces for each input node. This put me at ~17 hours to finish.  Now I'm seeing a steady increase in completion time for each complete loop of the 399 traces.

I'm at the 55th iteration right now and it's pushing 6 min now to complete each iteration.

this trend projects the 400th iteration to take 23 minutes!  So, on average, I estimate 12 min each * 400 iterations = 80 hours to run!

The chart below shows the completion time of each loop vs the loop number.  it seems to be a pretty steady increase.

The number of processes does not increase with each loop; it does the same thing, just for 400 different points.

I don't think I'm looking for an improvement to the code right now, but rather as a non-computer-scientist,  I'm wondering if this is simply expected functionality of how For loops execute.

Does it eat up memory that slows it down each time?

weight = 1
mains_cur =arcpy.SearchCursor(mains)
for main_cur in mains_cur:
    start_timeloop1 = datetime.datetime.now()
    ID = int(main_cur.getValue(("OBJECTID")))
    whereclause = ("OBJECTID" + ' = ' + repr(ID))
    whereclause_NOT =  ("OBJECTID" + ' <> ' + repr(ID))
    arcpy.SelectLayerByAttribute_management(mains_Lyr,"NEW_SELECTION", whereclause)
    arcpy.FeatureVerticesToPoints_management(mains_Lyr, r"in_memory\traceflag_A", "MID")
    arcpy.Merge_management([r"in_memory\traceflag_A",source],r"in_memory\traceflag")
    traceflag = r"in_memory\traceflag"
    arcpy.SelectLayerByAttribute_management(mains_Lyr,"NEW_SELECTION", whereclause_NOT)
    arcpy.MakeFeatureLayer_management(mains_Lyr, r"in_memory\NOTmains_Lyr", "","","")
    NOT_mains = r"in_memory\NOTmains_Lyr"

    # For each OTHER pipe, make a barrier at mid point

    not_mains_cur =arcpy.SearchCursor(NOT_mains)
        for NOT_main in not_mains_cur:
        ID2 = int(NOT_main.getValue(("OBJECTID")))
        whereclause = ("OBJECTID" + ' = ' + repr(ID2))
        arcpy.SelectLayerByAttribute_management(NOT_mains,"NEW_SELECTION", whereclause)
        arcpy.MakeFeatureLayer_management(NOT_mains, r"in_memory\not_main4point")
        arcpy.FeatureVerticesToPoints_management(r"in_memory\not_main4point",
                     r"in_memory\Not_main_pt", "MID")
        bar = r"in_memory\Not_main_pt"
        arcpy.TraceGeometricNetwork_management(geomnet, r"in_memory\outNet",
                     traceflag, "FIND_PATH", bar,"", "", "", "", "NO_TRACE_ENDS",
                     "TRACE_INDETERMINATE_FLOW", "", "","AS_IS", "", "", "", "AS_IS")
        for layer in arcpy.mapping.Layer(r"in_memory\outNet"):
            if str(layer)== "Water_Mains":
                calcstr = "["+ str(fieldname)+"]"  +r"+" + str(weight)
                arcpy.CalculateField_management(layer, fieldname, calcstr, "VB", "")
    end_timeloop1 = datetime.datetime.now()
    print('Duration: {}'.format(end_timeloop1 - start_timeloop1))

Edit:
in the above, this appears to be what really adds time to each step:

for layer in arcpy.mapping.Layer(r"in_memory\outNet"):
        if str(layer)== "Water_Mains":
            calcstr = "["+ str(fieldname)+"]"  +r"+" + str(weight)
            arcpy.CalculateField_management(layer, fieldname, calcstr, "VB", "")

I hypothesize that even though it is only adding +1 to the traced records, after a while, it's adding +1 to existing values of 100,000+  storing and calulating with progressively larger numbers may lead to more writes and time? the field it's calculating to is formatted as long.


texts:
Goal- Building python script to compute monthly average from 8 days MODIS data

I have MODIS data which is contain 8 days product. Some of month has 3 images and some of month has four imgaes. So how i can build a  Arc GIS python scrip to calculate monthly average from this data

I have listed below number of images in each month.

Jan - 4, Feb- 4, Mar- 4, Apr- 3, May- 3, June-4, July- 4, Aug - 4, Sept - 4, Oct - 3, Nov - 4, Dec - 4,

Total 45 images are in same foldar with 8 days interval name

e.g. MOD15A2.MRTWEB.A2005337.005.Fpar_1km, MOD15A2.MRTWEB.A2005345.005.Fpar_1km

After calculation i just i want to save it in month wise name.

I have Build one code but i can take only defined interval e.g 3 days or 4 days. not 3 days and 4 days at a same time.

Below i have shown my python code

import arcpy, os, sys
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
arcpy.env.extent = "MAXOF"
env.workspace = 'D:\MODIS-NDVI\FPAR-2005\MASKED-FPAR-05B'
rasters = arcpy.ListRasters()
out_ws = 'D:\MODIS-NDVI\FPAR-2005\AVG-FPAR-05'

rastersCount = len(rasters)
counter = 0
dek = 1
while counter < rastersCount:
    dekad = []
    for i in range(counter,(counter+2)):
        print i
        dekad.append(rasters[i])
        outCellStatistics = CellStatistics(dekad, "MEAN", "NODATA")
        outRasterName = out_ws + 'tmax_dek_{:03d}.tif'.format (dek)
        #outRasterName = "00" + str(dek) + ".tif"
        outCellStatistics.save(outRasterName)
    counter += 2
    dek += 1

Above can take only defined interval e.g 3 days or 4 days. not 3 days and 4 days at a same time.


texts: How to measure the sum of the areas of each class in QGIS?

texts: The project has a layer with several groups of polygons with different weights. Each group is associated to a weight. How can I get the sum of the area of each group?


texts: Google Earth Engine: cannot read "undefined" of undefined

texts:
Yesterday I was able to run the next code but today is giving me an error such as: "cannot read "undefined" of undefined"

//Choose country using GEE Feature Collection
var region = geometry;

// Collect data and filter by total dates
var modisNDVI = ee.ImageCollection('MODIS/MOD13Q1');

print(modisNDVI)

// Image collections for NDVI for all years, one value per month.
// Here is where you can change the dates.
var collection05 = ee.ImageCollection(modisNDVI.filterDate('2002-01-01', '2016-12-31'));
var collection01 = ee.ImageCollection(modisNDVI.filterDate('2008-01-01', '2008-12-31'));

//Clip to Specified Region
var clipped05 = collection05.mean().clip(region)
var clipped01 = collection01.mean().clip(region)

// Charts //
//Long-Term Time Series
var TS5 = Chart.image.seriesByRegion(collection05, small_region,  ee.Reducer.mean(),'NDVI', 500, 'system:time_start').setOptions({
          title: 'NDVI Long-Term Time Series',
          vAxis: {title: 'NDVI'},
});
print(TS5);

//Short-Term Time Series
var TS1 = Chart.image.seriesByRegion(collection01, small_region,  ee.Reducer.mean(),'NDVI', 500, 'system:time_start').setOptions({
          title: 'NDVI Short-Term Time Series',
          vAxis: {title: 'NDVI'},
});
print(TS1);

var mult05 = clipped05.multiply(0.0001)
var mult01 = clipped01.multiply(0.0001)

//Set Center of Map & Add Clipped Image Layer
Map.setCenter(-50, -10, 3); //long, lat, zoom
Map.addLayer (mult05.select('NDVI'))
Map.addLayer (mult01.select('NDVI'));

// Create a geometry representing an export region.
var roi = ee.Geometry.Polygon(
    [[[-68.675537109375,-16.799282416303967],
      [-68.79638671875,-17.45023229881414],
      [-68.291015625,-17.539296556236692],
      [-68.126220703125,-16.92019465044387],
      [-68.675537109375,-16.799282416303967]]]);
// Export the image, specifying scale and region.
Export.image.toDrive({
  image: clipped01,
  description: 'NDVI_MODIS_2008',
  scale: 250,
  region: roi
});

Does anyone know where is the mistake? It can be a server issue?


texts: Are there open-source alternatives to ArcGIS online that can be hosted on my own server?

texts:
I know this might be too broad topic. But I am looking for web-GIS portal similar to ArcGIS online that is self-installable (configurable) to be hosted in my private server. The main features should include:

Users can upload and interact with spatial data (shp, CSV, KML, JSON etc)
Users can rate, comment and download layers submitted by others

In my opinion, QGIS-server, QGIS-web-client and/or QGIS cloud could develop to be a full featured cloud service that can be deployed on any server.


texts: How I can incorporate results of histogram stretching into georeferenced image?

texts: I have one channel 8 bit raster and when I open it in QGIS with histogram stretching "on the fly" it looks fine. How I can recalculate original raster accordingly changed histogram?


texts:
I am new to Ubuntu, so forgive me if this a very standard question:
I installed QGIS 2.18 and all seems to work well. However, there is a major issue: the program uses a very strange font so I can't read anything? I tried to find something on the internet, but none seem to have the same issue.

As I do not now much about Ubuntu, my trouble-shooting-abilities are very low.

Is there a workaround to change the font in QGIS so I can get started working with it?


texts: ITopologicalOperator.Intersect, how does this work?

texts:
I'm trying to intersect all features of a layer with another's but I'm getting the web ADF application stuck while executing the "ConstructUnion" function

Here's the code I've written so far

// create a geometrybag for layer1
ESRI.ArcGIS.Geometry.IGeometry geometryBag = new ESRI.ArcGIS.Geometry.GeometryBagClass();
geometryBag.SpatialReference = m_firstLayerFirstfeatureIGeometry.SpatialReference;
ESRI.ArcGIS.Geometry.IGeometryCollection geometryCollection = (ESRI.ArcGIS.Geometry.IGeometryCollection)geometryBag;
// Cycle all features of first layer and add them to the collection
foreach (DataRow m_dr in m_firstLayerfeatureGraphics.Rows)
{
    // Takes this feature's geometry
    ESRI.ArcGIS.ADF.Web.Geometry.Geometry adfFeature =
        m_firstLayerfeatureGraphics.GeometryFromRow(m_dr) as ESRI.ArcGIS.ADF.Web.Geometry.Geometry;
    // convert to a igeometry object
    ESRI.ArcGIS.Geometry.IGeometry m_tempfeatureIGeometry =
    ESRI.ArcGIS.ADF.Web.DataSources.ArcGISServer.Local.Converter.ToIGeometry(adfFeature, ctx);
    System.Diagnostics.Debug.WriteLine("m_firstLayerfeatureGraphics +1 ");
    object missing = Type.Missing;
    geometryCollection.AddGeometry(m_tempfeatureIGeometry, ref missing, ref
    missing);
}
ESRI.ArcGIS.Geometry.ITopologicalOperator mtopo1 = new ESRI.ArcGIS.Geometry.PolygonClass();
// FOLLOWING LINE BLOCKS APPLICATION
mtopo1.ConstructUnion((ESRI.ArcGIS.Geometry.IEnumGeometry)geometryBag);

Any idea why?


texts: Calculate distances between series of points in postgis

texts:
I have a postgres/postgis table containing a series of sequential points. I want to add a new column to the table which contains the distance between each point and the first one, and so forth until all the points are processed (the first point should have a distance of 0). Can this be done with a query or do I need to use Python or equivalent to create a loop?
Many thanks


texts: I have a geodatabase with all of my data, but I would like to add more friendly aliases from a separate key I have stored in a CSV. How do I go about adding these aliases to my field names in my geodatabase?


texts: Simplify Geometries 'NoneType' Error in QGIS

texts:
I need to simplify the geometries of a large amount of polygons, yet every time I run the Simplify Geometries tool in QGIS, I get the error message below.

'NoneType' object has no attribute 'type' See log for more details

What is causing that error? How can I fix it?


texts: Preserving text when converting from .dxf to .shp?

texts:
I have tried using dxf2shp converter in QGIS but the text in dxf file is being lost during the conversion process, despite selecting export text labels.

How do I convert the .dxf to .shp file without losing data?


texts: Make a grid of 5*5km on a raster layer map in R?

texts:
I'm trying to make a grid of 5 by 5km on my raster layer map but not sure how to do this. I have read my shapefile in R and rasterized it and have an resolution  of 10 by 10 meter. And I have plants distribution coordinates on this raster map. Now I want to make a grid of 5 by 5 km to calculate the number of plants in each grid.

But stuck with how to make a grid on this raster layer.

I followed the suggested procedures but i managed to plot this new raster separately and could only plot the plant distribution coordinates on my base raster map and I could not overlay this new layer with 5 by 5km on top to extract or find out the plants distribution per area. I am not sure what i am doing wrong.
below is what i did:

landuse <- readShapeSpatial("~/Documents/vegetation_analysis/PNPP.shp",
                        proj4string = CRS("+proj=utm + datum=WGS84"))

p <- raster(extent(landuse)) #creats a raster with the extent of VEG in area PNPP
projection(p) <- proj4string(landuse) #uses the projection of the shapefile
res(p) <- 100.0 #sets the resolution of the raster

p1 <- rasterize(landuse, field="TYPES", r) #converts BHT to a raster
plot(p1) ##plot this one

#loaded my vegetation locations
veg_d <- SpatialPointsDataFrame(coords = cbind(VD$long, VD$lat),
                               data = VD,
                               proj4string = CRS("+proj=longlat +datum=WGS84")
)

veg_d <- spTransform (veg_d, CRS("+init=epsg:4240"))

#make a new occurance layer as suggested here
t <- raster(xmn=0, ymn=0, xmx=100000, ymx=100000, resolution=10,
            crs='+proj=utm +ellps=WGS84')
t <- setValues(t, sample(x=0:1, size=ncell(t), replace=T))
plot(t)

##then raster aggregation with aggregate
t_agg <- aggregate(x=t, fact=400, fun=sum)
t_agg
plot(t_agg)
points(veg_g,col="blue") ##does not make the points on the plot????

I also plotted:

plot(landuse)
points(veg_g, col="blue")
plot(t_agg, add=TRUE)  ##

but i don't plot the raster over this and i am not able to get it.

I am new to R.


texts:
In GeoNetwork/ISO 19115, how do you describe physical access parameters for a physical paper map (i.e. of type CI_PresentationFormCode==mapHardcopy)?

For example I want to be able to find the map through a room number, archive box and signature number associated to the map. Also interesting would be to store the size (A0, A1, ...) or weight.

Storing online resources works quite well with GeoNetwork and is good documented, for example you can enter a online location (URL) via CI_OnlineResource, But what about physical parameters of "hardware" maps? Are there other standards than ISO19115 that are better suited for physical maps?

Clarification: I'm aware that in ISO-19115, one can define that a data IS a paper map by setting CI_Citation.presentationForm to mapHardcopy (an element from the CodeList CI_PresentationFormCode). But: If my data is a paper map, how do I describe properties of
it that are very specific to paper maps:

weight
size
book shelf
material health
library number

The only thing I can imagine is the library number which is storable in CI_Citation.identifier. What about the others?


texts: Comparing table attributes in ArcGIS 10(ArcInfo)

texts: I'd like to compare attributes of two point layers . the layers have some similar records and i try to find these . ArcGIS 10.1 has a tool to comparing attributes but i could'nt find any similar tools in ArcGIS 10 . How can i compare two attribute and find similar records?


texts: Reproject numpy array based on bbox (extent)

texts:
I have a NumPy array image which is basically NDVI values.
The image looks like this:

I know the bbox coordinates of this image (e.g (bbox_coords_wgs84=[minx, miny, maxx, maxy]))
I'm looking for python functions to use in order to reproject easily this image based on the bbox coordinates. I sawnrasterio has functions but it looks a bit complicated for the situation I have already the bbox.


texts: Can I import a .gml file into Spatialite in QGIS?

texts:
I have some Ordnance Survey  Address Base data files which have a .gml extension. I have imported them to GGIS as vector layers along with some mapping data but would like to explore them with PostgreSQL. I see that QGIS comes with SpatiaLite installed. Can I import the gml files into SpatiaLite or for that matter PostgreSQL?

I suspect the format and content of my question says more about what I don't know so any comments or even answers would be appreciated.


texts: Selecting raster data by location: for raster and feature class data

texts:
I have rasters covering certain spots in my environment, I'm looking for a way to extract the values which occur at locations which are defined by a selection of feature layers. I've tried with the SelectLayerbyLocation tool (from ArcMap 10.2 (advanced license)), which seems to do the right thing, but not for this type of data. Does anybody have experience with this or recommendations for tools to use?

Currently it is set up like this:

arcpy.SelectLayerByLocation_manager('codist_{0}'.format(FIDa), "INTERSECT", out_flab).save(overlap) #find out if villages overlap with cost layers

where 'codist_{0}'.format(FIDa) is the raster (FIDa a list of values), and out_flab is the feature layer extracted from a table using a cursor in previous steps. Ultimately I need to find the minimum value within the rasters at the zones defined by the feature classe (so the values of the raster cannot be altered by this operation)


texts: Cut hole in raster with shapefile and sample at closest edges

texts:
I have a raster in which I'd like to cut a "hole" with a shapefile. I have points close to the boundary, but inside the hole, and I would like to find the closest point on the edge. To give a concrete example, I have a global dataset (1' resolution) and a shapefile of the coast of Africa. My data points are locations of cities near the coast, but some of these points are not exactly on the coast, so I'd like to find the closest point outside the shapefile, i.e. closest point over water.

Normally, I'd use: gdalwarp -of -GTiff -cutline Coast.shp -crop_to_cutline -dstnodata 0.0 Data.tif Clipped_Data.tif but this obviously only gives me what's inside the shapefile; I can't find an inverse clipping for this command.

I did some digging and found this, which allowed me to set everything inside the shapefile to 0, but when I sample using gmt grdtrack it just sets the values to 0 (which makes sense). Also, when I tried crop=True, invert=True I got: ValueError: crop and invert cannot both be True. So I reverted to using crop=False, invert=True.

So, how can I sample this raster at the closest point outside the shapefile? I think if I can get the hole to be set to NaN, I can use the -T flag in gmt grdtrack.


texts: How to extact very narrow water linear features from GeoEye-1 image

texts: I need to extract beaver channels from a GeoEye image (attached image). What is the best approach to automatically delineate such very narrow and not very sharp water features?


texts: Creating arrays and polylines from dictionary using ArcPy?

texts:
I have created a dictionary out of a csv values from an excel spreadsheet. I'm trying to use the keys and values in the dictionary to create arrays (one per 'Rhino') that I can use to create polylines.
The code I've writen so far has got me stuck when creating the arrays, since when I print the coordinate pairs in the array it gives me the same value repeated many times.

I am still in the beginning of learning Python for ArcGIS.

#reads through a text file with coordinates of rhino tracks,
#extracts the locations and adds pairs of coordinates to a polyline as an array
#to be automated with a function

import csv
import arcpy

#define a function to be used to add rhinos and their tracks to a dictionary
def rhinoTrackFinder(rhino, coordX, coordY, dictionary):
    coordPair = [coordX, coordY]
    if rhino in dictionary:
        dictionary[rhino].append(coordPair)
    else:
        dictionary[rhino] = coordPair

#open text file with rhinos and coordinates
trackFile = open("C:/Rafa/Lesson4PracticeExercises/RhinoObservations.csv")

#read the header line and create indices, then create empty dictionary
#that will store values and keys for rhinos and the coordinates of their   tracks
csvReader = csv.reader(trackFile, delimiter=",")
header = csvReader.next()

surveyorIndex = header.index("Observer")
XCoordIndex = header.index("X")
YCoordIndex = header.index("Y")
RhinoIdIndex = header.index("Rhino")
CommentIndex = header.index("Comments")

rhinoTrackDictionary = {}

for row in csvReader:
    rhino = row [RhinoIdIndex]
    CoordX = row [XCoordIndex]
    CoordY = row [YCoordIndex]
    #new seg - boolean to define if it is a new polyline or rhino
    rhinoTrackFinder(rhino, CoordX, CoordY, rhinoTrackDictionary)

#print rhinoTrackDictionary

#SO FAR SO GOOD#

#create an array to store lat and lon values of the dictionary as points

vertexArray = arcpy.Array()

#iterate through dictionary to create and add vertices as array points

for rhino in rhinoTrackDictionary:
    vertexArray = arcpy.Array()
    for pair in rhino:
        lon = CoordX
        lat = CoordY
        vertex = arcpy.Point(lon, lat)
        print vertex
        vertexArray.add(vertex)

        #def createPolyline (vertexArray,)

for pair in rhinoTrackDictionary:
    for vertex in pair:
        vertex = arcpy.Point(lon,lat)
        vertexArray.add(vertex)


texts: Extracting field values to Excel by iterating in ArcGIS ModelBuilder?

texts:
I have a database with a lot of feature classes containing shapefiles.
What I want to do is to use ModelBuilder in order to iterate through all the feature classes and for each feature class go through the point shapefiles to extract a field value. The field value is then to be exported to an Excel file.

You can see my workflow in ModelBuilder in the attached image.
An Excel file is created but it doesn't contain any data at all, only the column headers.

How can I fix this?


texts:
Given a vector line and symbology settings as a simple dashed line (in QGis v3.4.3):

As I understand the QGIS concepts there is a difference between a line (feature) contained on a vector layer and its representation as defined in the layer properties (under Symbology).
But I would like to get the outline of this dashed line (as configured in the symbology of the vector layer) as a set of polygons (may be multipolygons or similar) on a new created vector layer.

For a solid line I would use the tool buffer from the vector geometry toolbox to "convert" a vector line to a polygon representation of an outline with a certain width.

Is it possible to convert/transform/"bake" (whatever the name is) a symbology (out-)line into a vector layer like polygons or multipolygons?


texts:
I last used QGIS 2.14.9 Essen in September 2018 with no problems. Today 28/11/18 I tried to open Google Maps plugin (and Bing maps etc.) but I am getting the following error message - OpenLayers Plugin: Printing and rotating of JavaScript API based layers is currently not supported

I have reinstalled the Open Layers plugin to no avail.


texts:
With QGIS 2.14.3 I am attempting to convert a raster layer to a binary layer, using a certain threshold. With some layers I can use this formula in the Raster Calculator and get the expected output;

"raster@1" < 0.994 = 0 AND "raster@1" > 0.994 = 1

But when I try it with another layer the output is either all nan or all 0's. I've tried saving it in a different format (.img instead of .tif) as suggested in another answer here, but get the same result.

Is there a mistake I'm making with this?

Or an obvious reason it might give this output?

Or maybe the raster is just corrupted in some way?

(The imagery is the output from a long process, and confidential, so I can't share the actual file).


texts:
I am using liblas 1.8.1 and laszip 3.1.1 on Mac OS High Sierra. Both were installed using brew install <package>. I am looking to get information from a laz file using the following command:

lasinfo lazfile.laz

Which gives me the following error:

error: Compression support not enabled in liblas configuration

Assuming liblas allows working with laz files, how can I resolve this error and work with laz files in liblas?


texts:
I'm working in a urban biodiversity project and I have a csv file for a birds inventory with a column index called 'geom' composed by multypolygons as this:

        name    geom                                                id
0   Abubilla    MULTIPOLYGON(((421405.897134601 4585295.726121...   83
1   Abubilla    MULTIPOLYGON(((424405.897979997 4584795.702582...   83
2   Abubilla    MULTIPOLYGON(((424905.87219322 4581295.6934019...   83
3   Abubilla    MULTIPOLYGON(((424905.883573821 4582795.695714...   83
4   Abubilla    MULTIPOLYGON(((424905.906345894 4585795.700331...   83

the column 'geom' is not recognized as geometry.

when I run

birds_geo = birds.set_geometry('geom')

it returns the following error:

TypeError: Input geometry column must contain valid geometry objects.

I've tried this script which I found for a similar question:

from shapely.wkt import loads
birds['geom'] = birds['geom'].apply(lambda x: loads(x))
for index, row in birds.iterrows():
    birds.set_value(index, 'geom', loads(row['geom']))

but it seems to work for polygons but not for multipolygons, returning the following error:

AttributeError: 'MultiPolygon' object has no attribute 'encode'

Is there a way to make to turn multipolgons into a valid geometry index?


texts:
I want to write C# code for utility network analyst (downstream trace task).

My code is:

IUtilityNetworkAnalysisExt m_utilNetExt = new UtilityNetworkAnalysisExtClass();
ITraceFlowSolverGEN tfs = (ITraceFlowSolverGEN)new TraceFlowSolverClass();
IGeometricNetwork pGNetwork = frmNetwork.pGeometricNetwork;
INetworkAnalysisExt nax = (INetworkAnalysisExt)m_utilNetExt;
INetwork net = pGNetwork.Network;
INetSolver netSolver = (INetSolver)tfs;
netSolver.SourceNetwork = net;
INetworkAnalysisExtFlags naxFlags = (INetworkAnalysisExtFlags)m_utilNetExt;

#region set Flag display from feature select
IJunctionFlagDisplay junctionFlagDisplay = new JunctionFlagDisplayClass();
IEdgeFlagDisplay edgeFlagDisplay = new EdgeFlagDisplayClass();
IFlagDisplay flagDisplay;
IEdgeFlag edgeFlag=new EdgeFlagClass();
IJunctionFlag junctionFlag = new JunctionFlagClass();
INetFlag netFlag;
IEdgeFlag[] edgeFlags = new IEdgeFlag[1];
IJunctionFlag[] junctionFlags = new IJunctionFlag[1];

if (pFeature != null && pFLayer.FeatureClass.ShapeType == esriGeometryType.esriGeometryPoint)
{
    flagDisplay = (IFlagDisplay)edgeFlagDisplay;
    flagDisplay.FeatureClassID = pFLayer.FeatureClass.FeatureClassID;
    flagDisplay.FID = pFeature.OID;
    flagDisplay.SubID = -1;
    edgeFlagDisplay =(IEdgeFlagDisplay) flagDisplay;
    edgeFlag.Position = Convert.ToSingle(edgeFlagDisplay.Percentage);
    netFlag =(INetFlag) edgeFlag;
    netFlag.UserClassID = flagDisplay.FeatureClassID;
    netFlag.UserID = flagDisplay.FID;
    netFlag.UserSubID = flagDisplay.SubID;
    edgeFlags[0] = edgeFlag;
    tfs.PutEdgeOrigins(ref edgeFlags);
}

if (pFeature != null && pFLayer.FeatureClass.ShapeType == esriGeometryType.esriGeometryPolyline)
{
    flagDisplay = (IFlagDisplay)junctionFlagDisplay;
    flagDisplay.FeatureClassID = pFLayer.FeatureClass.FeatureClassID;
    flagDisplay.FID = pFeature.OID;
    flagDisplay.SubID = -1;
    junctionFlagDisplay = (IJunctionFlagDisplay)flagDisplay;
    //junctionFlag.Position = Convert.ToSingle(junctionFlagDisplay.Percentage);
    netFlag = (INetFlag)junctionFlag;
    netFlag.UserClassID = flagDisplay.FeatureClassID;
    netFlag.UserID = flagDisplay.FID;
    netFlag.UserSubID = flagDisplay.SubID;
    junctionFlags[0] = junctionFlag;
    tfs.PutJunctionOrigins(ref junctionFlags);

}

ITraceTasks traceTasks =(ITraceTasks) m_utilNetExt;
tfs.TraceIndeterminateFlow = traceTasks.TraceIndeterminateFlow;
IEnumNetEID resultJuncs, resultEdges;
esriFlowElements flowelements = traceTasks.TraceFlowElements;
tfs.FindFlowElements(esriFlowMethod.esriFMUpstream, flowelements, out resultJuncs, out resultEdges);

INetworkAnalysisExtResults naxResults =(INetworkAnalysisExtResults) m_utilNetExt;
naxResults.ClearResults();
naxResults.CreateSelection(resultJuncs, resultEdges);

#endregion

This code does not work, please help me?


texts: Matplotlib skip color for assigned interval in map plotted from GeoPandas DataFrame, if attribute don't have value in range of user defined bins

texts:
I have a GeoPandas GeoDataFrame which I am plotting against a specific column and the colour of the output/map is same for different assigned bins/intervals.
The issue is because, sometimes, there is no value in/between the range assigned in the user defined scheme for the column/variable which we are plotting, like here
I have given 4 ranges
0 - 0.25,
0.25 - 0.50,
0.50 - 0.75,
0.75 - 1.00

but the variable/column from dataframe which I am plotting does not have any value in the range 0.75 - 1.00.
So it is assigning the same colour to the legend of  the intervals 0.50 - 0.75, 0.75 - 1.00,

But I want to keep the bins values and legends static, as if, it should show the legends of 4 different colours only,

and it should ideally skip that interval/colour while plotting inside the polygons of the map if there is no value in that interval/range,
and it should still shows that interval/range in the legend.
Like, if here, there was no value between 0.75 to 1.00, then it may skip the green colour from the listed colour map assigned for that specific bound to plot inside the polygons,
but still shows that range/interval in the legend.
cmap = ListedColormap(["red", "orange", "blue", "green"])

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

Also, if you will see the script, and look closely at the legends of 1st figure in my question, the color sequence which I have given in cmap = is in the order "red", "orange", "blue", "green",
but when the value in range (0.75 - 1.00) is missing, it is skipping the orange color and jumping to the blue (3rd color in the sequence) instead of picking the orange (2nd color in the sequence), i.e., not following a proper trend.
Is there any way possible?
for day in days_list:
    map_attribute = "far"
    df2 = csv_file[["state_id", "state_name", map_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.replace(np.nan,0)

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

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

    df_new.plot(
    ax=ax,
    column=map_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
    )


texts: I would like to create a wind map like here: www.windytv.com. It doesn't need to be animated. However, an animation would be a nice final add-on. I am working with ArcGIS, QGIS, several python libaries and R. However I have no clue where to get started and how to accomplish the task.


texts: Why are the changes to my raster attribute table not saving?

texts:
I have a tiff file with an existing raster attribute table. I am adding a column in from a csv file. From the print statements, it seems to be adding it correctly; however, it's not saving for some reason. I don't see it when I load it in ArcMap.
On the first run of the script, I expect the first print statement to produce an error and the second print statement to work. Both of which happens as expected. On the second run, I would expect both print statements to work but the first one still doesn't work indicating that it's not saving.
This is gdal 2.4.2.
edit: I put a print statement after the "band.SetDefaultRAT(rat)" line and that isn't printing so I assume that line is producing some sort of error killing the script. But it's not showing me any error in the terminal either. Not sure what I'm doing wrong.
edit2: I tried removing basically everything in between getting the rat and setting the rat (with no changes in between) and it still seems to not do anything past that band.SetDefaultRAT(rat) line. So that line seems to be causing me issues. Not sure what the issue is though.
import gdal
import csv
import numpy

csv_file = "location"
tiff_file = "location"
col_fields = numpy.loadtxt(csv_file, dtype=str, delimiter=';', usecols=(6))
ds = gdal.Open(tiff_file, gdal.GA_Update)
band = ds.GetRasterBand(1)
rat = band.GetDefaultRAT()
print(rat.GetValueAsString(15,4)) #gives error for not existing (theoretically should work on 2nd run but it doesn't)
rat.CreateColumn("Name", gdal.GFT_String, gdal.GFU_Generic)
for index in range(rat.GetRowCount()):
    rat.SetValueAsString(index, 4, str(col_fields[index]))
print(rat.GetValueAsString(15,4)) #successfully prints
band.SetDefaultRAT(rat)
ds = None


texts: ArcCatalog Representations Options Greyed Out

texts:
I'm guessing this is a permissions error.   But everything I've tried enabling hasn't fixed the problem.  Does anyone know why this is greyed out or how to fix it?

This is in ArcCatalog.


texts: Preventing SLD visualisation rendered as RGB in Earth Engine

texts:
Following this answered question, I would like to visualise the rainfall forecast using specified threshold. I tried to use Styled Layer Descriptor (SLD). But unfotunately, but it doesn't work. All pixels are converted to RGB.
How to prevent that? Because standard visualisation is hard to see the rainfall category based on milimeters value. I really want to visualise it based on SLD and rendered as milimeters.
Is there any option that I can follow?
var region = ee.Geometry(Map.getBounds(true))

// Strip time off current date/time
var today = ee.Date(new Date().toISOString().split('T')[0])
var start_period = ee.Date('2015-07-01');

var collection = ee.ImageCollection('NOAA/GFS0P25')
  .select('total_precipitation_surface')
  .filterMetadata('forecast_hours', 'equals', 24)

// UI widgets needs client-side data. evaluate()
// to get client-side values of start and end period
ee.Dictionary({start: start_period, end: today.advance(1, 'day')})
  .evaluate(renderSlider)

function renderSlider(dates) {
  var slider = ui.DateSlider({
    start: dates.start.value,
    end: dates.end.value,
    period: 1, // Every 5 days
    onChange: renderDateRange
  })
  Map.add(slider)
}

// Define an SLD style of discrete intervals to apply to the image.
var visRainForecast =
  '<RasterSymbolizer>' +
    '<ColorMap  type="interval" extended="false" >' +
      '<ColorMapEntry color="#cccccc" opacity="1.0" quantity="1" label="No Rain" />' +
      '<ColorMapEntry color="#f9f3d5" opacity="1.0" quantity="10" label="1-10" />' +
      '<ColorMapEntry color="#dce2a8" opacity="1.0" quantity="20" label="11-20" />' +
      '<ColorMapEntry color="#a8c58d" opacity="1.0" quantity="30" label="21-30" />' +
      '<ColorMapEntry color="#77a87d" opacity="1.0" quantity="40" label="31-40" />' +
      '<ColorMapEntry color="#ace8f8" opacity="1.0" quantity="60" label="41-60" />' +
      '<ColorMapEntry color="#4cafd9" opacity="1.0" quantity="80" label="61-80" />' +
      '<ColorMapEntry color="#1d5ede" opacity="1.0" quantity="100" label="81-100" />' +
      '<ColorMapEntry color="#001bc0" opacity="1.0" quantity="120" label="101-120" />' +
      '<ColorMapEntry color="#9131f1" opacity="1.0" quantity="150" label="121-150" />' +
      '<ColorMapEntry color="#e983f3" opacity="1.0" quantity="200" label="151-200" />' +
      '<ColorMapEntry color="#f6c7ec" opacity="1.0" quantity="1000" label="&gt; 200" />' +
    '</ColorMap>' +
  '</RasterSymbolizer>';

function renderDateRange(dateRange) {
  var gfs1day = collection
      .filterDate(dateRange.start(), dateRange.end())
      .median(); // make a composite of the collection

  var vis = {min: 0, max: 200, palette: [
      'cccccc','f9f3d5','dce2a8','a8c58d','77a87d','ace8f8',
      '4cafd9','1d5ede','001bc0','9131f1','e983f3','f6c7ec'
    ]}
  var layer = ui.Map.Layer(gfs1day, vis, 'GFS 1-day')
  var layer2 = ui.Map.Layer(gfs1day.sldStyle(visRainForecast), {}, 'GFS 1-day SLD')
  Map.layers().reset([layer, layer2])
}

GEE code


texts: kml GroundOverlay href of image with two '=' within the url

texts:
I wish to provide a lightweight kml (to be zipped to a kmz)
containing GroundOverlays of the NASA MODIS subsets across the Alaskan Arctic waters.  Recipients must be able to receive this kmz over a very limited bandwidth, and then select the GroundOverlay of choice
(for their field study region; cloud-free day; and desired resolution).
For years my kml worked without a glitch, however last fall NASA reorganized their url to include two "=" symbols

(example: https://lance.modaps.eosdis.nasa.gov/imagery/subsets/?project=other&subset=BeringSea.2017079.terra.721.2km.jpg).

When I pass this as the target of my kml GroundOverlay href, I get an error message indicating a "not well informed (invalid token)" at the position of the second "=" symbol within the href url in the example kml below.

Please advise with strategies to restore this GroundOverlay function to the offending url's.

Example GroundOverlay kml:

<?xml version="1.0" encoding="UTF-8"?>
 <kml xmlns="http://www.opengis.net/kml/2.2"
 xmlns:gx="http://www.google.com/kml/ext/2.2">
 <Folder>
 <name>Sea Ice Data for field opperational use only
 <br>Updated 2017 Mar 20 12:39 </br>
 </name>
 <open>1</open>
 <LookAt>
 <longitude>-163</longitude>
 <latitude>71</latitude>
 <altitude>0</altitude>
 <range>2500000</range>
 <tilt>0</tilt>
 <heading>0</heading>
 </LookAt>
 <Folder>
 <Style><ListStyle><listItemType>radioFolder</listItemType></ListStyle></Style>
 <name>MODIS  BeringSea terra </name>
 <open>0</open>
 <Folder>
 <Style><ListStyle><listItemType>radioFolder</listItemType></ListStyle></Style>
 <name>MODIS  BeringSea 2km terra </name>
 <GroundOverlay>
<name>MOD2k_Ber_2017-Mar-21</name>
 <open>0</open>
 <visibility>1</visibility>
 <color>ffffffff</color>
 <Icon>
<href>https://lance.modaps.eosdis.nasa.gov/imagery/subsets/?project=other&subset=BeringSea.2017079.terra.721.2km.jpg</href>
 <refreshMode>onChange</refreshMode>
 </Icon>
<LatLonBox id="khLatLonBoxMODIS_BeringSea">
<north>70</north>
<south>58</south>
<east>-155</east>
 <west> -190 </west>
 </LatLonBox>
 </GroundOverlay>
  </Folder>
 </Folder>
 </Folder>
 </kml>


texts: How to make a VRT with an input raster containing multi bands?

texts:
I have 3 different rasters each containing 2 bands. Such as :

Input1.tif B1 B2
Input2.tif B1 B2
Input3.tif B1 B2

I want to build a vrt with the B2 of each raster. With gdalbuildvrt I always get this warning :
Warning 1: Input1.tif has 2 bands. Only the first one will be taken into account in the -separate case

for the moment, I launch :
gdalbuildvrt -separate output.vrt Input1.tif Input2.tif Input3.tif

How can I do this ?


texts: GRASS GIS: Computing center line of raster feature representing a river

texts:
My goal:

Convert a raster layer representing a river feature (see picture) to the center line of the river stored as a polyline feature.

Limitation:

The whole process has to be automated.
I can only use GRASS GIS for this application.

The initial input is a raster resulting from a scanner. From this scan I have reclassified the river.

I have tried the following so far:

Convert to polygon --> convert to line (out-lines e.g. river banks) --> convert river banks to point --> apply voronoi diagrams and clip this output to the extent of your river --> see output. I tried using v.clean to clean the output e.g. remove dangles but this has not worked yet.

r.grow to thin my raster or r.thin to thin it --> convert thinned raster to a line --> v.clean to clean the result. This result (black lines) contains "loops" which I did not manage to delete.

Especially here "v.generalize" with a high threshold yielded some better results but still not all loops are deleted. So, you can consider this tool as well in a possible solution!

Any recommendations on a different approach or to alter my previous attempts?
In an other post here (Finding the Center Line from a set of 3D Points) they mention a possible interesting approach suggested by Joseph.

EDIT: 28.05.2019 Will work on this tomorrow again and test again my approach which is also mentioned by Micha. Jons' solution seems to be the most robust approach however with my skill and time budget I am not sure if I can reproduce a algorithm which will work in GRASS GIS.


texts:
Is it possible to extrude an an angle other than vertical in ArcScene?

I am trying to visually display some mine adits that cut between two layers at an angle of 30 degrees.


texts:
I have a plugin with QDockWidget and QActions. When I first load (or reload) the plugin, the actions' size is set according to current settings.

However, when QGIS is closed and open, the actions are smaller. The size of my actions is "one step behind" QGIS settings unless the smallest icon size is set.

Here is a part of my code where dock widget is initialized and showed/hidden.

def initGui(self):
    self.dockWidget = dockwidget.DockWidget(self.iface)
    self.iface.addDockWidget(Qt.TopDockWidgetArea, self.dockWidget)

def unload(self):
    self.iface.removeDockWidget(self.dockWidget)

def run(self):
    if not self.dockWidget.isVisible():
        self.dockWidget.show()
    else:
        self.dockWidget.hide()

Is there a way to keep the actions' size according to QGIS settings?

The only idea I came up with is to keep track of the dock widget's current area and initialize the dock widget every time QGIS starts but that is not satisfying.


texts:
Is there a method to return the length of a Cesium JS polyline, in real-world units?

For example, if I create a polyline:

var orangeOutlined = viewer.entities.add({
        name : 'Orange line with black outline at height and following the surface',
        polyline : {
            positions : Cesium.Cartesian3.fromDegreesArray([-104.606667,50.454722,
                                                        15.71666666, 69.1 ]),
            width : 5,
            material : new Cesium.PolylineOutlineMaterialProperty({
                color : Cesium.Color.ORANGE,
                outlineWidth : 2,
                outlineColor : Cesium.Color.BLACK
            })
        }
    });

... is there a method to return the length of that line, like:

console.log(orangeOutlined.length)


texts: Skipping ArcGIS Online hosted MXDs using ArcPy?

texts:
I've previously posted a question about an ArcPy tool that walks through a specified area and scans all MXDs and associated Layers and reports them (as well as the inverse)
Search for a shapefile of any kind and report all associated MXDs.

The issue is if a shapefile is used within a MXD hosted on ArcGIS Online, the script halts and ArcGIS prompts the user for login info. Until that window is closed, the script will not continue.

Can ArcPy detect this?

Can it be written into an exception and skipped?

I am familiar with exceptions and how they work.


texts:
I need to update a building data (shapes) in ArcGIS Desktop.

One Shape is with the buildings of a city from 2012 and the otherone with the buildings of the same city from the year 2019 (perhaps 300.000 buildings). I need to update the shape from 2012 while changing the use of some buildings and while removing buildings, which doesn't exist anymore and add those which are new. I need to keep the 2012 shape because it has attached much more information than the 2019 shape.
The problem is, that there is no field in the attribute table to connect both shapes, there differences between the ID and all other fields. The second problem is, that the buildings are not exactly at the same place. Some yes but most of them are just a little bit beside (already checked: it's not a question of projection).

How can I update the 2012 shape without having to manually edit every single building?


texts:
I have several SVGs, the simplest of them is the following SVG:

<?xml version=‘1.0' encoding='utf-8'?>
<svg viewBox="-337 43 2223 1812"
    xmlns="http://www.w3.org/2000/svg">
    <path d="M1508 43h-1528v1528h1528v-1528z" glyph-name="exclam" name= "foreground" fill="rgba(252,133,33,255)" transform="scale(1,-1) translate(0,-1855)" unicode="!" />

</svg>

When I try to render then in GeoServer they all fail to draw and resort to the default marker instead of the orange square I'm expecting (for the simple example).
I get a ton of messages "InterpolationConverterFactory can be applied from string to Interpolation  only", and every once in a while i get "creating defaultMark".

If I convert the SVGs to PNG it displays fine (and changing the rule format to image/png from image/svg)

My GeoServer version is 2.11.0

The SVGs display okay in Chrome.


texts: Enabling "Show all labels for this layer (including colliding labels)" in PyQGIS

texts:
I am using PyQGIS 3.10.
I am working on a project that will display about 30 cities and their forecast weather variables (temperatures, winds, etc.), and some of the labels are a bit close together and technically collide.
I can click this checkbox in the UI, but I haven't been able to find a way to code this in the Python console.
Any ideas?
Below is my code that I've been working with.
for layer in QgsProject.instance().mapLayers().values():
    layer_settings = QgsPalLayerSettings()
    text_format = QgsTextFormat()
    buffer_settings = QgsTextBufferSettings()
    shadow_settings = QgsTextShadowSettings()

    text_format.setFont(QFont("Arial Rounded MT Bold", 12))
    text_format.setSize(12)
    text_format.setColor(QColor("white"))

    buffer_settings.setEnabled(True)
    buffer_settings.setSize(0.5)
    buffer_settings.setColor(QColor("black"))

    shadow_settings.setEnabled(True)
    shadow_settings.setOffsetAngle(135)
    shadow_settings.setOffsetDistance(1)
    shadow_settings.setBlurRadius(1.5)

    text_format.setBuffer(buffer_settings)
    text_format.setShadow(shadow_settings)
    layer_settings.setFormat(text_format)

    layer_settings.fieldName = "City"
    layer_settings.placement = 1
    layer_settings.quadOffset = 7

    layer_settings.enabled = True
    layer_settings = QgsVectorLayerSimpleLabeling(layer_settings)
    layer.setLabelsEnabled(True)
    layer.setLabeling(layer_settings)
    layer.triggerRepaint()


texts:
I’m sorry, I may be missing something very basic; but I wanted to fill in some missing lat/long values in my attribute table (using QGIS 1.8). I found this:
How do I calculate the latitude and longitude of points using QGIS?
However, it seems to produce UTM figures rather than degrees—even if I change my Project Properties from NAD 83 16N to WGS 84 for the exercise. Can you please explain where I’m amiss?

Thank you

TVZ


texts:
Is there any way to differentially correct a GPX file using free or open-source tools? Does a GPX file contain the necessary information for differential correction (I can put the timestamp in it accurate to 1s)?

As an aside, do I need to disable WAAS on my receiver for this to work?


texts:
I'm trying to learn how to work a bit with Qgis for a part of a university work, and I'm facing some difficulties with an old map, and setting it up with the right coordinates.

I have some old maps from a part of Portugal, and when I added them in Qgis all together, they fitted well and were all well ordered, so I started to mark some spots on them.

After, when I tried to add a Openmap over it, I saw that they didn't match, and the old maps were over Africa (and maybe with the wrong scale).

Is there any way of transposing the old maps and the marked points (a .shp file) to the right coordinates? I've spent quite some time marking all those spots, and I needed to mark them over the old maps, and have them fit in the right spot of the OpenStreetMap.

I've tried multiple CRS configs, moving them... But I can't seem to figure it out how to do this. Anyone can give me some lights on how to set the old maps and the .shp file on the right spot?

I can provide the files, and all the info needed


texts: I have successfully built my multi tab gui in Qt Creator but I don't know how do I implement it in the code. I already have my plugin files. I generated them using Plugin Builder. So how do I implement it from there?


texts: Creating border at a particular distance from country border

texts: I am new to GIS platform. I am using QGIS. I want to create lines inside map similar to border of country at a particular distance (say 100 KM). My data is in geographic coordinate system.


texts: Using ArcPy to add all unique values from TIFs

texts:
I am looking for a way to use ArcPy to add all unique values to a layer.

My tool converts ASC files to TIF files and then adds the layers to my mxd.

All works well: I have created layer files to store color attributes for each layer.

The problem is that some layers do not always have the same amount of unique values.

To explain:
if I convert a project A, it could be that layer BEST_SERV_MO99.TIF has only 3 unique values stored. If I convert an other project B, that layer BEST_SERV_MO99.TIF has now 12 unique values.

Is there a possibility to add for this layer the unique values by ArcPy?
If I use the layer file (.lyr) it will show for Project B only 3 colors, even if that layer would have 8 colors...

Here is a small snippet out of my converter:

filename_asc = i
tif_file = filename_asc.replace(".ASC", ".TIF")
tif2_file = filename_asc.replace(".ASC", "")
lyr_file = filename_asc.replace(".ASC", ".lyr")

outRaster = os.path.join(path_now, tif_file)
symbologyLayer = os.path.join(self.lyr_dir, lyr_file)

addLayer = arcpy.mapping.Layer(outRaster)
arcpy.mapping.AddLayerToGroup(self.df, targetGroupLayer, addLayer, "BOTTOM")
layers = arcpy.mapping.ListLayers(self.mxd)

self.pyramid(os.path.join(path_now, tif_file))
for layer in layers:
    if layer.name == tif_file:
        if os.path.isfile(symbologyLayer) == True:
            arcpy.ApplySymbologyFromLayer_management(layer,symbologyLayer)
        else:
            self.err("Farbschema konnte nicht angewendet werden: " + symbologyLayer.encode('utf-8'))

        layer.name = tif2_file

edit:
if I try to see the symbologyType, it tells me always: OTHER

self.log("layer.symbologyType: " + layer.symbologyType) # output --> layer.symbologyType: OTHER


texts: What is the best approach to using IQueryFilter.WhereClause in ArcObjects?

texts:
What is your preferred / best approach to using IQueryFilter.WhereClause statements in ArcObjects (e.g. when selecting features by attributes)?

You need to check the source type (shapefile / personal / file geodatabase etc.), then, check field types and so on, before you can actually construct the WhereClause statement and find the proper SQL syntax. I find it rather easy to get lost in all those apostrophes, quotation marks etc, especially if multiple fields are used in the query.

What is your best way to easily construct the WhereClause statement that would follow the valid SQL syntax..?


texts: GeoJSON: PolyhedralSurface geometry type not supported

texts:
I am trying to make a 3D web application which should be able to render 3D buildings (I will be using Cesium, but that is not that important for now).

I basically have a database consisting of the 2D geometry of the buildings and  the height, amongst other things. In that way I should be able make simple 3D building blocks (LoD1). I am using PostGIS (PostGIS 2.2.0 with SFCGAL extension) to extrude (st_extrude) the 2D geometry to a 3D geometry. No problem so far...

The problem I am facing now is that st_extrude gives me a polyhedral surface. However, I would like to convert my 3D geometry to GeoJSON (which can easily be loaded into my Cesium viewer). St_AsGeoJSON gives me an error saying that 'PolyhedralSurface' geometry type is not supported.

Do you guys know if I can work around this problem? Are there any similar functions in PostGIS like st_extrude which will not leave me with a polyhedral surface?

I'm looking into using st_Dump to try and transform my polyhedral surface into a multipolygon, because I will be able to make a GeoJSON type with a multipolygon.

As far as I can tell by looking at the documentation, using st_dump will give me six "facets" (say my building is a standard cube). Does someone know if I will be able to leave out the side-facets of my cube and just keep the top and bottom? Will I be able to go to a multipolygon from there? Just an idea I had today...


texts: Writing a final output to geodatabase in Model Builder

texts:
I wanted to write a Model Builder but I am not sure how to do save the final output as a geodatabase.

I am working on to classify each imagery (about 83 of them) and then sent it to a polygon and then rename it to the same name from the original file to save to the geodatabase.

It is working but I could not get the final output to change the pathname to the geodatabase I wanted.


texts: Is there a name for the latitude/longitude pair 0,0?  On a graph, (0,0) is referred to as the "origin" but I'm not sure if that term also applies to the location 0,0 in cartography.


texts: Given a top left location and center location find bottom right location

texts:
If I have a rectangle, in which top left point corresponds to a specific location (lat1, lng1) and center corresponds to another location (lat2, lng2), how can I find to which location corresponds bottom right point?

I will try to illustrate now what I mean. I have a rectangle (by rectangle I mean an actual map area, a randomly zoomed-in part of Mercator projection) (say x in width and y in height) with L point in top left, M point in the center and R point in the bottom right:

L----------------------+
|                      |
|                      |
|                      |
|          M           |
|                      |
|                      |
|                      |
+----------------------R

L(0, 0) or L(-x/2, y/2) (whichever coordinate system you prefer) corresponds to (lat1, lng1) and M corresponds to (lat2, lng2) which are known. Is it possible to find (lat3, lng3) of R from the given data? It is safe to assume that lat1 > lat2 and lng1 < lng2 always.

I have tried adjusting this formula (which is basically the inverse of my task) to solve the problem, but could not come to solvable equation, and the resulting system did not look easily solvable.


texts:
I am trying to load a KML file to Google Maps using OpenLayers 2, below is the code:

var wms = new OpenLayers.Layer.WMS(
          "OpenLayers WMS",
          "http://vmap0.tiles.osgeo.org/wms/vmap0",
          { layers: 'basic' }
      );

var sundials = new OpenLayers.Layer.Vector("KML", {
    projection: this.application.map.map.displayProjection,
    strategies: [new OpenLayers.Strategy.Fixed()],
    protocol: new OpenLayers.Protocol.HTTP({
        //url: "http://dev.openlayers.org/examples/kml/lines.kml",
        //url: this.application.CommonLocation + "scripts/Testing.kml",
        //url: 'data:application/vnd.google-earth.kml+xml;base64,' + this.tmpKml,
        url: 'data:application/vnd.google-earth.kml+xml;base64,' + "PD94bWwgdmVyc2lvbj0iMS4wIiBlbmNvZGluZz0idXRmLTgiPz4NCjxrbWwgeG1sbnM9Imh0dHA6Ly93d3cub3Blbmdpcy5uZXQva21sLzIuMiI+DQogIDxEb2N1bWVudD4NCiAgICA8U3R5bGUgaWQ9ImJvdW5kYXJ5TGluZUNvbG9yMSI+DQogICAgICA8TGluZVN0eWxlPg0KICAgICAgICA8Y29sb3I+ZmYwMDAwMDA8L2NvbG9yPg0KICAgICAgICA8d2lkdGg+Mzwvd2lkdGg+DQogICAgICA8L0xpbmVTdHlsZT4NCiAgICAgIDxQb2x5U3R5bGU+DQogICAgICAgIDxmaWxsPmZhbHNlPC9maWxsPg0KICAgICAgPC9Qb2x5U3R5bGU+DQogICAgPC9TdHlsZT4NCiAgICA8Rm9sZGVyPg0KICAgICAgPG5hbWU+U295YmVhbiBZaWVsZDwvbmFtZT4NCiAgICAgIDxQbGFjZW1hcms+DQogICAgICAgIDxuYW1lPlNveWJlYW4gWWllbGQgKDEpPC9uYW1lPg0KICAgICAgICA8c3R5bGVVcmw+I2JvdW5kYXJ5TGluZUNvbG9yMTwvc3R5bGVVcmw+DQogICAgICAgIDxQb2x5Z29uPg0KICAgICAgICAgIDxvdXRlckJvdW5kYXJ5SXM+DQogICAgICAgICAgICA8TGluZWFyUmluZz4NCiAgICAgICAgICAgICAgPGNvb3JkaW5hdGVzPi05MS4zMzQyNjYzNzU5MTU0LDMyLjUwMzQ3NjY1MDA5MjcNCi05MS4zMzQyNjY3MDg0ODMxLDMyLjUwMzQ1ODYxNTU1NDgNCi05MS4zMzQyOTg2MzI1MTcxLDMyLjUwMjg4MTc5MjEyMzINCi05MS4zMzQ3NTI4NjQzMDQ2LDMyLjUwMTMzNjQ1NjYwNTgNCi05MS4zMzUzNDkwODMzOTI2LDMyLjUwMTMyNjMwODc5NTMNCi05MS4zMzcxNzk5NzIwMzEzLDMyLjUwMTMxNDQ0NTU3Mw0KLTkxLjMzODExNjY5ODMyMSwzMi41MDEzMDg3ODUxMTIyDQotOTEuMzM4NjQ4NzQwODEwMiwzMi41MDEzMTU4MTI3NTQNCi05MS4zMzg2NjYzNzQyMjEyLDMyLjUwMTUxNDQ3NDIzODgNCi05MS4zMzg3MjYyOTUxNTEyLDMyLjUwMjg4NjIyNjcyMTYNCi05MS4zMzg3NDA2MTIxNzY2LDMyLjUwMzI2NTIzMzk2MTQNCi05MS4zMzg3ODA1Nzg0MjYyLDMyLjUwNDU2NDU2Njc1MTENCi05MS4zMzg5MTMxMzgwMDk2LDMyLjUwODkzMTc0MzUwMDQNCi05MS4zMzg5MjE4MTY0Nzk2LDMyLjUwOTYxNzMzODAwNDgNCi05MS4zMzg5ODU5MzIwNTc2LDMyLjUxNTM5MDY0NDM1MTgNCi05MS4zMzg5ODMyNzc4MjU3LDMyLjUxNTUzNDkyMDcyNTkNCi05MS4zMzg4MTEwMDY0ODk5LDMyLjUxNTY0MDg3ODk3NjENCi05MS4zMzc4NzQxMzEyNTQsMzIuNTE1NjQ2NTM5NzQwNg0KLTkxLjMzNTkzNjUyNjIyMzQsMzIuNTE1NjU2OTk1NDQyMg0KLTkxLjMzNDk1NzA4MTUyOCwzMi41MTU2NjIwNzExOTI1DQotOTEuMzM0NjU5MDkwOTUxMiwzMi41MTU2NTgxMjU3MDgNCi05MS4zMzQ1OTYyMzM4MDc4LDMyLjUxNTYwMzE3NjYzNzINCi05MS4zMzQ1Nzk5Mzg2NjkxLDMyLjUxNTMzMjM3NzE4NDgNCi05MS4zMzQ1MTA0MzkzMDU4LDMyLjUxMzMyOTEzNzM4MDMNCi05MS4zMzQ0OTY4MDU2NzQ1LDMyLjUxMjkxNDA2MTc0MzINCi05MS4zMzQyNjYzNzU5MTU0LDMyLjUwMzQ3NjY1MDA5Mjc8L2Nvb3JkaW5hdGVzPg0KICAgICAgICAgICAgPC9MaW5lYXJSaW5nPg0KICAgICAgICAgIDwvb3V0ZXJCb3VuZGFyeUlzPg0KICAgICAgICA8L1BvbHlnb24+DQogICAgICA8L1BsYWNlbWFyaz4NCiAgICA8L0ZvbGRlcj4NCiAgPC9Eb2N1bWVudD4NCjwva21sPg==",
        format: new OpenLayers.Format.KML({
            extractStyles: true,
            extractAttributes: true
        })
    })
});

However, Chrome throw the following error:

XMLHttpRequest cannot load data:application/vnd.google-earth.kml+xml;base64,PD94bWwgdmVyc2lvbj0iMS4wIi…AgICA8L1BsYWNlbWFyaz4NCiAgICA8L0ZvbGRlcj4NCiAgPC9Eb2N1bWVudD4NCjwva21sPg==. Request header field X-Requested-With is not allowed by Access-Control-Allow-Headers.

I have followed the step in:

https://stackoverflow.com/questions/12409600/error-request-header-field-content-type-is-not-allowed-by-access-control-allow

But still unable to solve it. Can someone please help?


texts: Multiple centroids of complex outline polygons by recursive buffering

texts:
How can I algorithmically determine multiple centroids for overlapping (dissolved), very irregular regions representing canopy outlines?

This is effectively a forestry issue. The context is that I have canopy outlines I generated with concave hull from a raster. What I want to do now is get centers for each 'bulb' of the complex outline (imagine many overlapping circles).

My thought is to use a negative buffer recursively until I have many small polygons that I can find the centroids of. The challenge is that each time I buffer I will get polygons which should be excluded because they are now independent, small, and should not be buffered again. This needs to be checked on every buffer run.

The raster in the second image is from a separate workflow -- an attempt to extract a point cloud based on coloring the points with NDVI.

I am working in QGIS 3.10, have very limited ability with Python, strong experience in ArcMap and decent experience with ModelBuilder but have never used the modeler in QGIS.


texts: ArcGIS Model isn't accepting my SQL Query

texts:
I'm running a model that breaks a dataset into several dbf tables based on a query. However, only half of the query is being honored.

The PopDenZone field is not being accepted and I am getting other values in the output, but not all of the total, maybe 50-100 records are being thrown into the output from other Density Zones.

Anybody know what could be causing this behavior? The portion of the query where I need specific Priority Descriptions is working correctly.

Using ArcGIS 10.1


texts: Is there a non-ArcPy way of joining a table to an existing raster attribute table using Python?

texts: I have an existing raster with attribute table and would like to join a csv to it. Is there an open source way of doing so using Python libraries like rasterio, pygeoprocessing, etc.?


texts: QGIS Grass-plugin output does not match the CRS of input files nor project CRS

texts:
My project CRS and raster layer are EPSG:31467. However, when I use for example r.to.vect to create a polygon from the raster, the created vector layer has an unknown CRS and is a little bit off (does not match the source raster). When I manually set the CRS to EPSG:31467 the vector layer fits the source raster. But, when in a next step I create a line vector from the polygon using v.to.lines it happens again - although the input vector is EPSG:31467 the resulting line vector layer has an unknown CRS and is off.

Where is the problem? Can it be that the original raster file is actually not EPSG:31467? Because when I load it, the CRS is also unknown at first. I have just been told that it should be EPSG:31467 and therefore set it manually.

I am running QGIS 3.12.0.


texts: find out those points which are intersecting each other using spatial filter in C# arcobjects

texts:
Possible Duplicate:
  Determine overlapping points using a spatial filter in ArcObject

I have a point layer having thousands of point, now what i have to do using Spatial filter i have to find out those points which are intersecting each other. And then objectid of intersecting point should be display in a list box.
in C# Arc Objects


texts: Writing value from attribute table to text element using ArcPy?

texts:
Is there a way, using ArcPy, to read a field from an attribute table and write the contents of that field into a text element?
I have a point feature class with a field name called "Cities". I will make separate MXDs for different subsections of cities but I would like the contents of the cities field to write into a text element.

Using the following blog, I can write the values to a notepad. But rather than a notepad, can I write them to a text element in the MXD?


texts:
I'm looking for sources containing data on measured river flow speeds, any / all locations worldwide, in any format or language.

Googling didn't help much.


texts:
I am using gdal_rasterize to create a raster of a huge file containing polygons. I managed to generate the raster file with the right color and with the minimum size to tile it after with the gdal2tile.

gdal_rasterize -of GTiff -co COMPRESS=DEFLATE
-b 1 -b 2 -b 3 -burn 136-burn 46 -burn 17
-co PREDICTOR=2
-co ZLEVEL=9
-ot Byte
-a_nodata 255
-a Niveau
-ts 500000 500000
-l layer W:/layer.shp
W:/layer.tif

Tile will be used to overlay Google Maps so the polygon must be transparent.

Has anybody experimented with this usage of gdal_rasterize and succeeded in setting the opacity of the output raster?


texts: Does anyone know an equivalent package for the velox package for fast zonal statistics in R? This package was by far the best performing tool when working with large raster objects but has been recently removed from CRAN...


texts: Ignoring certain data values in graduated maps

texts:
I'm working and creating a large amount of maps using the latest version of Qgis with vector files and corresponding data sheets. Usually, there will be missing values which is usually easily avoided by assigning these data values to "NULL" and adding an additional layer below the graduated one with a constant color, however when there are "several missing" values i run out of ideas of how to easily handle these.

If say i have two missing values, one with "NULL" and one with "-999", the "-999" values will mess with Qgis' automatic interval determination tools since it regards the -999 value as relevant - Is there any way to tell Qgis to ignore this value?

It seems as though the Feature Subset feature it is greyed out in my vector file - even though the vector file is in edit mode. This seems to be connected with the fact that my vector file is joined with an external data file - Is it really not possible to filter joined data?

From other questions here at stackexchange i learned that the query builder sadly doesn't work with joined layers - the only solution was to join layers, save as a new layer and apply the subset features here.


texts: What view of sky is required for RTK GPS

texts: When using RTK or differential GPS technologies, considering the Base has a clear unobstructed view, what view of the sky is required to be able to get an accurate fix with the rover when is moving (boat on water)?


texts: Geometry workflow from Shapely to GeoJSON

texts:
I'm undertaking a fairly simple workflow all within a python script with the goal of passing urlencoded GeoJSON in a URL to the Mapbox static API:

create a point from DD string in EPSG:4326:
convert to UTM for easier units and buffer the point generating a polygon
transform back to WGS84 for the API
generate valid, URL encoded GeoJSON string.

My hangup is the last part. I've used Shapely and pyproj to conduct 1-3 but using:

json.dumps(mapping(<shapely.geometry.polygon.Polygon at 0x106a22690>)) (and several permutations of this) and encoding it using urlencode or url_quote generates an encoded URL with JSON that returns 'Invalid Geojson' from Mapbox.

I am optimistic that the geojson library may solve this, but I can't seem to find a simple way to convert a shapely Polygon object to a GeoJSON polygon object.


texts: Filling attribute table automatically but with specific values in QGIS 3.10

texts:
I'm currently trying to built a database for a cocoa plantation. I would like to know how to fill my database easily with the field calculator. For example I would like to have a specific value for my field "Varieté", but not for all the column.

Is it possible to put the value "Upala Silva" only for the "ID" 1 to 70?


texts: How to know if two projections are the same?

texts:
I am using R but also tried with QGIS, visually the two maps look the same but the projection info it's not. QGIS provides the same info.

first map

+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

second map

 +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84
    +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

First I would like to know what these abbreviations mean, some are rather intuitive but not all of them. lat_ts? no_def? towgs?

How do I check with R that these two projections are the same?


texts: How to calculate tortuosity of a coastline?

texts:
I would like to calculate the tortuosity (or fractal dimension) of a coastline within a 1km buffer surrounding a point on that coastline. I have hundreds of points and therefore need a tool or scripts that will do the calculation on a mass scale.
Can someone help me please?

I am using ArcMap10


texts: I have a DEM of a city with numerous hills and a shapefile of the city streets. I want to find the five paths from a low point A to a high point B on the street network that involve the least elevation gain, that is, that have the fewest sections of downhill. Is it possible using QGIS/GRASS?


texts: Does WFS support what OGR calls "OLCFastSetNextByIndex"?

texts:
Background

I'm working on an application that has multiple threads getting features from the same OGR layer at the same time, but with different spatial filters. This is problematic because it causes collisions. To solve this I have the worker threads ask a different thread to get the features for them and place them in a queue for the workers. This introduces a different problem.

Since I have no control over how many features that might exist on a layer I have to put a limit on how many features to get at a time so as to not let the memory usage run away uncontrollably (I've currently set this limit to 50 per worker). That means that for every worker the feature getting thread has to know how many features it has received and start from there. In other words, it has to execute this code before it starts getting the next 50 features from the layer:

layer.SetAttributeFilter(aFilter);
layer.SetSpatialFilterRect(anEnvelope.MinX, anEnvelope.MinY, anEnvelope.MaxX, anEnvelope.MaxY);
layer.SetNextByIndex(index);

That layer.SetNextByIndex is a massive time sink. OGR does not have an optimized version of that function implemented in its WFS layer class so what it does is reset the reading to zero and do GetFeature() index number of times. layer.TestCapability(Ogr.OLCFastSetNextByIndex); always returns false, as can be seen in ogrwfslayer.cpp.

Question

What I'm wondering is whether OGR doesn't have that function implemented due to a hiccup or laziness or if WFS simply doesn't have that capability.

I'm using GeoServer for the WFS. Is it possible for me to implement the OLCFastSetNextByIndex functionality on my own?

Sources:
http://www.gdal.org/classOGRLayer.html#ac0b0869ff6fcadda35cbac1f3edf8eae
https://github.com/OSGeo/gdal/blob/master/gdal/ogr/ogrsf_frmts/wfs/ogrwfslayer.cpp
https://github.com/OSGeo/gdal/blob/master/gdal/ogr/ogrsf_frmts/generic/ogrlayer.cpp


texts: Using Sort_Management for multiple tables

texts:
I like to convert tables (100+) living in access database to a geodatabase tables, sorting them and deleting unnecessary fields in each table at the same time. I figured I use the Sort_management method to convert and sort the tables followed by DeleteField_management method to delete the fields. Here is the code that I have written,

import arcpy

from arcpy import env

env.workspace = "E:\Database Project\ACS Estimate 2011\ACS 2011 Tables.mdb"

tableList = arcpy.ListTables()

for table in tableList:

    out_dataset = "E:\Database Project\ACS Estimate 2011\Data\To Database\Seq Tables.gdb"

    #Process: Sort
    arcpy.Sort_management(table, out_dataset, "GEOID ASCENDING", "UR")

    #Process: Delete Field
    arcpy.DeleteField_management(out_dataset, "FILEID;FILETYPE;STUSAB;CHARITER")

Here is the error message that I am getting,


texts: Debug target is missing - how to set in VS 2010 Express?

texts: I updated ArcGIS from 10.1 to 10.2 and recompiled all Addins with ArcObjects 10.2.  I am having an issue with 1 addin that accesses a custom tool in a toolbox and want to debug the project.  However, when i hit run I get a "...the debug target '' is missing" erro message.  The path is still set to the Desktop10.1 target - HOW DO I RESET IT to Desktop10.2 in VS 2010 EXPRESS. I only found instruction on how this can be dine in the full version if VS but this option is not avaialble in the express version.


texts: OWSLib getmap - how to encode layernames with non-ascii charaters?

texts:
Whenever i think I understand encodings, I don't.. How do I correctly send a layername to the getmap function?

running python 2.7

from owslib.wms import WebMapService
wms_url = "http://ext-geodata.lansstyrelsen.se/arcgis/services/WMS/LST_WMS_riksintressen_3_5_st/mapserver/wmsserver?"

wms = WebMapService(wms_url)
my_layer = list(wms.contents)[1]

#output = u'ST Riksintresse K\xe4rnomr\xe5de  MB3kap5'

img = wms.getmap(layers=[my_layer],
                 styles=['default'],
                 srs='EPSG:3006',
                 bbox=(7000000,577000,7113000,690000),
                 size=(500, 500),
                 format='image/png',
                 transparent=True
                 )

out = open('testarea.png', 'wb')
out.write(img.read())
out.close()

When I run this, I get a UnicodeEncodeError. I understand I need to encode the layername correctly, but I fail to understand how.

l = my_layer.encode('utf-8')
img = wms.getmap(layers=[l], #etc etc...

yields

ServiceException: Parameter 'layers' contains unacceptable layer names.

How do I parse values from the contents list to layername correctly?


texts: I use Mean coordinates tool for a layer where the UID is a string and weights are real numbers. I find it very slow (literally hanging QGIS for few minutes) when processing datasets of more than few hundred points. Is there a fix or a workaround?


texts: Show GIS data on website using automated qgis2web process

texts:
I managed to export my GIS map to my website through Qgis2Web. My map shows a daily view on a specific infestation of parasites in the poultry industry and I want this updated once a day, say at midnight CET.
Is there a way to do this?


texts: Automatically changing the angle of a symbol in QGIS?

texts:
I am working with an iron age cemetery in QGIS. I marked every grave in a point layer with several columns; amongst others, burial rite (inhumation, cremation etc.). I inserted a second column for the orientation of the inhumations. I categorized the vector symbols according to the burial rite column and created one simple vector symbol for the inhumations, now I would like it to rotate automatically according to the orientation column.

Is that possible using one layer, or do I have to put all the inhumations in a second layer.


texts: Coloring polygon based on date using TileMill/CartoDB?

texts:
How would you color a polygon/line/dot based on a date? The colors would be a pre-set list. I need the area to update its color based on the current (or given) date. I would be using tilemill or cartodb.

Example:

Polygon A:
Color would be white after Jan 1,
Color would be grey after Feb 15
Color would be brown after March 1, etc.

Line B:
Color would be blue after March 15,
Color would be white after Dec 1, etc.