texts: Geometry Generator - Straight lines not curved

texts:
I have QGIS 3.8.2 installed and started to look at the Geometry Generator.  I have a layer of airline tracks and trying to create a curved line linking start and end points.

I did a lot of searching and found the followng code but when I apply this I get straight lines not curved lines like I have seen elsewhere.

make_line(
    start_point($geometry),
    centroid(
        offset_curve(
            $geometry,
            length($geometry)/-5.0
        )
    ),
    end_point($geometry)
)

Hopefully its something simple that I am doing wrong.


texts: Using OpenLayers example in OpenStreetMap Wiki?

texts:
I'm trying to work my way through this wiki page to add OSM into a OpenLayers map I am working on, however I cannot even get the first example to work.

I have downloaded OpenLayers V3, copied the code into a text file (saved as html), then copied "os.js" into the same folder, not "openlayers.js" as it says in the example, this doesn't appear to exist so i think i copied the correct one. I'm also not sure what it means as the "theme directory" as well.

Anyone managed to get this to work, if so then how? And is there any other ways this could be done?


texts: Freehand editing without holding down mouse button in QGIS?

texts:
I am trying to delineate wetlands using satellite imagery in QGIS 2.10.1-Pisa.

I'm using the Freehand Editing plugin.

Is there a way to draw polygons without having to hold down the mouse button, like you can with the freehand tool in ArcMap?


texts: Different ST_AZIMUTH results for ST_PROJECT points along the same original bearing

texts:
I'm trying to develop some way of determining line of sight (LOS) between two geographical points about 1km apart or less using ST_3DINTERSECTS or something similar.

BASIC PROBLEM:

Along the way, I've found that if I calculate the ST_AZIMUTH between two points and then use that bearing in a series of ST_PROJECT statements to create additional points along the line connecting the two original points, I get a close but slightly different ST_AZIMUTH when calculating from the original start point out to each of the new points along the line.

I believe this slight difference in bearings between the lines is causing my ST_3DINTERSECTS function call to return false for all scenarios, even when terrain elevation should be blocking LOS between the two original points.

DETAILS:

To begin, I have a start point (SP) and an end point (EP) both represented as geographies.  I calculate the original bearing: (OB) = ST_AZIMUTH(SP, EP).  Then I use my STEP size (STEP = 30 meters), SP and OB to create points along the line connecting SP and EP in order to pull elevation samples along the way: ST_PROJECT(SP, STEP, OB).  In this manner, I now have roughly 30 geography points with elevations along the line connecting SP and EP.

In all cases, elevation data is retrieved using my raster data in conjunction with these geography points casted as geometries e.g.: ELEV = ST_VALUE(myraster, SP::geometry).  I also extract the X and Y coordinates from the geography points casted as geometries e.g.: X = ST_X(SP::geometry) and Y = ST_Y(SP::geometry).  Using this information, I create a series of 3-dimensional points e.g.: ST_MAKEPOINT(X, Y, ELEV) in order to represent the original SP and EP as well as all sample points along the line connecting SP and EP.

Using these 3D points, I create two lines e.g.: ST_MAKELINE(my3dpoints).  One is a simple line connecting just the SP and the EP in 3D.  The other is a segmented line comprised of all the sample points taken every 30 meters along the straight line connecting SP and EP.  When I use ST_3DINTERSECTS on both of these lines, it always returns a false result, even if I know the terrain should be obstructing LOS.

QUESTION:

Can anyone verify if the slight bearing difference is what's causing my lines not to intersect when they should?  Is there a better way to do this?

SETUP:

Example SP = POINT(-74.714919, 41.282384), EP = POINT(-74.7147106468, 41.2913868143398)
SP is always 3 meters above the actual terrain elevation.  EP is always 2 meters above the actual terrain elevation.  These two points should not have LOS due to terrain.

My elevation data source is: http://www.viewfinderpanoramas.org/dem3/K18.zip loaded using raster2pgsql with SRID 4326.

My postgis_full_version() is: POSTGIS="2.1.1 r12113" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.0, released 2013/04/24" LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER


texts: Getting ERROR 030098 from Generate Service Areas tool of ArcGIS Network Analyst?

texts:
I'm attempting to run the 'generate service areas' tool. I've just created a network following the instructions here using the data from the Utah Automatic Geographic Reference Center here. I wish to create service areas for a series of points within Utah. There are many points of which i present the lat-lon for 3:

40.5361, -112.025
40.6874, -111.915
40.7095, -112.022

Procedure
Download it from AGRC, deleted the 'test', deleted the 'network', and 'network_junctions', leaving on the 'roads' file within the 'network' feature dataset. To reduce the time required to build the network, I selected the road network by location, buffering my  apartment locations by 2 miles, and exported that as a new file in the existing feature dataset as 'roads_2miles'. I then followed the procedures for creating a network according to the tutorial here, for ArcGIS 10.3. I omitted steps 39, 40, and 41 as irrelevant to this task. Building the network took about 2 minutes.

I then attempt to use the 'generate service areas' tool, using my apartment location as the 'facilities' and my 'Network_2miles' network dataset as the Network Dataset. I receive the following error:

  ERROR 030098 The network dataset must have at least one time-based and
  one distance-based cost attribute.

What am I doing wrong here?

Potentially Relevant:
http://desktop.arcgis.com/en/arcmap/10.3/tools/network-analyst-toolbox/make-service-area-layer.htm
http://desktop.arcgis.com/en/arcmap/10.3/tools/network-analyst-toolbox/make-service-area-layer.htm


texts: Adding more basemaps to Add Basemap window of ArcMap?

texts:
How can I add more types of basemaps to the existing Add Basemap window?

I am using ArcMap 10.4 on Windows 10


texts:
What am I doing wrong if my raster image becomes completely black after polygonising the tif file.
I tried converting the following raster image into vector

And got this as output


texts: OpenLayers 5 Vector Tiles Styling and ESPG:4326 projection issues

texts:
I have a OL5 map using a projection of EPSG:4326. I am adding a Vector Tile layer with a projection of ESPG:4326.

I am using Open Layers Mapbox Style (https://github.com/boundlessgeo/ol-mapbox-style) to apply a style to my Vector Tile I am seeing weird issues with how the styling is being rendered. I am a rookie with this stuff but I am assuming there is a disconnect with the Vector Tile Styling data vs the projection my map supports.

Any idea on why the below code does not seem to style the Vector Tile correctly?

Is this an issue with the plugin or an issue with the way the layer is configured?

Is it because "if the vector tile layer was built on a map with a projection of EPSG:3857 it can not be used/styled on a map with a different projection"?

// Create vector tile layer
var layer = new ol.layer.VectorTile({
    source: new ol.source.VectorTile({
    format: new ol.format.MVT({ featureClass: ol.render.Feature }),
    url: "https://basemaps.arcgis.com/v1/arcgis/rest/services/World_Basemap/VectorTileServer/tile/{z}/{y}/{x}.pbf",
   projection: ol.proj.get("EPSG:4326")
   })
});
// Create map
var map = new ol.Map({
    target: 'map',
    view: new ol.View({
       projection: 'EPSG:4326',
       center: [0, 0],
       zoom: 2
    }),
    layers: [layer]
});

// Get style
$.get("https://www.arcgis.com/sharing/rest/content/items/291da5eab3a0412593b66d384379f89f/resources/styles/root.json")[![enter image description here][1]][1]
.then(function(res) {
    olms.applyStyle(layer, res, "esri");
 });


texts: Rasterio showing single band images wrong

texts:
I'm using rasterio in Python at the moment, and I have no problem when I'm dealing with RGB satellite images. However, if I try using singleband Sentinel 2 images, or a singleband Sentinel 1 image, something happens. As this figure here shows, you can see that the image on the left, is suddenly green. Whilst the exact same image on the right in QGIS, looks at it should. Any suggestions?

I use the following code:

dataset = rasterio.open(project_dir + '/data/satellite_data/S1/merged/S1HH.tiff')
img=dataset.read()
fig=show(img)


texts:
I'm trying to get the cursor to iterate through rows and fields looking for blank or null values and change them to 999999 but it's not changing anything.

fc = "C:\Users\\bbrock\Documents\ArcGIS\Ports.shp"

# Create a search cursor
#
rows = arcpy.SearchCursor(fc)

# Create a list of string fields
fields = arcpy.ListFields(fc, "", "String")

for row in rows:
    for field in fields:
        if field.type != "Geometry":
            if row.getValue(field.name) == '':
                row.setValue(field.name, '999999')
                print "%s: Value = %s" % (field.name, row.getValue(field.name))

            if row.isNull(field.name):
                row.setValue(field.name, '999999')
                print "%s: Value = %s" % (field.name, row.getValue(field.name))


texts:
Does anybody know, how I can update the QGIS Server without changing all files.
I have a Version from January and want to update it.

Which files I have to change ?


texts: geometry point and varchar multicolumn index in PostgreSQL

texts:
Is it possible to create a multicolumn index on a geometry point column and a varchar column? For example, say I have a geometry point column called "location_point" and a varchar column called "zip". I have tried this...

CREATE INDEX contacts_80_zip_location_index ON contacts_80 USING GIST (location_point, zip);

This did not work. I received this error...

 ERROR:  data type character varying has no default operator class for access method "gist"
HINT:  You must specify an operator class for the index or define a default operator class for the data type.
SQL state: 42704


texts: Using TIGER/Line shape files, how do I get address coordinates?

texts:
I've downloaded addressfeature shapefiles from TIGER/Line on the Census Bureau website. I converted one file to CSV using ogr2ogr. I get the following sample output:

WKT,TLID,TFIDL,TFIDR,ARIDL,ARIDR,LINEARID,FULLNAME,LFROMHN,LTOHN,RFROMHN,RTOHN,ZIPL,ZIPR,EDGE_MTFCC,ROAD_MTFCC,PARITYL,PARITYR,PLUS4L,PLUS4R,LFROMTYP,LTOTYP,RFROMTYP,RTOTYP,OFFSETL,OFFSETR
"LINESTRING (-98.730002999999982 37.659596999999991,-98.730001 37.659660999999986,-98.729996999999983 37.659855999999991,-98.729995999999986 37.65992099999999)",161042451,202265803,202266201,400110711124,,110200839772,Stout St,1409,1499,,,67124,,S1400,S1400,O,,,,I,,,,N,N

I understand the LINESTRING (etc etc etc) portion are the coordinates for the address range (in this case, 1409-1499 Stout St).

How can I get the exact coordinates for an address, such as 1431 Stout St?


texts: Adding PostGIS raster to QGIS via database manager

texts:
When I attempt to add a PostGIS raster to QGIS I get a GDAL error saying I need to set PGHOST. I have set PGHOST with

export PGHOST=localhost

but still get the error. How can I resolve this problem?

Using an Ubuntu 14.10 machine with Postgresql 9.4, PostGIS 2.1.3 and QGIS 2.6.


texts: Network distance computations using GRASS

texts:
I'm an ArcGIS user who's interested in learning to use GRASS. Can anyone point out resources that will help me in learning how to obtain network distances between (1) residential locations (points) and (2) service locations (points also)? While I do wish to obtain the shortest paths, I need multiple paths between each residential location and a given number of services (say, the shortest paths to the 5 closest stores for each residential locations).

Help and suggestions will be much appreciated.

MP


texts:
I am creating a GIS application using OpenJUMP. I would like to show Google Maps as the base layer on top of which other spatial data is being rendered using PostGIS database.

Kindly let me know how can I go about integrating Google Maps in OpenJUMP application.


texts: Is it possible to make satellite imagery viewable through a marker?

texts:
I am working on the project using Mapbox.js and jQuery and the client has asked for something that I'm not sure is possible.

A screenshot of what I currently have is below. The client wants the circular markers to appear on a vector base map (as they currently do), but inside of the vector markers, they want satellite imagery to be shown (like a mask).

Is this even possible? If so, how would I go about doing this? I tried loading both tile sets at the same time, but I haven't been able to get that working. Currently the circles are being generated as markers, but based on Chris W.'s comment, it seems like I might be able to make these markers into masks?


texts: Changing Color Ramp in ArcGIS Online Heatmap?

texts:
I have a series of points that are being visualized as a heatmap in ArcGIS Online. The issue is that my team does not particularly like any of the 4 default color ramps offered by ESRI in the "symbols" tab (see the image).

Is it possible to alter these color ramps without requiring additional resources--ArcGIS for Server/ etc? (I'm aware of how to do it in ESRI leaflet). It just seems strange that there would be such a limited palette in AGOL.


texts: vb.net arcobjects enable customized context menu item only if a feature is selected

texts:
The solution may be very simple, but I haven't made it work.

I want to enable a data view context menu item only if a feature (in polyline layer) is selected, all other places the item should in grey. May I say what I need to do is check whether a feature is selected? I searched for this but don't find a proper solution. Or any other concerns?

I am using ArcObjects 10.2.2 and VB.NET.


texts: Filter features from a polygon shapefile in QGIS

texts:
I want to filter features from a polygon shapefile in QGIS using a list of numbers that occur only in one column of the attribute table.

For example:
My attribute table has 20 columns. One of them is named "telephones". I have a txt file that only contains one telephone number per row. How can I filter all features that assign to each number?


texts: animate polygons in series of kml/geojson/shp

texts: I have a series of KML files (though they could easily be converted to geoJSON, SHP, or anything else that ogr2ogr can handle), each file containing a set of polygons that correspond to the alpha shape (concave hull) of my movements for that day. I'd like to animate the polygons so that I can see how my personal alpha shape changes over time. Is there a particularly straightforward way to do this? I've looked at what Google Earth can offer but I can't help but think there's an easier way.


texts:
Does anyone know good resources to learn GeoMedia Professional?

I am looking for resources which can serve from a future teaching perspective, but other kinds of resources are also welcome.


texts: Understanding AutoCAD to ArcGIS conversion data limits?

texts:
I'm working to transition some of our data from CAD to GIS by bringing it in by adding the .dwg file to ArcMap itself and remapping using the DWG as a guide. I was wondering if there was a limit as to how many features ArcMap can handle? Is it determined by the software or hardware? I had initially tried bringing in about ~1400 features from AutoCAD, which caused crashes any time I tried to edit a feature. It seems like ~650 features worked okay.

Edit: I'm using Civil 3D 2016 and ArcMap 10.4.
Pretty basic features from Civil 3D, just points and lines that contain X,Y,Z.

I just added the some 13k features from CAD through on my personal computer with the GIS data I was using and it seemed to run just fine. Mine is spec'd a bit better than my work machine so I would assume it's a hardware issue.

Short of getting a new machine or doing some hardware upgrades, are there any settings on the software side that might ease the load?


texts: ArcMap Layer Time Export Using A Time Window

texts:
I currently have a layer with around 17,000 entries with time data formatted as YYYYMMDDHHMMSS as a string. ArcMap can correctly display the time and and I have applied an hour time window to show all of the available data points from a one hour window (for example: show all points in the range of 1:00 PM - 2:00 PM). ArcMap will correctly play this animation in the time slider by clicking play.

The issue arises when I go to export this as either a AVI or sequential images. Arc begins the export and ignores the time window and only shows 1:00 PM instead of the one hour window from 1-2. Because there is almost no data falling on the exact hour, the export fails to produce what just played in the time slider before exporting.

Screenshot:

http://imgur.com/a/Ay4EO


texts:
I am using a graduated legend, using the default format "%1 - %2". The numbers are formatted using decimal points, like "1.1 - 1.2".
How could I change this to using decimal commas, like "1,1 - 1,2"?
I tried the suggestion below, changing the locale to Dutch (interface language still American English. That does change the format of the values, but not the legend (see screenshot), which is what I need.


texts: Identifying coordinate system of Wetterzentrale maps?

texts:
I'm trying to georeference two maps from wetterzentrale in ArcGIS Desktop.

In what coordinate system are those maps created?


texts: Adding 'Add Clause' in ArcPy like with ModelBuilder?

texts: I am trying to set up an input parameter that is a drop down menu which is an input layer field. I can easily do that with ModelBuilder (Add Clause).


texts: Can't create a graduated symbology on a vector layer

texts:
I am trying to display a graduated symbology for a vector layer based on UK postcodes, but having major issues.
I have the postcode vector layer, and my data set is UK Postcode Lending data, where I've taken only the last column of data (Q2 2019) and saved this (along with the postcode sector as column 1) as a csv and uploaded to QGIS 3.12.3 Bucaresti. I'm running a Macbook Pro 2016, Mac OS Catalina 10.15.5. So far so good.
I've then joined my attribute only table of postcode lending data with the vector layer of UK postcodes, using postcode sector as the join. The join goes through, but the lending data gets saved as a Qstring text layer, not as an integer, no matter what I do. This is where the problem starts, as QGIS obviously can't make graduated symbology from a text attribute type.
I've tried to run Refactor Fields from the processing toolbox, but this doesn't change the attribute type for some reason, even though that's what it's supposed to do.
I've also tried reformatting and importing the CSV file so that QGIS imports it as an integer type off the bat - again, no dice.
These are the only two solutions I can find on the message boards, and I don't know why it's not working.


texts: How to find the Minimum value for the Dark Object Subtraction in ERDAS?

texts: I have a Landsat 7 Image with MTL file showing minimum value as 1 and the maximum value is 255. The water pixels show a range of 63 to 74. In the image, there are a few pixels that have a value of 1, but the histogram starts from around 63. I am trying to do Dark Object Subject using model maker in Erdas Imagine. What should I consider as a minimum value for DOS/haze correction?


texts: Manually moving layers in QGIS legend in print composer

texts: I have a legend in my QGIS print composer and would like to manually move a layer to position them where I want them to be. In the following example I would like to move the layer "area" down to the same postition as "Messtationen". Is there a way to do that? Ideally, a drag & drop function but other ways are also fine.


texts: Using keyboard shortcuts to invoke custom ArcGIS Base commands using ArcObjects and C#?

texts:
I developed multiple ArcGIS Commands added in a toolbar. I want to invoke these custom commands based on the keyboard keys pressed by the user.

How to implement this?


texts: Chart - trendlines - slope - NDVI - visibleInLegend: true - GEE

texts:
I am using the following code in GEE to calculate the trend line:

ndviChart.setOptions({
    title: 'NDVI (Sentinel 2)',
    vAxis: {title: 'NDVI'},
    trendlines: {type: 'exponential', 0: {color: 'CC0000'}, visibleInLegend: true},
    hAxis: {title: 'Data', format: 'MM-YY', gridlines: {count: 7}},
  });

The above code generates the trendline formula as a legend in the chart:

nd = 73.163e^-1.151E-4*system:time_start

where the slope is indicated in the caption as 73.163e ^ -1.151E-4

How do I extract the subtitle slope and store this value in a variable?


texts: How can I unselect a feature in Openlayers?

texts:
I am able to identify the selected features on my map using the following line:

vector_layer.selectedFeatures

I would like to be able to deselect this feature using a standalone function (not an event-handler). I've tried using the functions unselect() and unselectAll() unsuccessfully and can't find any examples of this approach.

I am aware that clicking on the selected object can be used to unselect. In this case I want a programmatic solution.


texts:
I have several tile WMS layers with about 5000 earthquakes and I want to query their features by clicking: example
That works fine for me, but I would like to query individual features, like only the date, the coordinates and the magnitude and not the whole attribute table.

So is there a way to query individual features of a tile WMS?

My recent research has shown that I should send the URL as JSON and than parse it.


texts:
Create Feature and store to different FGDB
I'm creating a feature and I want to store it to a specific geodatabase. At the moment I'm using feature.Store() but it's not storing it to where I want. How do I specify a location?

I want to store it to C:\College Stuff\Dissertation\2009 stuff\Geo\Geo\Full_3D.gdb

I've been told that if I am looking to store in a particular database I will need to open the feature class from that database.

I've looked at connecting to a geodatabase using this code

    Public Shared Function FileGdbWorkspaceFromPath(ByVal FGDB_WORKSPACE As String) As IWorkspace
    Dim factoryType As Type = Type.GetTypeFromProgID("esriDataSourcesGDB.FileGDBWorkspaceFactory")
    Dim workspaceFactory As IWorkspaceFactory = CType(Activator.CreateInstance(factoryType), IWorkspaceFactory)
    Return workspaceFactory.OpenFromFile(FGDB_WORKSPACE, 0)
End Function

This is my CreateFeature sub

 Public Sub CreateFeature(ByVal layername1 As String, ByVal x As Double, ByVal y As Double, ByVal FeatureClass1 As IFeatureClass)

    Dim cPoint As IPoint = New ESRI.ArcGIS.Geometry.Point
    ' To individually assign a value to the coordinates, use cPoint.X = x and cPoint.Y = y.
    cPoint.PutCoords(x, y)

    ' Ensure the feature class contains points.
    If FeatureClass1.ShapeType <> esriGeometryType.esriGeometryPoint Then
        Return
    End If

    ' Build the feature.
    Dim feature As IFeature = FeatureClass1.CreateFeature()
    feature.Shape = cPoint

    ' Commit the new feature to the geodatabase.
    feature.Store()

End Sub

To create the feature do I need to connect to the geodatabase in the CreateFeature sub? How do I do that? I'm modifying someone elses code and I am new to this so that's why I am very unsure

Here is the GetFeatureClass function, this where the CreateFeature gets the feature from. GetFeatureClass gets a feature class from the active map. The feature class is a single shapefile, I want it to get a feature class that is in the active map but in a geodatabase

Public Function GetFeatureClass(ByVal layername1 As String) As IFeatureClass
    Dim mxDocument As IMxDocument
    mxDocument = m_application.Document
    Dim map As IMap = mxDocument.FocusMap
    Dim check As Boolean = False

    Dim layerCounter As Integer
    Dim FeatureClass1 As IFeatureClass = Nothing
    Dim featureLayer As IFeatureLayer
    Dim currentLayer As ILayer
    For layerCounter = 0 To map.LayerCount - 1 Step layerCounter + 1
        currentLayer = map.Layer(layerCounter)

        If (TypeOf currentLayer Is IFeatureLayer) And (currentLayer.Name = layername1) Then

            featureLayer = currentLayer

            FeatureClass1 = featureLayer.FeatureClass

            check = True

            layerCounter = map.LayerCount
        End If
    Next

    If check = False Then
        MsgBox("cannot find InputStops layer")
        Return Nothing
    End If

    Return FeatureClass1
End Function


texts: Unable to list Sentinel l1c file using AWS CLI

texts:
Following Setting for AWS CLI using AWS configure

AWS Access Key ID : **********
AWS Secret Access Key : ********
Default region name : eu-central-1

which I get from IAM user
After setup I try to run following commands
aws s3 ls s3://sentinel-s2-l1c/ --request-payer requester

getting error as

An error occurred (AccessDenied) when calling the ListObjects operation: Access Denied

when I try to run
 aws s3 ls --request-payer requester s3://sentinel-s2-l1c/

getting same error
And when I try to copy the folder I get
aws s3 cp --request-payer requester s3://sentinel-s2-l1c/tiles/43/R/EQ/2017/7/productInfo.json pi.json

Error like this

Unknown options: --request-payer,pi.json

When I try to use wget
wget -O test.tar.gz "http://sentinel-s2-l1c.s3-website.eu-central-1.amazonaws.com/tiles/43/R/EQ/2017/7/7/0/B01.jp2?AWSAccessKeyId=********&Signature=**********"

I get error 403
How can I do this?


texts: How to create a MultiLineString feature with a postgis layer in QGIS

texts:
I have troubles handling with MultiLineStrings in QGIS.
I added the table where I defined a column as geometry(MultiLineString, 31468) as a postgis layer in QGIS. Now I want to create features with multiple linestrings.
When I finish to create the feature (with right click), then I have to insert an id (which is a value for another column in the table). So I have to create the next linestring feature and to assign the same id, to have a kind of a multiple linestring.
id  |  geom
1   |  "linestring1"
1   |  "linestring2"

I think this is not a proper way to work with multilinestrings!?
I want to have this:
  id  |  geom
  1   |  "linestring1" and "linestring2"

How can I create a MultiLineString in QGIS/postgis like this?


texts: Add z value from point layer to vertexes on polyline layer

texts:
I'm trying to get my z values from my 3d point layer over to my vertex' on a 3d polyline layer, so I can have lines with slopes in 3d.

I've first done a spatial location by intersect to filter out lines that aren't touching points. I've converted the points to 3d by attribute using their Z-values. I've done the same with the line but with a 0 as a value.

Now I would like for each end of the line to have that points Z-value, so I can get a sloped line in 3d. I've tried seeing if a tool like Add Geometry Attributes and feature to 3d by attributes works, but can't really figure out the proper proces.

I have access to ArcGIS Desktop Advanced 10.5.1, Spatial Analyst and 3D Analyst.

Edit:

After BERA's answer I've tried his code, modified to my workspace:

import arcpy

arcpy.env.workspace = r'D:\ArcGIS Geodatabase\Default.gdb' #Change
polylinefc = 'Vandledning_med_kotepunkt_Sp3' #Change to input line name, have to be in workspace
pointfc = 'Vandledning_med_kotepunkt_Sp1' #Change to input Point name, have to be in workspace

pointzdict = {k:v for k,v in arcpy.da.SearchCursor(pointfc,['SHAPE@XY','SHAPE@Z'])} #Create dictionary of points with XY as key and Z as values

with arcpy.da.UpdateCursor(polylinefc,'SHAPE@') as cursor: #Update Z of polyline vertices where a match is found in point dictionary
    for row in cursor:
        linelist = []
        for part in row[0]:
            partlist=[]
            for pnt in part:
                if (pnt.X,pnt.Y) in pointzdict:
                    partlist.append(arcpy.Point(pnt.X,pnt.Y,pointzdict[(pnt.X,pnt.Y)]))
                else:
                    partlist.append(pnt)
            linelist.append(partlist)
        row[0] = arcpy.Polyline(arcpy.Array(linelist), arcpy.Describe(polylinefc).spatialReference, True)
        cursor.updateRow(row)

However I get the following error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 14, in <module>
TypeError: 'geoprocessing point object' object is not iterable


texts: Split polygon using points and their directions with python

texts:
I`m new for the python but want to split Voronoi polygons using points and their directions as input.
Voronoi polygons were created by in-built tool: Vector>Geometry Tools>Voronoi Polygons from layer with these fields:

X;Y;ShapeID;LocationID;Azimuth
37.430027;55.518814;5;1;115
37.430027;55.518814;4;1;325
37.430027;55.518814;6;1;235
37.44712;55.535125;3;3;250
37.44712;55.535125;2;3;100
37.44712;55.535125;1;3;355
37.450652;55.520555;9;2;240
37.450652;55.520555;8;2;120

Could someone advise what are the steps to split each polygon into parts with python using
location coordinates and azimuths of shapes?
I plan those steps:

Merge two layers (Voronoi polygons + shapes).
For each shape draw the lines (start point is LocationID coordinates, end point is intersection with line with direction =
Azimuth +/- 360/(2 * shape_quantity_in_location)).
Split each polygon by lines calculated in step 2.

So end points must be calculated precisely to avoid errors.
The result should looks like this:

Black lines - Voronoi polygons, Violet lines - shapes azimuths, Green lines - lines from step 2. Black points - locations.
I`m stacked with this task.
I'm using QGIS 2.14.0-Essen.


texts: Huge x,y,z dataset into oracle, possible?

texts:
I have 200 .xyz files containing 1500000 points each. Is this possible to import this dataset to oracle?

It makes a total of 300.000.000 points!!!

Background is I want to improve a FME script that cuts out a desired region of elevation data. Today it uses a .shp grid to find the .xyz files intersected. This is a done with two linked fme workspaces, very complicated procedure. If I can get it into oracle, an oracle reader with bounds would do the same job more straightforward..

Any input regarding elevation data in databases or FME retrival of elevation data is much appreciated.

Best regards!


texts: pgsql2shp - Number of shapes does not match number of records

texts:
I have a MultiLineString postgis layer in a db.  The layer is created with the following...

-- Table: public.pipe_clip_102737

-- DROP TABLE public.pipe_clip_102737;

CREATE TABLE public.pipe_clip_102737
(
  geom geometry(MultiLineString,102737),
  gid integer NOT NULL DEFAULT nextval('pipe_clip_102737_gid_seq'::regclass),
  CONSTRAINT pipe_clip_102737_pkey PRIMARY KEY (gid)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE public.pipe_clip_102737
  OWNER TO jcfcqrso;

-- Index: public.pipe_clip_102737_geom_idx

-- DROP INDEX public.pipe_clip_102737_geom_idx;

CREATE INDEX pipe_clip_102737_geom_idx
  ON public.pipe_clip_102737
  USING gist
  (geom);

The layer displays correctly in both ArcGIS and QGIS when connecting directly to the db.

However will not display when I try to export to shp via cmd line or via the PostGIS DBF and Shapefile Loader utility.

pgsql2shp -f package/pipe_clip3.shp -h xxxxxxxxxxxxxxxxxxx -u xxxxx -P xxxxxxxx routing "SELECT gid, geom FROM pipe_clip_102737"

Both methods dump 1727 records as a shp.  But cannot be displayed. ArcGIS displays:

  "Error opening feature class. Number of shapes does not match number
  of records."

and QGIS seems to be aware of the extents of the shp but will not draw features and table is empty.

Why can I successfully display db table but not export a usable shp using either method?

EDIT: The issue seems to be related to the cloumn name gid. As mentioned by @tinlyx renaming the gid column while running the cmd seem to fix the problem.

pgsql2shp -f myshape.shp -h xx.xx.xx.xx -u myuser -P xxxxx mydb "SELECT gid as myid, geom FROM mytable;"


texts:
I am trying to do the rule based styling in openlayers 3, here is my code:

    CONFIG

I am trying trying to get the style from symbolizer, but I am getting error. Like openlayer2 I used the filters for AND operation but I stucked. Can someone help me on this??


texts: Access is Denied error on ProjectEx call

texts:
I need to be able to convert a lat/long coordinate into a point coordinate and am having trouble getting my code to work.  I created a test app isloating the code that is giving me problems.  I'm running ArcSDE 9.3.1 SP2 on my local machine and also have it running on another box on the network, Jester. If I run my test code with passing in "localhost" to be the host, it works fine.  If I pass in "Jester", I get the error: System.UnauthorizedAccessException: Access is Denied (Exception from HRESULT: 0x80070005 (E_ACCESSDENIED)) and it occurs on the line where I do geom.ProjectEx.  The S4_SPATIAL geodataService is created from the same database on both boxes.

In the code this snippet came from, I use the same server context created from the same geodata service on Jester and it works fine ( I can get to all my tables, feature classes and locators, I can query them, etc.), except when I try to do geom.ProjectEx. I've looked online for this error and haven't been able to find anything helpful.  Apparently its a pretty generic message.

Does anyone have any idea why it would work on my local server but not on the other?  What am I missing?

static void Main(string[] args)
{

    double latitude = 25.73338;
    double longitude = -80.3207233333334;

    IGISServerConnection conn = new GISServerConnection();

    string host = args[0];;

    conn.Connect(host);
    IServerObjectManager som = conn.ServerObjectManager;
    IServerContext gdbContext = som.CreateServerContext("S4_Spatial", "GeodataServer");

    IPoint pt = (IPoint)gdbContext.CreateObject("esriGeometry.Point");
    pt.Y = latitude;
    pt.X = longitude;

    // Create spatial reference for the latitude and longitude point
    SpatialReferenceEnvironment spatialReferenceEnv = (SpatialReferenceEnvironment)gdbContext.CreateObject("esriGeometry.SpatialReferenceEnvironment");
    ISpatialReference latlongRef = spatialReferenceEnv.CreateGeographicCoordinateSystem((int)esriSRGeoCSType.esriSRGeoCS_WGS1984);
    pt.SpatialReference = latlongRef;

    // Get spatial reference of the map
    ISpatialReference outgoingCoordSystem = spatialReferenceEnv.CreateGeographicCoordinateSystem((int)esriSRGeoCSType.esriSRGeoCS_NAD1983HARN);

    NADCONTransformation datumConversion = new NADCONTransformation(spatialReferenceEnv.CreateGeoTransformation((int)esriSRGeoTransformation2Type.esriSRGeoTransformation_NAD_1983_TO_HARN_FL));

    IGeometry2 geom = (IGeometry2)pt;
    geom.SpatialReference = latlongRef;

    // Project the latitude and longitude to map coordinate system
    geom.ProjectEx(outgoingCoordSystem, esriTransformDirection.esriTransformForward, (IGeoTransformation)datumConversion, false, 0.0, 0.0);

    IPoint newPt = (IPoint)geom;
    Console.WriteLine(newPt.X + " , " + newPt.Y);
}


texts: Changing GP parameter value type dynamically for ModelBuilder?

texts:
I'm building custom geoprocessing tool with ArcObjects where output value can be feature class or table with GPFeatureSchema as default schema. But the type of an output dataset depends on some option. In tool's UpdateParameters method I change schema of the output parameter to GPTableSchema or to GPFeatureSchema according to this option. But in ModelBuilder any tool after the my one always thinks that my tool produces feature classes.

Is it possible to dynamically change type of schema with ArcObjects so chained tools get right output information from my tool?


texts: QGIS Crashes using Smooth Geometry Generator on Lines from GPKG

texts:
I'm trying to use the geometry generator in symbology for a line layer. Specifically I want to smooth out the lines using smooth($geometry,5).
Unfortunately QGIS crashes every time I try to implement this. The data is stored in a GPKG. The GPKG contains other layers, but I have exported this particular layer to a separate GPKG file for testing.
It works fine if I export the lines to shapefile, but the data needs to be stored in a GPKG so I need to get this working.
Interestingly, the original GPKG file also contains a layer of polygons and I have successfully implemented the smooth($geometry,5) for these without any issues.
Happy to send through the GPKG file and see if someone can figure out what the problem is. I've run the geometry checker (Check Validity tool) on it and doesn't show any issues.
I'm running version 3.12.3 of QGIS.

Crash Report:
Crash ID: 86223684825a81f552f396a7391a484f166dfa07

Stack Trace
QVector::QVector qvector.h:358
QgsLineString::QgsLineString :
qgis::make_unique qgis.h:425
smoothCurve qgsgeometry.cpp:3091
QgsGeometry::smoothLine qgsgeometry.cpp:3177
QgsGeometry::smooth qgsgeometry.cpp:3046
fcnSmooth qgsexpressionfunction.cpp:2463
QgsStaticExpressionFunction::func qgsexpressionfunction.h:447
QgsExpressionFunction::run qgsexpressionfunction.cpp:102
QgsExpressionNodeFunction::evalNode qgsexpressionnodeimpl.cpp:901
QgsExpressionNode::eval qgsexpressionnode.cpp:29
QgsExpression::evaluate qgsexpression.cpp:373
QgsGeometryGeneratorSymbolLayer::render qgsgeometrygeneratorsymbollayer.cpp:200
QgsSymbol::renderUsingLayer qgssymbol.cpp:716
QgsLineSymbol::renderPolyline qgssymbol.cpp:1967
QgsSymbol::renderFeature qgssymbol.cpp:915
QgsFeatureRenderer::renderFeature qgsrenderer.cpp:130
QgsVectorLayerRenderer::drawRenderer qgsvectorlayerrenderer.cpp:328
QgsVectorLayerRenderer::render qgsvectorlayerrenderer.cpp:284
QgsMapRendererParallelJob::renderLayerStatic qgsmaprendererparalleljob.cpp:353
QtConcurrent::MapKernel::iterator,QtConcurrent::FunctionWrapper1 >::runIteration qtconcurrentmapkernel.h:69
QtConcurrent::MapKernel::iterator,QgsImageOperation::ProcessBlockUsingPixelOperation >::runIterations qtconcurrentmapkernel.h:78
QtConcurrent::IterateKernel::iterator,void>::forThreadFunction qtconcurrentiteratekernel.h:256
QtConcurrent::ThreadEngineBase::run :
QThreadPoolPrivate::reset :
QThread::start :
BaseThreadInitThunk :
RtlUserThreadStart :

QGIS Info
QGIS Version: 3.12.3-Bucure?ti
QGIS code revision: 8234261527
Compiled against Qt: 5.11.2
Running against Qt: 5.11.2
Compiled against GDAL: 3.0.4
Running against GDAL: 3.0.4

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

Happy to supply any other information as needed.


texts: Is it possible to use Python to perform network analysis on only part of a network?

texts:
I have a large network dataset and want to perform an analysis on it which would output the start/end point, the distance between them and whether they're first, second or third (etc.) nearest neighbours. There are approximately 30,000 nodes in the dataset, so I thought it might be possible to process 1,000 at a time (still comparing to all 30,000 nodes, but only doing the analysis on part of the dataset). I have an idea of how to do it (using a for loop, creating new files for each part of the analysis and putting it together) but I'm not certain as to how to create a for loop to use only part of the data and to go through and eventually finish analysis on all data.

Does anyone know if this is possible?

Edit: I've been having thoughts about the partitioning by rowID and was wondering then if I were to be making a service area analysis layer if then the layer would be impossible to create due to the partitioning or otherwise incorrect because of it?


texts: "Does not exist" error when loop through directory

texts:
I need to write a script that loops through shapefiles in a folder and adds fields to each shp. It works for one file at a time, but the added for loop go through the folder is throwing an error that I have included after the code. The print statement prints the name of the first file in the folder, then the next line (fieldsShp02 = arcpy.ListFields(file)) throws the error, which says it does not exist.

Why is it making that complaint?

>>> import arcpy
... import arcgisscripting
... import os
... import sys
... import types
... import locale
... #Input FC
... sourceFC = r"G:\Backup.gdb\muni_provided_utilities_points_copy"
... rootdir = r"G:\Orland Park\CH EAST ORLAND^WESTORLAND2"
... fieldsShp01 = arcpy.ListFields(sourceFC)
... newFields = []

... for file in os.listdir(rootdir):
... #Output shp
...     if file.endswith(".shp"):
...         print "The file does exist:", file
...         fieldsShp02 = arcpy.ListFields(file)
...         for field01 in fieldsShp01:
...             exists = False
...
...             for field02 in fieldsShp02:
...                 if (field01.name == field02.name and field01.type == field02.type and field01.length == field02.length):
...                     exists = True
...                 if (exists == False):
...                     if (field01.name != 'OBJECTID' and field01.name != 'SHAPE' and field01.name != 'CHECKED'):
...                         newFields.append([field01.name, field01.type, field01.length])
...             exists = True
... for field in newFields:
...     fieldnames = field[0]
...     trimmed_fnames = fieldnames[0:10]
...     arcpy.AddField_management(file, trimmed_fnames, field[1], "", "", field[2], "", "NON_NULLABLE","NON_REQUIRED","")

...
The file does exist: AbandinedSanitarySewerForceMain.shp

Runtime error
Traceback (most recent call last):
  File "", line 17, in
  File "j:\apps\arcgis\desktop10.4\arcpy\arcpy__init__.py", line 1138, in ListFields
    return gp.listFields(dataset, wild_card, field_type)
  File "j:\apps\arcgis\desktop10.4\arcpy\arcpy\geoprocessing_base.py", line 346, in listFields
    self._gp.ListFields(*gp_fixargs(args, True)))

IOError: "AbandinedSanitarySewerForceMain.shp" does not exist


texts: Why do merged rasters in QGIS look 'patchy'?

texts:
Poor title for this question, but I can't think of a better way to explain it.

I have merged two rasters. I want the joins to be seamless, but there is an obvious line where they meet...

I've tried playing around with the histogram buttons, but don't seem to be able to make them any better. Any ideas how I should change the rasters before merging?


texts: How to replace pixel values in a single band DEM?

texts: I have a single band DEM as a GeoTiff. I would like to replace values above or below a certain threshold with 0. Can I do that with GDAL?


texts: In DEM data there the elevation ranges from - 156 to 3877. What does this mean?

texts: The elevation of that (project) particular landform cannot be negative (I think). Does this mean the DEM data is wrong?


texts: The Filter geometry function in T-SQL is not precise and does not return valid results

texts:
I have two tables, edges (a set of polyline features):

and nodes (set of points):

I want to find the "id" of those points that have an intersection with the first point of any edges. I used the STIntersects() function in SQL and wrote down the below script :

select EdgeTable.id, NodeTable.idj,
from EdgeTable inner join NodeTable
ON NodeTable.GeomJ.STIntersects(EdgeTable.StartPointG) =1

This code works correctly, but the problem is with their performance. As a matter of fact, EdgeTable and NodeTable have about 2 million records and the execution time of the above script is about 55 minutes which is not at all suitable for my work.

Thus, I found Filter geometry function which seems to be a good solution for improving the performance of this task. In this regard, firstly I created the spatial index on my tables and then use this function as below:

select EdgeTable.id, NodeTable.idj,
from EdgeTable inner join NodeTable
ON NodeTable.GeomJ.Filter(EdgeTable.StartPointG) =1

But this function returns completely wrong results. For example, I tested the script for the edge 1 (idj = 1) which the start point of this edge intersects with node 1 , this is the result of STIntersects() function :

and this for Filter function:

Due to description of Filter geometry function, it remarked that this method is not deterministic and is not precise. Hence, why the results of these two functions are different and how can I come up with it?


texts: calculate distance of points from reference points along a line

texts:
I have a points layer and a multi-part non-overlapping line layer (a line divided into segments). All the points fall exactly on the line. I have the coordinates of the start and endpoint of each segment stored in the segment layer. I want to calculate the distance of each point from the starting point of the segment on which it lies.

In the attached figure, the orange line is the track which is segmented at black lines. The number in red is the FID of each segment. All Green points fall on the line. I want to find the distance of each point from the start/end of a segment (any direction will do).


texts:
I have two INEGI files from Mexico (1990 and 2010) that I want to use the same projection. The 2010 file works fine (see below) and I would like the 1990 file to match the 2010 one. The shapefiles are of all of Mexico and have detailed information about various cities.  Right now, the 1990 prj file reads:

"GEOGCS["GCS_ITRF_1992",DATUM["D_GRS_1980",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]"

I'd like the file to match my 2010 file, which reads:

Projected Coordinate System:    Conica Conforme de Lambert
Projection: Lambert_Conformal_Conic
False_Easting:  2500000.00000000
False_Northing: 0.00000000
Central_Meridian:   -102.00000000
Standard_Parallel_1:    17.50000000
Standard_Parallel_2:    29.50000000
Scale_Factor:   1.00000000
Latitude_Of_Origin: 12.00000000
Linear Unit:    Meter

Geographic Coordinate System:   ITRF92
Datum:  D_GRS_1980
Prime Meridian:     Greenwich
Angular Unit:   Degree


texts: Show string that contains specific number?

texts:
For example, how to show all the strings that contain '841' on this layer?

I tried to use "NAME" LIKE '841' but somehow it did not just show 841_1Q and 841_1, instead it still showed every point on this layer just highlighting those two containing 841.

I tried using a wildcard but the "Output preview" shows 0 result.


texts: How do I represent the whole Earth as a Polygon?

texts:
This is related to determining the interior of a Polygon.

My difficulty is that if I specify a Polygon covering the Earth going from East to West, it will still have a boundary along the antemeridian, which would show up in azimuthal projections. For example, in GeoJSON:

{
  "type": "Polygon",
  "coordinates": [
    [[-180, -90], [-180, 90], [180, 90], [180, -90], [-180, -90]]
  ]
}

Is there a standard way to specify a whole-sphere Polygon, with no boundary?

I've thought about using an empty exterior ring:

{
  "type": "Polygon",
  "coordinates": [
    []
  ]
}

The logic being that I could specify further rings as holes to be punched into the whole-sphere Polygon.

Unfortunately, GeoJSON stipulates that LinearRings must have at least 4 points, so this would be straying from the specification.  However, I'm still interested in understanding whether this approach is taken anywhere else.


texts: ArcGIS Online attribute hyperlink is not reading the space in URL, how to fix?

texts:
We are using Arc GIS Online to display some point data.  In ArcGIS Online we configured the attributes to point to a hyperlink that is located in our gdb. Our GDB contains a field that points to where we are hosting an image of this point (pic2)  For some reason this link is ending at the first space even though the image name is in the format NAME STREET NUMBER.tiff with spaces.  This format has worked perfectly in the past and arc was able to pull through the spaces without a problem and link to the tiff files, but now it has stopped working.  We tried to configure the custom attribute display to link to {HYPERLINK} {HYPERLINK%20} and other combinations but it is still not pulling the proper link through.

Why has this stopped working now, is there a solution short of either changing the file names to remove any spaces or changing the links in HYPERLINK field in the GDB to something like WESTMONT%20ST%20003.0.tiff .  Is there something setting we can change or a way to change it in the custom attribute field?

(source: gislab.net)
!

edited for grammar


texts: Is there free GIS reader that supports transparency of JPEG layers?

texts:
I need to make from my ArcMap layers (which are georeferenced JPEG files) a file that can work on a free and easy to use program.

The program needs to be able to open multiple layers so we can put them on top of each other and to be able to change the transparancy of the layers.

Does anybody know a solution?


texts: Downloading the info of a QGIS WMS Url for reading it as a spreadsheet file

texts:
I've recently started using QGIS and learned how to connect my map to a server through a WMS connection. Is there a way to download this data as a more "spreadsheet-like" way?
The URL is:
http://acervofundiario.incra.gov.br/i3geo/ogc.php?tema=imoveiscertificados_privado_go


texts: What is the correct way to show 10000 maps/tiles on my contact address list web page?

texts:
I have a webpage with a grid of contact list (around 10K people).
I want to put in one of the columns a small map/tile of the address from a different column.
Do I need to generate 10k maps object and zoom to geocode address or is there an alternative solution to retrieve an image of the address coords and show it as a tile ?


texts: What's a good toolchain to reproject some PNGs in plate carre to UTM?

texts:
I have some PNGs produced from Canadian DEM data. They're 1201px by 1201px heightmaps. Each one represents one quarter of a degree latitude and one quarter of a degree longitude.

But I want to work in UTM, and it's looking like it's a huge pain to convert these. Ideally I'd love to use a tool that's a bit simpler than something like GRASS GIS, but if a big suite like that is the tool for this job, maybe I could ask for a walkthrough on what settings to be using?

Thanks.


texts: Can FileGDB SDK be made to support custom coordinate systems?

texts:
In May 2012, a question was asked about File Geodatabase (*.gdb) support in QGIS?

In the comments there, it was pointed out by @MattWilkie that the FileGDB SDK does not support custom coordinate systems.

Has there been a workaround developed for custom coordinate systems?


texts:
I need your help. I'm trying to calculate the global irridiation for a specific area. I'm interested in the results for rooftops in special.
I have a DEM, calculated slope, aspect and horizons and started r.sun for every day with a monthly changing Linke Factor. All calculated maps (slope, aspect, horizons) seem to have believable values.
My Problem is: On the yearly sum (addition of all 365 days) I get values for total global irridiation of about 1900 KW/m². Which is 1,9times more then expected. Especially when the solar constant is only 1368 kW/m²....
Could you imagine what went wrong? I'm lost.


texts: Is there a way to calculate zonal entropy index for a raster?

texts:
I have an elevation raster and a grid shapefile:
What I want to do is calculate zonal statistics of the raster based on the grid features. The "statistic" I want is the entropy index (Shannon entropy) and it is not included in the tools of most popular software (ArcGIS and QGIS). I searched for a python solution and I found the user-defined statistics of rasterstats library: https://pythonhosted.org/rasterstats/manual.html#user-defined-statistics.
I attempted the following:
from rasterstats import zonal_stats
from scipy.stats import entropy

def my_entropy(x):
    return(entropy(x,base=2))

result=zonal_stats("...path/to/grid.shp","...path/to/elevation.tif", add_stats={'my_entropy':my_entropy})

The code runs, but the result is odd in terms of dimensions. While the length of the result is the same as the number of features of the grid shapefile, the calculated entropy of each feature is an array of values, and not a single value:
Any help is appreciated. A solution based on software (with a model or arcpy) is welcome too, although after making some experiments, I couldn't figure out something efficient (it requires the production of histogram tables for each grid feature etc.)


texts:
I need to get the intersection of two polygons and both sides use the same dataset (titlenumberpolygon - UK National polygons). PostGIS is used in the analysis.

First a spatial index was created which took 27 mins. There are 30 million polygons.

CREATE INDEX geoid ON titlenumberpolygon USING GIST(geom)

The following code is used and for both the tables the polygon data is coming from titlenumberpolygon table.

SELECT tb2.UPRN,tb1.titleno
INTO AB_VOA_Onwer1
FROM (
    SELECT t1.titleno, tp.*
    FROM Ownership as t1
    INNER JOIN titlenumberpolygon as tp
    ON t1.titleno=tp.title_no) as tb1
INNER JOIN(
    SELECT AU.UPRN,AU.Postcode,tp.*
    FROM AB_VOA_UPRN as AU
    INNER JOIN titlenumberpolygon as tp
    ON AU.inspiredid=tp.poly_id) as tb2
 ON ST_Intersects(ST_MakeValid(tb1.geom),ST_MakeValid(tb2.geom));

It takes over 2 hours to run for one postcode and there are 500,000 postcodes in the table.

How can I optimise this code? Maybe a different spatial indexing method as I am using the same table in the intersection.


texts:
I have a task to get the parcel numbers along a specific way.
The way is a vector layer with specific points, and the parcel numbers are supported via WMS server.
The task is splitted in to two tasks:
First is to write down all the parcels that are crossed by the way,
second for given parcel write out all of points.

EXAMPLE:
First task: way goes
135/3; 135/5; 503; 333/7; 333/6; 332/2; 499; 613; 498; 303; 302/2
Second task: parcels containing points
135/3 - 19/L-975; 333/7 - 20/L-975; 333/6 - 21/L-975; 332/2 - 22/L-975; 613 -  23/L-975 & 24/L-975

The manual writting down all of the parcels takes much time. And for a large regions is problematic.
Is it possible to create some kind of rule or script to extract this data automatically?


texts: How can units be defined in QGIS for the Rasterize tool?

texts:
I am using QGIS 2.4 (Chugiak) via the kyngchaos build for OS-X. My question is--when using the Rasterize tool to convert a point vector layer to a raster and using the "Raster resolution in map units per pixel" option, how are "map units" defined?

So far as I can tell thus far, "map units" are decimal degrees regardless of other settings. The "Canvas units" setting under Project Properties --> General does not change this. Changing the coordinate system under Project Properties --> CRS does not change this. "Map units" are not dependent on the CRS of the input layer, nor on the CRS of the first layer added to the project. Is decimal degrees the only option and not modifiable in any way, or is there some other setting I have yet to find that defines "map units" in this context?

The reason I am curious is that if you have a point shapefile you wish to rasterize, and that shapefile uses a regular grid that defined in a projected coordinate system (e.g., a grid with points at 1000 m in NAD83 Z13) rather than a geographic coordinate system, there is no "Raster resolution in map units per pixel" that will be appropriate except over rather small areas if decimal degrees are the only map units that can be used. What you want to tell QGIS in this context is "rasterize this vector layer with a pixel size of 1000 m by 1000 m in NAD83 Z13" not "rasterize this vector layer with some pixel size of X by Y degrees". It seems there must be a way to get QGIS to do this, but I have yet to find it.


texts: Arcpy InsertCursor.insertRow crashes ArcGIS after many (thousands) insertRow calls

texts:
I am looping over a quite large feature class (~5million records)
While looping using a SearchCursor, I insert new geometries in a FileGeodatabase FeatureClass using arcpy.da.InsertCursor.
The inserted geometries are simple poygons (squares).
The code works flaylessly for about 2k loops when ArcGIS suddenly crashes.
Geometries fed to cursor.InsertRow() before crash are available in target FeatureClass.
I felt like this was some kid of buffer issue so I sought for flushing commands but I could not find any in InsertCursor.

What can it be?

********** UPDATE *********

Adding script for clarification.
My source FeatureClass is a Point FC with 2 text attributes: 'cellx' and 'celly'
My destination FeatureClass isa Polygon FC with geometry plus one text attribute named 'gid'

hc = 0.00833333 # halfcell
# Open an InsertCursor toinsert new geometries
insertCursor = arcpy.da.InsertCursor('meteogrid_polygon', ['SHAPE@','gid'])
i = 1
with arcpy.da.SearchCursor("meteogrid_points_w_gid",['SHAPE@','celly','cellx']) as readCursor:
for row in readCursor:
    point = row[0].firstPoint
    x = point.X
    y = point.Y
    # Create a square polygon geometry (point(X, Y) around source point
    array = arcpy.Array([arcpy.Point(x+hc, y+hc),arcpy.Point(x+hc, y-hc),arcpy.Point(x-hc, y-hc),arcpy.Point(x-hc, y+hc)])
    polygon = arcpy.Polygon(array)
    # Create string to put in 'gid' attribute
    cellname = str(row[1])+"_"+str(row[2])
    # Log to screen for debugging
    print str(i)+" "+cellname
    i = i + 1
    insertCursor.insertRow([polygon,cellname])


texts:
I am trying to do a raster calculation in the QGIS 3.7 python console, using the following code:

import gdal

bf_thickness = "D:\Minnesota data\Data files\clipped\Buffalo_thickness.tif"

bf_capacity = "D:\Minnesota data\Data files\clipped\Buffalo_capacityPyqgis.tif"

bf_thickness_template = gdal.Open(bf_thickness)

driver_tiff = gdal.GetDriverByName ("GTiff")

bf_capacity_new = driver_tiff.CreateCopy(bf_capacity, bf_thickness_template, strict=0)

band_base_1 = bf_thickness_template.GetRasterBand(1).ReadAsArray()

band_new = (band_base_1 * 100 * 30 * 4 * 3.14)/(exponential (0.25 ** 2 * 0.001/(4 * 100 * 30)))

bf_capacity_new.GetRasterBand(1).WriteArray(band_new)

bf_capacity_new.GetRasterBand(1).ComputeStatistics(0)

Here is the error message that I receive:

Traceback (most recent call last):
  File "C:\OSGEO4~1\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 20, in <module>
NameError: name 'exponential' is not defined


texts:
I work in a project in QGIS where I have some layers with points. All these points are stations. For each station I have a table in database with time-series measurements.
I want to find a way to click on each point and have some type of menu from where I can plot a graph based on stored measurements I have in database.
Is there any plugin that can help me with this functionality?


texts:
Is there in arcobjects the ability to store some data at the layer in table of content sth. like

ILayer.Tag or an Interface which I can misuse for this?


texts: Connecting to WFS - no feature types available in parameters

texts: I've installed the interoperability connection folder on arcgis 10.2 and I'm going through the steps to view a WFS. I've got to the WFS Parameters and I'm trying to select the Feature Types from the Constraints box but there are no feature types listed? Would this be because the URL is wrong? Any help please?


texts: Confusion over arcpy.Describe() properties?

texts:
Can someone clarify the meaning of the 'fullPropsRetrieved'and 'metadataRetrieved' properties of a Describe object within ArcPy.

In particular my confusion is regarding the definitions provided on the ESRI website of each of the properties:

fullPropsRetrieved: "Indicates whether full properties have been retrieved"
metadataRetrieved: "Indicates whether the metadata has been retrieved"

What is the difference between these two properties?


texts: Converting feature class with domains to SHP file

texts:
I have a feature class which has coded value domain associated with it. The code and the description are different. For example, I have coded value domain for pipe diameter. For pipe diameter of 8" I have assigned coded value 1. Now, when I convert the feature class to shapefile I understand that it assigns the coded value not the description.

I know that I can convert those values using field calculator but I have 20 coded values. It will take considerable amount of time to do this operation.

Is there a script/tool that I can use which can assign description not the coded value for that attribute column when converting from FC to SHP?


texts: I have a number of satellite image rasters covering my study area. I want to mosaic and clip the output raster to my study area. Note that the images overlap each other and do not join end to end. How can I perform this task?


texts:
Lately I've been using the OGR projection classes that come with ogr/gdal, but pyproj was recommended to me, so I thought I'd give it a try. To help me decide if I should make the switch, I did a speed test. The following is a small (almost) reproducible example that I came up with to test the two. I'm not sure if this test is totally fair, so comments and recommendations are welcome!

Imports first, to make sure we start off with a level playing field:

from pandas import Series # This is what our geometries are stored in
from shapely import wkb
import functools
from shapely.geometry import shape
from osgeo import ogr
# The following two imports are the important ones
from pyproj import Proj, transform
from osgeo.osr import SpatialReference, CoordinateTransformation

Because I'm storing the shapely geometries in a pandas 'Series', the
functions need to work with Series.apply(). Here I define two functions
(one using 'ogr.osr' and one using 'pyproj') to perform coordinate
transformations inside a call to Series.apply():

def transform_osr(geoms, origin, target):
    target_crs = SpatialReference()
    target_crs.ImportFromEPSG(origin)
    origin_crs = SpatialReference()
    origin_crs.ImportFromEPSG(origin)
    transformer = CoordinateTransformation(origin_crs, target_crs)
    def trans(geom):
        g = ogr.CreateGeometryFromWkb(geom.wkb)
        if g.Transform(transformer):
            raise Exception("Ahhhh!")
        g = wkb.loads(g.ExportToWkb())
        return g
    return geoms.apply(trans)

def transform_pyproj(geoms, origin, target):
    target = Proj(init="epsg:%s" % target)
    origin = Proj(init="epsg:%s" % origin)
    transformer = functools.partial(transform, origin, target)
    def trans(geom):
        interface = geom.__geo_interface__
        geom_type = interface['type']
        coords = interface['coordinates']
        result = apply_to_coord_pairs(transformer, coords)
        return shape({'coordinates':result, 'type':geom_type})
    def apply_to_coord_pairs(fun, geom):
        return [not all([hasattr(y, "__iter__") for y in x]) and \
                fun(*x) or apply_to_coord_pairs(fun, x) for x in geom]
    return geoms.apply(trans)

Each of these functions takes an EPSG code as input for the origin and
destination coordinate reference systems. Both libraries offer alternative
ways to define projection information, but EPSG codes keeps the code pretty
simple for now.

Here are the results, using the %timeit magic function in ipython:

In [1]: %timeit transform_pyproj(l, 29903, 4326)
1 loops, best of 3: 641 ms per loop

In [2]: %timeit transform_osr(l, 29903, 4326)
10 loops, best of 3: 27.4 ms per loop

Clearly the 'ogr.osr' version is faster, but is it a fair comparison? The
'pyproj' version is done on individual points, and is mostly run in Python,
whereas the 'ogr.osr' version operates directly on the OGR geometry object.
Is there a better way to compare these?


texts: How to transform location data from Lambert (EPSG:31370) to GPS (EPSG:4326)?

texts:
Sorry if my way of writing is wrong, I'm a data scientist and I have little to none experience in GIS.

I would like to use an open dataset (from http://opendatastore.brussels/fr/ to be precise). It contains mobility data and some locations are encoded in EPSG:31370 (Belgian Lambert 72). I'd like to transform these locations in GPS coordinates (EPSG:4326 if I'm not mistaken).

I'm not sure to understand how you could derive a formula to automatically transform data from EPSG:31370 to EPSG:4326.

Does it exist Python libraries (or R or any other languages) in order to automate these location transformations? I found an interesting website (https://epsg.io/transform) to do it manually but I have +/- 250.000 transformations to do.


texts: ArcGIS JavaScript choropleth map with custom data

texts:
I am using ArcGIS JS API for creating a different type of maps. I am pretty much clear on how do I create point maps using Lat/Long fields. But I don't seem to get a way to create choropleth map like this.

https://developers.arcgis.com/javascript/latest/sample-code/visualization-vv-color/index.html

I want to create maps similar to this. It will be limited to the USA but it can be by state or county or zip. From where I am supposed to get geometry information for these regions so that I can use them with my own data and then host a feature layer on ArcGIS online to use in my application? All the examples provided by ArcGIS uses some hosted layer like https://services.arcgis.com/V6ZHFr6zdgNZuVG0/arcgis/rest/services/counties_politics_poverty/FeatureServer/0
but I don't know what do I have to do to host something similar with my own data to create a feature service.

By now I am very familiar with the ArcGIS API, the only thing I don't understand is from where I am supposed to get geometry information for say USA states.


texts: Has anyone used QGIS to plot 3D data (temperatures profiles at different depths at sea)?


texts: resample etopo.grd with gdal_warp gives many null values

texts:
I'm having a problem with resampling ETOPO NetCDF GRD to a bigger pixelsize for using it with Generic Mapping Tools.

I'm using gdal_warp for the resampling (just for testing with pixelsize of 1°):

gdalwarp -overwrite -s_srs EPSG:4326 -tr 1 1 -r near -of netCDF input_etopo.grd etopo_resampled.nc

but the output grd gets many null-values. If I use the geotiff of etopo the output is the same.

Does anybody had the same problem and knows how to prevent it?
Different resampling methods doesn't solved the problem.

Etopo-Input:

Etopo-Resampled:


texts: Draw rectangles with Shape Digitizing Toolbar and using a snapping function in QGIS 3

texts:
I draw rectangles with the Shape Digitizing Toolbar. In QGIS 3 there is a Snapping Toolbar which provides a snapping function. But unfortunately it is not possible to use both together. So whenever I want to draw a rectangle I can not snap it onto a line which leads to incorrect results.

Here is the same problem but it occured in QGIS 2.

Any solutions for QGIS 3?


texts: ArcGIS Raster Calculator fails with Error 000539

texts:
I have been using Raster Calculator in ArcMap 10.2, today it started returning error 000539 for all actions. Previously it had worked fine. I have not made any change ot my ArcGIS, Python, or Windows 7 installation.

This is the error message I get:

Messages
Executing: RasterCalculator ""Carbon"  + 2" C:\Users\mem48\Documents\ArcGIS\Default.gdb\rastercalc6
Start Time: Tue Mar 17 15:48:30 2015
Raster(r"Carbon")  + 2
ERROR 000539: Error running expression: rcexec()
Traceback (most recent call last):
  File "<expression>", line 1, in <module>
  File "<string>", line 2, in rcexec
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\__init__.py", line 24, in <module>
    from arcpy.toolbox import *
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\toolbox.py", line 351, in <module>
    from management import Graph, GraphTemplate
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 22, in <module>
    import _management
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\_management.py", line 14, in <module>
    import _graph
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\_graph.py", line 27, in <module>
    import numpy
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\numpy\__init__.py", line 137, in <module>
    import add_newdocs
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\numpy\add_newdocs.py", line 9, in <module>
    from numpy.lib import add_newdoc
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\numpy\lib\__init__.py", line 4, in <module>
    from type_check import *
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\numpy\lib\type_check.py", line 8, in <module>
    import numpy.core.numeric as _nx
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\numpy\core\__init__.py", line 5, in <module>
    import multiarray
ImportError: DLL load failed: %1 is not a valid Win32 application.

Failed to execute (RasterCalculator).
Failed at Tue Mar 17 15:48:33 2015 (Elapsed Time: 3.37 seconds)

All the files mentioned in the message are present. I've tried using different inputs, output locations, and expressions. It seems that Raster Calculator has just failed completely.

What is wrong and how do I fix it?


texts: Catching invalid expression errors in Python script tool?

texts:
I have a script that creates a new feature using arcpy.GetParameterAsText, but I'm having trouble creating an error message for incorrect user input.  The input for a parcel # should be something like '1500300000029', but what if 'sdfjfh' is entered?  As it stands no error is thrown and a empty feature is created.  I've also tried inside PythonWin by assigning a bad entry for 'input' and it runs without error.

import arcpy
import sys

arcpy.env.overwriteOutput = True
try:
    Parcels = "..path to data"
arcpy.env.workspace = "C:/Data_new/Temp/Default.gdb"

#input = arcpy.GetParameterAsText(0)
input = "dfdsgdg"

expression = "PARCEL_ID = " + "'" + str(input) + "'"
arcpy.MakeFeatureLayer_management(Parcels, "parcels_lyr")
arcpy.SelectLayerByAttribute_management("parcels_lyr", "NEW_SELECTION", expression)

arcpy.CopyFeatures_management("parcels_lyr", "output")
except:
    print arcpy.GetMessages()

I've also tried

except StandardError, ErrDesc:
    arcpy.AddWarning(input)
    arcpy.AddWarning(ErrDesc)

except:
    arcpy.AddWarning(input)
    arcpy.AddWarning(ErrDesc)

What is the appropriate way to handle this?


texts: Importing QGIS script parameters from another script

texts:
I have 2 scripts within QGIS processing toolbox and I can't seem to work out how to access the parameters of 1 script from the other. Since QGIS sets up the processing parameters (by 'looking for' the ## notation) before the rest of the script is run, the parameters are already set up by the time the second script is called. Has anyone found a work around for this?

script a:

#script a
import imp, os.path

def import_(filename):
    (path, name) = os.path.split(filename)
    (name, ext) = os.path.splitext(name)

    (file, filename, data) = imp.find_module(name, [path])
    return imp.load_module(name, file, filename, data)

foo = import_("path to b")
foo.getParameters()

and script b...

#script b

def getParameters():

    ##parameter_1=string
    ##parameter_2=number
    return

The reason I'm trying to do this is because QGIS makes a local copy of the script to run when a user clicks 'add script from file'. So, if I later want to change the input parameters, I would need to inform the user that the script is updated and that they need to 'add script from file' again to get the most recent version. If I can create a script for local copies that just calls the main server script, I can make updates without having to get all the script users to manually change to the latest version.

I know plugins are a more appropriate way to control this kind of workflow, but in the corporate I'm working on I don't have (and can't get) Qt installed to make the GUI.


texts: Difference between Viewshed and Observer Points?

texts: Can someone explain to me in layman's terms the difference between Viewshed and Observer Points?


texts: Deleting attribute fields through SQL Server Management Studio

texts:
I have a problem. After appending annotations into the feature class on the SQL Server database I need to delete extra attribute fields that where made in the process. Usually I would just simply use Delete Field tool. However because people are working on this database I can't do this because of unique lock. I also have Python script that disconnects everyone and delete these fields during the night, but I want to find the way to do this when people are connected. So I had an idea that I could delete these fields through SQL Server Management Studio but I'm afraid that this may cause problems to database or text feature class.
Is it safe to delete attribute fields through SQL Server Management Studio?
I'm using ArcGIS Desktop 10.6.1 and SQL Server Management Studio 18.


texts: Using python script in ArcGIS field calculator?

texts:
I am using ArcGIS 10.3 version in desktop. I wanted to assign landuse names (Landuse_04 column) based on the codes those are in my attribute table as the column ID_MAP. once I am done with editing the script in python, field calculator, this dialogue box is appearing! I do not have any idea about python or programming just worked on a prepared script.

Does anyone know why is this occurring?


texts: Plotting Geo Points with "by arc centered at" Instructions

texts:
A recent Supreme Court decision redefined the legal boundaries of California along the Pacific Ocean: http://www.theverge.com/tldr/2014/12/15/7397085/supreme-court-california-boundaries-map
I am confused how you would plot the points as they are described in the Supreme Court decision.

BEGINNING AT 482577.890 3599275.555
BY ARC CENTERED AT 488133.576 3599216.475
TO 482623.800 3599931.673

I am assuming these are bezier curves in a CAD system but I am not certain how to do this calculation in QGIS or OGR using a spatial reference system (EPSG:26910).
Any thoughts?


texts: Attribute Assistant-Dynamic Value Table-TableName Values- Expressions ValueInfo

texts:
I am having a problem populating an ID attribute "GISOBJID" which is a string in a versioned SDE data set. The AUTONUMBER I want to combine "ASSETID" works fine for everything , but I want to combine that with an ALPHA prefix (FM-123) into the "GISOBJID" attribute. I have tried several ways and no success. Only the AUTONUMBER works. I have tried using the full SDE name "UTILITIES2.CCU.ssPressurizedMain" in the DynamicValue Table and the shotened version "ssPressurizedMain" with no results. I have also tried an EXPRESSION "FM-" & Str( [ASSETNUM] ) and GenerateID with SSPRESSURIZEDMAIN|0|FM-[seq] still no success.

[ 2

[ 3


texts:
I work with QGIS 3.14 and should do an infographic map with pie charts in black and white (grayscales) only. The final result should look something like this.
I tried to use the built in pie chart rendering from layer properties. However, in these pie charts, you can only define a color, but not a pattern fill or similar styling. As I have about 12 different categories, styling them in different shapes of gray makes no sense as the eye can't distinguish the difference between the categories. Each pie, however, has just about 3 or 4 sectors represented (the rest will be styled as "others").
I looked at this workaround, but hope that in the meantime, there is a better solution.
Any idea how I could get a grayscale pie chart?


texts: The MassGIS parcel standard (http://www.mass.gov/mgis/ParstndrdVer1_5_1.pdf) utilizes a concatenation of the whole number portion of x- and y-coordinates to create a unique id (LOC_ID) for features.  I'm considering doing the same for a point feature class.  I like the consistent methodology, but perhaps I'm overlooking something.  Is there a standard or best practice for creating unique ids for point features?


texts: Moving multiple features to snap to nearest edge

texts:
I have thousands of long rectangular polygons which I need to align to a specific/nearest edge.  The screenshot below demonstrates what I need to achieve:

Yellow and white dotted features need to maintain original width (1m) but move and snap to the nearest edge.

I am working with ArcGIS 10.3.

I have looked at all of the Alignment Tools, but I'm just not clear on the best method.  I tried Align to Shape but this adjusts the width of the polygons.

Any clues?


texts: How to disable a base layer in OpenLayers based on zoom level

texts:
Possible Duplicate:
  How can I switch baselayers as a user zooms in to the map in OpenLayers

I have a map that has two base layers and several overlays.  One of the base layers only has 5 zoom levels.  I would like to disable that base layer in the layer switcher when the map is zoomed outside of the range of that base layer.  I have found how to do that with overlays in this forum at this link.  I changed dataLayers to baseLayers.  My problem seems to come when the base layer that I am trying to limit is the base layer selected when I zoom outside the range of that base layer.  If I then select the other base layer, it appears as if the zoom control is locked up.

Any suggestions are greatly appreciated.

Thanks.


texts: Using ModelBuilder to join table to shapefile in ArcMap?

texts: I have more than 40 shapefiles and tables. I want to join the table of attribute of these shapefiles with related  tables permanently. For small number of shapefiles and tables, I normally join shapefiles table of attribute with related table and then export the shapefile and save it as a new name.But how can i do it in ModelBuilder for large number of shapefiles?


texts: Finding largest number of observations within a given radius

texts: I have the following problem: Given the location of insured properties (say houses with fire insurance), i have to find the largest sum of insured values which fit into a 300m radius. Now, my first approach would be to divide the map into grids of appropriate size(e.g. 55m), assign the properties to the grids and calculate the sum of the grids in the first step, and then aggregate blockwise (10x10 grids), and finally search for the largest block. But I have a fealing that there must be a more elegant approach.


texts: gdal translate resampling -r flag has no effect?

texts:
Whenever I add a resampling parameter (-r) to my gdal_translate commands, it seems to have no effect.  For example the following two commands produce a binary exact result.
How can I get my image to have bilinear resampling?  I am trying to smooth out my SRTM 4.1 (DEM) data.
GDAL 3.0.4
gdal_translate -of GTiff -r nearest D:/srtm_41_03.tif D:/bulgaria-nearest.tif

gdal_translate -of GTiff -r bilinear D:/srtm_41_03.tif D:/bulgaria-bilinear.tif

Source tile available here (srtm_41_03).


texts: Error importing UTM delimited text layer into QGIS

texts:
I'm trying to import a small delimited (tab separated) text into QGIS. The coordinates are in UTM format.

Sample data:
id  N   E
PT-001  7237229 146872
PT-002  7237230 146897

The first column contains labels, the other columns (E and N) contain coordinates (numeric integers). I've specified the UTM zone in the Geometric CRS (EPSG:31983), so it seems I don't need it inside the file (there's not even an input field for it).
I'm getting an error message:

Layer is not valid: The layer file:///home/....csv?type=csv&delimiter=%5Ct&maxFields=10000&detectTypes=yes&xField=E&yField=N&crs=EPSG:31983&spatialIndex=no&subsetIndex=no&watchFile=no is not a valid layer and can not be added to the map. Reason:

The fact that it gives no reason is just a pain. What is possibly happening?


texts: OpenStreetMap export changes depending on browser?

texts:
I created an export in OpenStreetMap for an address in Melbourne, Australia with a zoom level of 14: http://www.openstreetmap.org/export#map=14/-37.8593/145.0390. The export coordinates changed depending on whether I view it in Firefox, Chrome or Microsoft edge. Why does this variation happen?

I assumed that the export would be the same given the same coordinates and the same map size/zoom level (i.e. 14). Is there a way to achieve that?

Chrome:

FireFox

Edge


texts: How to round up/down an pixel value of a raster to a pixel value with 2 decimal places?

texts:
I have a raster (DTM) with floating points as pixel value. Now i want to round up/down the pixel values (to its next nearest value) to 2 decimal places.

I use ArcGIS (with extensions) and tried the tools "int" and "round up/down" to do this but i haven't found any option to determine the number of decimal places. My next idea was to modify one of these tools within the raster calculator but unfortunatelly i don't have (enough) python skills - also with the help of google.

Do you have any easy solutions to solve this problem?


texts: Symbolizing values when GeoServer published GeoTiff is all black?

texts:
I am new to GeoServer. I'm using the latest version as of time of post on windows 10 running as a windows service.

I have published a single geotif. Each pixes has a value between 0 and 7, and I want to see each value as a color in the Raster. When I add the geotif to ArcMap, somehow ArcMap knows what the colors are.

How do I tell geoserver to do symbology on a Raster (Geotif file on my local disk)?


texts: Definition Query to Python script on Oracle 9i and ArcGIS 9.3.1

texts:
SUBTYPECD=1 AND SERIESNAME NOT LIKE '%-%' AND OPERATIONALSTATUS IS NULL OR SUBTYPECD=1 AND SERIESNAME NOT LIKE '%-%' AND OPERATIONALSTATUS = 976

The above is my def query in ArcMap against ArcSDE 9.3.1 on oracle 9i with python 2.5.  I want to be using the gp.makefeaturelayer with the sql parameter filled in with a rendition of the above.

ps I started to ask this question because it seemed daunting and I am in a crunch..then figured it out so I thought I'd finish the post for those to come.


texts: Geometry Transformation in Spatialite 4.1.1


texts: Aggregating multiple rasters within a polygon in PostGIS

texts:
I have several raster tables in PostGIS (four in fact but in this case let us assume there are just two) - same cell sizes, same alignment but with different types of value (e.g 3d_pot, 3d_area). Occasionally there will be -9999 values where data is missing but this might not be the case in all three datasets for the same cell (e.g. there could be 3d_pot but not 3d_area).

I want to calculate the aggregate statistics for all the cells that are within a polygon (e.g. a building footprint) but only for those cells where 3d_pot is greater than a value (say 100).

For one value I need to calculate a combination of the cell values at the same location (3d_area * 3d_pot).

I'm struggling to find a fast way to effectively join the rasters based upon their spatial location, having intersected with the polygon layer and ignore all cells if one of them is -9999.

What is the simplest SQL way of doing this on the fly (without generating new rasters or rasterizing the building polygons?)

Here is what I currently have but this feels wrong and excessive to do multiple Intersections. Can you join one raster to another just using the raster geometry?

WITH solarpot AS (
SELECT fid,
    (ST_Intersection(T.wkb_geometry, ST_Clip(ST_MapAlgebra(AREA.rast,POT.rast,'([rast2])',NULL,'INTERSECTION'), 1, T.wkb_geometry, True))).val AS potval,
    (ST_Intersection(T.wkb_geometry, ST_Clip(ST_MapAlgebra(AREA.rast,POT.rast,'([rast1])',NULL,'INTERSECTION'), 1, T.wkb_geometry, True))).val AS areaval,
    (ST_Intersection(T.wkb_geometry, ST_Clip(ST_MapAlgebra(AREA.rast,POT.rast,'([rast1]*[rast2])',NULL,'INTERSECTION'), 1, T.wkb_geometry, True))).val AS sumpotval
   FROM solar_rasters.solarp_3d_area AREA
   JOIN os_gla.topographicarea T ON ST_Intersects(T.wkb_geometry, AREA.rast)
   JOIN solar_rasters.solarp_3d_pot POT ON ST_Intersects(T.wkb_geometry, POT.rast)
   WHERE T.fid = 'osgb5000005156542071'
)
SELECT fid, COUNT(potval), SUM(potval), SUM(areaval) AS sumareaval, SUM(sumpotval) AS sumpotval, AVG(potval), stddev(potval), MIN(potval), MAX(potval)
FROM solarpot
WHERE potval IS NOT NULL AND potval >= 100 AND areaval IS NOT NULL
GROUP BY fid


texts: Contour Plugin generated contours do not have enough nodes

texts:
I am using the contour plugin to generate intermediate contours.  I have a vector file with 5 foot contours and want to interpolate 1 foot contours from it.  The procedure I used was:

Extract nodes from my 5 foot contours vector file.

Create intermediate contours from nodes extracted from original 5 foot contours.

Nodes from intermediate contours.

The problem is that the intermediate contours are not very good in some of the open areas. I assume that this is due to the nodes of the intermediate contours the open areas are too far apart.

Is there a way to get the contour plugin to generate more nodes in these open areas to make the generated intermediate contours in these open areas better?


texts:
while clicking on the button on.click event is not fired and getting error "Cannot read property 'on' of undefined",how to solve this issue.

    Simple Map

      html, body, #map {
        height: 100%;
        width: 100%;
        margin: 0;
        padding: 0;
      }
      body {
        background-color: #FFF;
        overflow: hidden;
        font-family: "Trebuchet MS";
      }

      var map,navToolbar;

  require(["esri/map",
           "dojo/dom",
           "dojo/on",
           "dijit/form/Button",
           "dojo/domReady!"], function(Map,dom,on,parser,registry) {
      map = new Map("map", {
      basemap: "streets",
      center: [72, 25.75],
      zoom: 5
    });
    //navToolbar = new Navigation(map);
      //on(navToolbar, "onExtentHistoryChange", extentHistoryChangeHandler);

    dijit.registry.byId("clickMe").on("click",showMessage);
    function showMessage(){
     alert("Hello")
    }
});

</script>

     Click Me


texts: Creating polygon with esrieditor.editor in C#?

texts:
How can I make a polygon using esrieditor.editor in ArcObjects?

I tried this method:

workspaceEdit.StartEditOperation();
IFeature feature = featureClass.CreateFeature();
Polygon p2 = featureClass.GetFeature(0).Shape as Polygon;

feature.Shape = p2 as IPolygon;

ISubtypes subtypes = (ISubtypes)featureClass;
IRowSubtypes rowSubtypes = (IRowSubtypes)feature;

rowSubtypes.InitDefaultValues();
feature.Store();


texts: How can I preserve ObjectIDs in a Featureset returned from an SOE?

texts:
I'm returning a featureset from an SOE.

To create the FeatureSet, I'm first creating a Recordset, inserting features into it, and then serializing it as json using SOESupport.  The client then deserializes it as a featureset.

The resulting featureset has different OIDs.  Is there a way to preserve the OIDs of the features I'm inserting into the recordset?

When I perform a query using the REST API, the featureset returned has the original objectID's.

IMapServer3.QueryData might do this, but I don't see a way to pass it an array of OIDs as I can with IFeatureClass.GetFeatures.


texts: QGIS 2.14/2.18 refuses to work after a fresh install

texts:
I am having trouble getting QGIS to work on my computer. I was working off 2.14.14 and I thought it would be a good idea to get hold of 2.14.17. Thereafter, I have not been able to fire up QGIS at all. When I tried reinstalling thereafter, I kept getting a message saying that there was QGIS already installed on my machine whereas I had already uninstalled it prior to the installation process. The desktop icons do not get created and there is no qgis.bat file in the bin folder. If I click on qgis-bin.exe in the bin folder, it returns an error saying qgis_app.dll is missing. I am completely stumped. I have uninstalled, rebooted, reinstalled and repeated this cycle at least 6 times in the last 3-4 hours but to no avail. I thought installing version 2.18 would solve the problem. No, it does not. The only thing that was different was that it did not say that 2.18 was already installed on my machine but everything else is the same.

Any suggestions?


texts:
I downloaded and loaded osm data into ArcGis Catalog using OpenStreetMap plugin.

But there's no field named "NAME" in the attribute table. And therefore I cannot display the name on the map.

Did I do something wrong?

Following are the steps and results in screenshots.


texts:
I would like to aggregate polygons that are within 2 meters of each other, like the Aggregate Polygons (Cartography) tool in ArcMap, but using GDAL with Python.

To add some detail in case the problem can be approached earlier in the pipeline, these polygons were created using Polygonize in GDAL from PNG binary masks with the following code:

dst_layername = blockname + "_Piles"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_layername + ".shp" )
dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )

gdal.Polygonize( srcband, None, dst_layer, -1, [], callback=None )

To demonstrate what I mean using the ArcMap tool, here are my polygons before:

And here are my polygons after running the Aggregate Polygons tool:

I wish to do this using non-ESRI tools that I can call from Python. I have been using osgeo and GDAL for my project thus far.


texts: Is the inclusion of Active Directory Groups in ArcGIS Server User management possible?

texts:
I am currently installing ArcGIS server 10.3.1 to replace 10.0. I would like to base my user management on our organisation's Active Directory user list. It seems that I receive the total list of Active Directory users, but cannot filter according to user groups defined through an RDN Attribute.

Can anyone tell me if they have succesfully filtered users through a group in the Active Directory?


texts: Aligning layers with questionable metadata in SW Ontario

texts:
After making my best efforts to read/try everything I can on this forum and elsewhere, I unfortunately have to contribute to the wide literature of "my layers aren't aligning" questions.

I am working in SW Ontario & have some somewhat reliable reference data in NAD 1983 CSRS (Canadian Spatial Reference System). Then I have two shapefiles of points:

1). One that I created from a raw spreadsheet of lat-longs. The source reports (somewhat half-heartedly) that it should be in NAD 1983, so that's the geographic coordinate system I've defined for it.

2). One that was sent to me already as a shapefile, projected into NAD 1983 Great Lakes Basin Albers. I think it's quite possible this was not the data's "original" coordinate system, but I don't know. I've gotten the original raw data and can try to plot it myself. Ultimately I want everything to be in Great Lakes Basin Albers.

Between these 2 datasets, certain points are supposed to be redundancies. When I hone in on such points, they have the same lat-longs (in the fields that were used to plot them), but the point from #2 is always 13-23 meters NE of the point from #1. I say my reference data (streams) is "somewhat reliable" b/c, while I know it's projected right, the points aren't necessarily going to be right on top of the stream lines, though the points from #1 often are, while #2 is offset. Thus, I've assumed my streams and #1 are okay.

This kind of looks like a central N. America NAD83 / NAD27 offset to me, but when I reproject #2 to a geographic coordinate system of NAD27, it only gets further away to the northeast. Throw in the datum transformations needed when the CSRS version of NAD comes into play, and I'm quite unsure how to diagnose my problem in a systematic way or fix it, i.e. to get these points with identical lat-longs to actually be on top of each other, preferably at the location of #1.

EDIT: The #2 points appear to have an attribute for "Geodetic Datum" specifying either NAD 1927, NAD 1983, or WGS84 (& very often blank). They're museum records covering roughly 1930-2010, so this kind of makes sense. Does this mean their lat/long attributes are all based on different reference points & the shapefile can't be fixed in a uniform way?


texts: Interpolating logarithimic values - groundwater pH

texts:
I am using Spatial Analyst to interpolate surfaces for several groundwater parameters based on point observations at monitoring bores, including groundwater level and pH.

It suddenly occurred to me however, that pH itself is a logarithmic measure of concentration of hydrogen ions. I can't seem to clarify in my head if this requires special attention.

In the most basic conceptual example, if point A is pH=4.0 and point B is pH=6.0, it could be assumed that a point C directly in the middle of the two (equidistant to both) might be ph=5.0 (left-hand mark on graphic attached). If we were looking at actual hydrogen ion concentration as the parameter, the value at point C might be closer to 5.7-5.8 (right-hand mark on the graphic attached).

Could someone please help clarify the assumptions I might be missing here and the correct way to handle this when interpolating?


texts:
I successfully installed QGIS 2.0 on my macbook running OSX Mavericks, but when I attempt to open the program, I get the following message:  "QGIS cannot be opened because of a problem. Check with the developer and make sure QGIS works with this version of OS X. You may need to reinstall the application. Be sure to install any available updates for the application and OS X."

As far as I can tell, everything is updated and I reinstalled QGIS twice with no luck.

Problem details and system configuration:
Process:         QGIS [1087]
Path:            /Applications/QGIS.app/Contents/MacOS/QGIS
Identifier:      org.qgis.qgis2
Version:         ???
Code Type:       X86-64 (Native)
Parent Process:  launchd [132]
Responsible:     QGIS [1087]
User ID:         501

Date/Time:       2014-01-06 17:36:06.120 -0500
OS Version:      Mac OS X 10.9.1 (13B42)
Report Version:  11
Anonymous UUID:  FC49D1E1-1069-09AE-FE37-E5ADF06A6277

Sleep/Wake UUID: 03863BF4-29D6-4A54-A415-DD0D62750EEF

Crashed Thread:  0

Exception Type:  EXC_BREAKPOINT (SIGTRAP)
Exception Codes: 0x0000000000000002, 0x0000000000000000

Application Specific Information:
dyld: launch, loading dependent libraries

Dyld Error Message:
  Library not loaded: /Library/Frameworks/PROJ.framework/Versions/4/PROJ
  Referenced from: /Applications/QGIS.app/Contents/MacOS/QGIS
  Reason: Incompatible library version: QGIS requires version 8.0.0 or later, but PROJ provides version 7.0.0

Binary Images:
    0x7fff6b7f7000 -     0x7fff6b82a817  dyld (239.3)  /usr/lib/dyld
    0x7fff87216000 -     0x7fff87280ff7  com.apple.framework.IOKit (2.0.1 - 907.1.13)  /System/Library/Frameworks/IOKit.framework/Versions/A/IOKit
    0x7fff8f37e000 -     0x7fff8f37efff  com.apple.ApplicationServices (48 - 48) <3E3F01A8-314D-378F-835E-9CC4F8820031> /System/Library/Frameworks/ApplicationServices.framework/Versions/A/ApplicationServices
    0x7fff90e77000 -     0x7fff9105cff7  com.apple.CoreFoundation (6.9 - 855.11)  /System/Library/Frameworks/CoreFoundation.framework/Versions/A/CoreFoundation

Model: MacBook5,1, BootROM MB51.007D.B03, 2 processors, Intel Core 2 Duo, 2.4 GHz, 2 GB, SMC 1.32f8
Graphics: NVIDIA GeForce 9400M, NVIDIA GeForce 9400M, PCI, 256 MB
Memory Module: BANK 0/DIMM0, 1 GB, DDR3, 1067 MHz, 0x80CE, 0x4D34373142323837344448312D4346382020
Memory Module: BANK 0/DIMM1, 1 GB, DDR3, 1067 MHz, 0x80CE, 0x4D34373142323837344448312D4346382020
AirPort: spairport_wireless_card_type_airport_extreme (0x14E4, 0x8D), Broadcom BCM43xx 1.0 (5.106.98.100.22)
Bluetooth: Version 4.2.0f6 12982, 3 services, 23 devices, 1 incoming serial ports
Network Service: AirPort, AirPort, en1
Serial ATA Device: TOSHIBA MK2553GSX, 250.06 GB
Serial ATA Device: MATSHITADVD-R   UJ-868
USB Device: Built-in iSight
USB Device: BRCM2046 Hub
USB Device: Bluetooth USB Host Controller
USB Device: Apple Internal Keyboard / Trackpad
USB Device: IR Receiver
Thunderbolt Bus:


texts: Deactivating layer on loading map using OSM Buildings in Leaflet?

texts:
I'm trying to use OSMBuildings in my Leaflet map.

<!doctype html>

  <meta charset="utf-8">
  <title>OSM Buildings Test</title>
  <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7.2/leaflet.css"/>

All is working but in my L.control.layers, I note that the layer "OSM Buildings" starts activated and I'd like to have it deactivated.

How do I do it?

You can test it


texts: Calculate side of a polygon that touches another feature

texts:
How can I measure the length of the side of a polygon that touches another feature (calculate from A to B on the image below)?

I'm creating and measuring lines to figure out the lenght of the 'intersection' and I wish that was a easier method using a geoprocessing tool.

How can it be done?

I don't have programing skills and I am using ArcGIS 10.2.2.


texts: Ogr: How to add layer description to dgn?

texts:
I have created a dgn-file using Ogr, but can't find a way to fill out the description for this layer programatically. (something like _dgnLayer.Description = "Highways")

opts.Add("SEED=" + seedfile);
opts.Add("COPY_WHOLE_SEED_FILE=YES");

// Prepare
Initialization.PrepareGDAL();

// Initialize output DGN
using (var dgnDriver = Ogr.GetDriverByName("DGN"))
{
    _dgnDataSource = dgnDriver.CreateDataSource(filename, opts.ToArray());
    _dgnLayer = _dgnDataSource.CreateLayer("", null, wkbGeometryType.wkbUnknown, null);
    _dgnFeatureDefinition = _dgnLayer.GetLayerDefn();
    _dgnFeatureDefinition.AddFieldDefn(new FieldDefn("ErrorTest", FieldType.OFTString));
}


texts: Where can I find the viewing properties on ERDAS?

texts:
The instructions state to click on vector and then viewing properties yet I cannot find the viewing properties.

When I click on vector in the main menu there is no option for viewing properties.


texts: Coordinate transformation/reprojection using DotSpatial?

texts:
The following piece of code worked with an older version of DotSpatial.

public static Geometry Project_EPSG25832_To_EPSG3857(byte[] wkb, int SRID)
    {
        NetTopologySuite.IO.WKBReader reader = new NetTopologySuite.IO.WKBReader();
        Geometry geom = (Geometry)reader.Read(wkb);

        var epsg25832 = new DotSpatial.Projections.ProjectionInfo();
        var epsg3857 = new DotSpatial.Projections.ProjectionInfo();
        epsg25832.ParseEsriString(ESRI_EPSG_25832);
        epsg3857.ParseEsriString(ESRI_EPSG_3857);

        Geometry transformedGeom = (Geometry)DotSpatial.Projections.GeometryTransform.TransformGeometry(geom, epsg25832, epsg3857, NetTopologySuite.NtsGeometryServices.Instance.CreateGeometryFactory());

        return transformedGeom;
    }

I had to upgrade the DotSpatial suite to the latest version and now the above code breaks. It seems the "DotSpatial.Projections.GeometryTransform" class doesn't exist any more.

I can't seem to find anything similar in the new version of DotSpatial except for this method

Reproject.ReprojectPoints(xy, z, pStart, pEnd, 0, 1);

...but it takes an array of points...

Anybody know of a nice and easy way of coming from a WKB to an array of points, or maybe there is another library that is easier to use in .NET/C#??