texts: Exporting Layout View as GeoTIFF controlling resolution and size with Python?

texts:
I have two data frames ("Final" and "OverviewMap1") on one Layout View and I can't export them as a whole in one GeoTIFF. What is wrong with the following script:

import arcpy
mxd = arcpy.mapping.MapDocument(r"C:\Users\gnowak\Desktop\Fredropol_FINALv3.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Final", "OverviewMap1")[0]
arcpy.mapping.ExportToTIFF(mxd, r"C:\Users\gnowak\Desktop\Fredropol.tif", df,
                           df_export_width=1346,
                           df_export_height=1003,
                           geoTIFF_tags=True)
del mxd

Edit:

The following script works but I want to export GeoTIFF layout and want to setup resolution, width and height of tif.

import arcpy
mxd = arcpy.mapping.MapDocument(r"C:\Project\Project.mxd")
arcpy.mapping.ExportToPDF(mxd, r"C:\Project\Output\Project.pdf")
del mxd


texts: Why can't I see v.net in QGIS GRASS tools

texts:
I have noticed answers such as Is there an open source tool to planarize graphs? stating you can use GRASS v.net from within QGIS.  However I can't see it - my QGIS processing GRASS toolbox contains 160 algorithms, none of which are v.net.

Am I looking in the wrong place, or has it been removed?


texts: Selecting raster class with highest count and attributing it to field polygon using ArcGIS Desktop?

texts:
I have a 5m satellite image that I classified in ENVI. I brought this classified map into ArcGIS 10x as a .dat image and extracted the raster class for every pixel centre (5m) as a point shapefile. Now I want to intersect the point shapefile with a polygon shapefile: e.g. I have a 160 ac quarter section polygon and have 100's of 5 m pixels within it. My goal is to attribute the field polygon with the dominant raster class.
I am thinking this is a matter of clipping the points within the field polygon, counting the number of raster classes, and selectiing the raster class with the highest count and attributing it to the field polygon.

I am wondering if there is anyway to automate this in ArcGIS?


texts: GeoServer SLD: fill a polygon with color sampled from gradient

texts:
I have a "geometry" and a "rating" column. I want to symbolize the geometries in an SLD (PolygonSymbolizer) but fill it with a color value that is interpolated between BLUE-PURPLE-RED based on the rating value, in range -100 to +100.

So the polygon would have a single fill color, but the color should be dynamic, based on the value.

I tried superimposing two symbolizers with scaled alpha values for the fill color: -100 to 0.0, +100 to 1.0 for red and the negated value for blue, but that doesn't work correctly because I get different results depending on whether I do red first and then blue or the other way around.


texts: libtiff.net how to read elevation data from non-tiled GeoTIFF files C#

texts:
I am using libtiff.net to read elevation data from GeoTIFF files and it works fine with tiled files but I can't obtain the elevation data from either stripped or scanline-based files especially the ones downloaded from NASA.
Any idea how to do that?


texts:
Has anyone any luck subsetting a small part out from NAIP?

I seek your feedback on pros or cons of how to subset an image from NAIP despite they are in mrsid and their uncompressed file size are hefty size.

What software do you use and why do you use it when you subset them?


texts: How to format coordinates as preparation for transformation to WGS1984 in ArcMap 10.8?

texts:
I have a pointshape with a lot of points, containing X and Y values projected in krasovsky_1940_albers. I want to transform the projection to WGS1984 and tried to use the Tool "Convert Coordinate Notation" but it only response with the Warning 001293:coordinate conversion failed.
I found a thread with a similar problem (Transforming coordinate system in ArcMap gives Warning 001293?) where  FSimardGIS suggests to format my coordinates before the conversion. Since I am very new to arcmap, I don't have any idea how to to such things. Is there a tool or a method to use the field calculator for such things? I would like to it not point by point because it would take to much time. I tried:
[POINT_X].SetFormat( [ + | - | E | W ] <DD.dd> [ + | - | E | W ]).AsString
in field calculator, but it did not work.
Here is a Screenshot of my CCN input:

In addition I tried:
Convert a raster to point (change the reference system), export (and inport) the attribute table of the point shape, then set the coordinate system to WGS1984 under DISPLAY X,Y, Export the Events as layer and inport the new Point shape, and it also shows me that the reference system is WGS1984, but when I calculate the coordinates by calculate field geometry in decimal degree, I only get 0 as an outcome.
How can I do this?


texts:
I thought you could give me some advice regarding my problem. I have a working WFS Layer and I want to be able to type in some text which then should be used as the filter value. Afterwards the map should be refreshed, however only with the features selected through my input text. For Example I want to see Restaurants, so I type it in and it refreshes my map.

This is my code so far. In advance - thanks for your help.

     <title>OpenLayers Demo</title>

    <form name="form" id="form">
        <input type="text"  name="Search" id="searchText1" value="" />
    <input type="submit"  />
      </form>


texts:
I am trying to create a map using Leaflet with circles that indicate distances  on the map at current zoom level (map scale). See illustration below:

I know how to calculate sizes and update the circles themselves, but I am struggling to add the text to the circles.

I have tried the leaflet.label plugin but this wasn't an option because it doesn't support static labels for Vectors.

The alternative seems to be to use a divIcon, but this seems to be a very clunky way to solve my problem, as I will have to edit the options of the divIcon, inside a marker on each zoom event.

What would you advise?


texts: Can I save MapInfo thematic*.thm files to a database?

texts: I have an MS-SQl database. Is it possible to save mapinfo *.thm files to a database as I would save mapinfo *.tab files?


texts: Create MultiDiGraph from Shapefile using Networkx

texts:
I have a .shp with a road network that I want to represent in a Graph with NetworkX. I have used it for further analysis such as shortest distance etc. The problem is the function in Networkx converts shapefile to DiGraph only. The Digraphs are not suitable for representing multiple edges to same nodes

Is there any way to create a MultiDiGraph so I don't lose any edges, or a way around it?


texts: Determine what is activated by QAction.activate(QAction.Trigger)

texts:
I'm trying to build a plugin for Qgis where I want to use some functions of a default plugin. I would like to control the options that the user are supose to fill in while using the Interpolation plugin.
So far I have only managed to start the plugin in my python code with:
from PyQt4.QtGui import QAction
for item in qgis.utils.plugins['processing'].iface.rasterMenu().actions():
    if "polation" in item.text():
        item.menu().actions()[0].activate(QAction.Trigger)

Now I want to fill in this form with python code:
But I can't figure out what I activated with item.menu().actions()[0].activate(QAction.Trigger), I've searched through the children of item.menu().actions()[0] but I can't find any with the title "Interpolation plugin" or anything else on the picture :(
I'd tried item.menu().actions()[0].dumpObjectTree() and got nothing in return..


texts: Combining data with partially overlapping polygon layers? (QGIS)

texts:
I am new to QGIS.
I have a dataset with Baltimore neighborhoods, and then a series of datasets for Baltimore census tracts dating back to the 1930s. My goal is to look at how population and the amount of units has varied over the years for each neighborhood. However, neighborhoods and census tracts overlap in imprecise ways so that the following situations occur:

Most neighborhoods contain multiple census tracts.
Some census tracts overlap multiple neighborhoods.
Some census tracts contain more than one complete neighborhood.

I think the best way to handle this will be to aggregate population counts by the weighted ratio of a given census tract area's containment in a given neighborhood. As some examples (see image as well)
In the 1940 census, Reservoir Hill (neighborhood) contains 100% of CT 0132 (pop 5,070), 80% of CT 0131 (pop 4,795), 20% of CT 0133 (pop 6263), and a tiny sliver (let's say 2%) of CT 0134 (pop 3,934). So 5070 + (4795*.8) + (6263*.2) + (3934*.02)=10,237.
Parkview/Woodbrook is an entire neighborhood contained within the CT 0134 in 1940. Since  maybe 10% of 0134 is contained in the neighborhood, pop 3,934*.1=393.
I know this will introduce some aggregation bias, but I think it's likely a lot more accurate than just finding the center point of each census tract and then assigning the total population to whatever neighborhood contains the center point.
Is there a reasonably feasible way to conduct this type of analysis and calculations w/in QGIS?
Here's sort of what I'd envision as an output of quasi-arrays in the "Neighborhood" shapefile:

Neighborhood
1940s_CT
1940s_CT_Full Pop
1940s_CT_Area_Percentage
1940s_CT_Estimated_Pop

Reservoir Hill
G240510000131,G240510000132,G240510000133,G240510000134
4795,5070,6263,3934
.8,1,.2,.02
10237

Parkview
G240510000134
3934
.1
393.4

This is the most related issue I've found so far, though I've never used PostGIS.


texts:
I want to export layers in QGIS as a S-57 basefile, but can't get it to work.

Is there any documentation available regarding the processing parameters?

I get error message when i try to export, regardless of geometry type.
Example:

 Export to vector file failed.
 Error: Creation of layer failed (OGR error: CreateLayer() not
 supported by this dataset.)

I mainly want to know where I can find more documentation about the workings of this export feature.


texts: Linux gives different result from Windows after import with osm2pgsql

texts:
I've done a blank Django project. Then shared it via git on gitlab.
Then downloaded a very small region from geofabrik
Then:

in Windows I've imported it via:
osm2pgsql -c -U xx -d yy -H localhost -W -E4326~/myfile.osm.pbf
in Linux I had to change 4326 to 3857 to make it work, otherwise coordinates in the admin interface didn't show up properly.

From there I've made my models:

Here's are my models (all the same fields for Line, Polygon, Point and Roads, just one is enough):

class PlanetOsmLine(models.Model):
    class Meta:
        db_table = 'planet_osm_line'
        managed = False

    osm_id = models.BigIntegerField(primary_key=True)
    # ...tons of other field (unused)...
    way = GeometryField(default=None)

Now, for the same osm_id, this requests works under Windows, and gives me many results but under Linux I have an empty array (!!):

source = PlanetOsmPolygon.objects.get(osm_id="717295797").way
print([a.osm_id for a in PlanetOsmPolygon.objects.filter(
    way__contains=source
)])

What am I doing wrong, what are the tests I should do to make it work/look for the problem? I dont know where to start.


texts:
I am writing a script that creates a new geodatabase, copies the feature classes from one geodatabase to this new geodatabase, and then writes the name of the feature class, the original geodatabase, the destination geodatabase, the date, and the time it was copied to a text file.
I have all of that working correctly except for the time field. For some reason, the time stamp stays exactly the same despite copying the feature class over at different times. I tried nesting the time stamp within the for loop, but still they are all the same even though they are not copied instantaneously.
If I leave the time stamp above the for loop, it stamps once, but if I nest it within the for loop, it still does the same thing.  Has anyone encountered something like this before? I expected to have at least a few seconds between each copy.

for fc in fcs:
    arcpy.CopyFeatures_management(fc, output_geodat + "\\" + fc)
    fc_cop = str(fc)
    Hour = str(nowutc.hour)
    min = str(nowutc.minute)
    if len(str(min))<2:
       Min = "0"+str(min)
    else:
        Min = str(min)
    sec = str(nowutc.second)
    if len(str(sec))<2:
        Sec = "0"+str(sec)
    else:
        Sec = str(sec)
    Time2 = str((int(Hour)-4)) + ":" + str(Min) + ":" + str(Sec)
    print fc_cop +"," + geod + "," + output_geodat + "," + Date + "," + Time2
#Step 9: Log feature classes and date of copy.
    log("\n" + str(fc_cop) +"," + str(geod) + "," + str(output_geodat) + "," + str(Date) + "," + str(Time2))

print "Finished"


texts: Implementing custom OAuth2 login for ArcGIS Web AppBuilder developer?

texts:
I have a web mapping application on ArcGIS Web AppBuilder developer edition. I downloaded the app and published on IIS.
I want to validate users using OAuth2 like IdentityServer, Google or Facebook.
The proxy.config file contains some settings that used OAuth2.
<serverUrl url="http://route.arcgis.com"
    clientId="6Xo1d-example-9Kn2"
    clientSecret="5a5d50-example-c867b6efcf969bdcc6a2"
    rateLimit="600"
    rateLimitPeriod="60"
    matchAll="true">
</serverUrl>

Can I change these settings for use IdentityServer or Google auth?


texts: QGIS print composer - Legend columns not aligning when title is deleted

texts:
I want a simple legend without a title and with 2 items side by side. I can't find an option to turn the legend off so have just deleted text from the Legend title box. I select Columns (2) - Equal widths - Split layers.

What I get is misaligned entries (see attached).

I've tried changing the Title space to 0 and played around with the other spacing options but cannot get the writing to align.


texts: Looking for succesful implementaion of PPDM in form of geodatabase

texts:
I am have difficulties in implementing spatial version of Professional Petroleum Data Management (PPDM) for our GIS database.

Does anyone have a reference to a sucessful case of implementing PPDM in a geodatabase?

Or how to build a spatial PPDM database?


texts: Create function for filtering date by quarter in Google Earth Engine

texts:
I'm trying to figure out a function to filter image by quarter (January - March, Arpil - July, etc) of Landsat Image Collection to calculate NDVI for each three-month period between 2013 and 2018. My purpose is to adjust the code this link https://code.earthengine.google.com/f590c95bd17454789cde0eb8d1be4850 to compute NDVI in quarter instead of five-year period.

I am planning to add this code below, however i don't know exactly it could be possible or not.

var start = ee.Date('01-01-2013');
var end = ee.Date('31-03-2013');
var date_range = ee.DateRange(start,end);
var fnc = function(d1) {
  var S1 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
    .filterDate(date_range)
    .filterBounds(table)
  return(S1.first())
}
var quarter1 = ee.ImageCollection(dates.map(fnc));
print(mt)


texts: Vector tile selection layer (canvas) drawn outside the map div. Bug?

texts:
I loaded a PBF vector tile layers into my map, and when I select a feature, OpenLayers 6.3.1 creates the selected layer canvas outside the map.
Is it a general problem, or just mine?
I don't attach any personal code, because I encounter the same error on the official site:
https://openlayers.org/en/latest/examples/vector-tile-selection.html


texts: E-mail ArcGIS Online GeoForm submissions

texts:
I have made a Geoform for folks to report incidences; however, it's useless unless someone is notified via e-mail that an incident has been reported on the map.  Does anyone know how to link the submit button to an e-mail address in addition to it's regular function?  I opened up the code but I'm not totally sure what I'm looking at.

Suggestions or is there another way to accomplish the same end goal?


texts: Issue with Weighted Raster Overlay Service for Suitability Modeler Widget

texts:
I'm attempting to build a hosted Weighted Raster Overlay Service (WRO) in ArcGISPro to use in AWAB's Suitability Modeler Widget. I've downloaded the Weighted Raster Overlay Service toolbox via THIS source and prepared raster layers as described in the documentation (same cell size/integer/LZW compression etc). I'm running the "Create Weighted Overlay Mosaic" tool from the toolbox on 10 input GeoTIFFs from a folder, in order to create a Mosaic Dataset. The rasters in TOC are either individual integers representing a given theme, or grouped into classes via "Symbology" tab in ArcGIS Pro. Every time I run the tool it seem to complete with warnings, which seem to prevent me from going through the next steps of sharing the mosaic as a WRO Service. For each layer, I'm getting the following warnings:

"Removed extension tif from LAYER 1,2,3,......"
"Adding rasters to the mosaic
LAYER1.tif is missing layer data.
LAYER1.tif will not be inserted into the mosaic."

This is repeated for all of my input layers. The rasters are present in the Mosaic Dataset, however I'm not able to run either of the two other tools ("Update WRO Layer Classification"/"Update WRO Layer Info") as the tool does not recognise the input mosaic as a valid dataset (Error: UnboundLocalError: local variable 'label_list' referenced before assignment). Consequently I cannot publish the mosaic to a WRO Service.

Has anyone come across this issue before? I was not able to find any reference regarding these errors and hoping someone may be able to offer some advice here.


texts: How to create regional map tiles for Leaflet

texts:
I would like to create the map tiles in leaflet, which will be used in regional scale. I mean, that once I provide sth like this:

  L.tileLayer('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png?{foo}', {foo: 'bar', attribution: 'Map data &copy; OpenStreetMap contributors, CC-BY-SA'}).addTo(map);

I will receive the tile map i.e for 1 or 2 cities covered instead of a whole world (non covered areas will remain as per tilelayer underneath).

I found a few possible solutions:

QGIS Qtiles - https://www.qgistutorials.com/en/docs/creating_basemaps_with_qtiles.html
Slippy Map
OpenMaptiles (with Python)

I am quite new in this issue. So far I found a few solutions mentioned here:

How can I convert an image into map tiles for leafletJS
that led me to the:
How to handle a geotiff map together with points with OpenLayers (or other libraries)?

Generating Tiles with QGIS?

https://leafletjs.com/plugins#non-map-base-layers

https://support.maptiler.com/i25-use-maps-interactively-with-javascript

Is there some quick solution to grasp it?


texts: Formatting SQL expression for where clause in Search layer by attribute with ArcPy?

texts:
I am trying to pass a variable into a where clause in the "Search Layer by Attribute" arcpy function. I cannot seem to figure out the correct format for the SQL expression so that the function reads the value of the given variable, and doesn't just read the variable name as a string.

The variable I am trying to pass is called "currentRaster" and is a list of dates in the 'YYYY-MM-DD 00:00:00' format.

The expression I am having trouble with is:
whereClause = """"DATE" = date 'currentRaster'"""

If I run the python script and tell it to "print whereClause", the result is:
whereClause = """"DATE" = date 'currentRaster'"""

instead of, for example:
whereClause = """"DATE" = date '1998-01-31 00:00:00'"""

If I substitute one of the dates for the variable name, e.g. 1998-01-31 00:00:00, and run the script, the expression works.

Does anyone know the proper way to format this SQL expression so that it reads the variable as the date it contains?

The rest of the script it pasted below.

#Interpolate points from temporal ArcGIS rasters

#import modules
import arcpy
from arcpy import env
import datetime
from datetime import timedelta

#environment settings
env.workspace = "Z:\\sharks\\Models\\Summer2013\\Data_From_Jason\\compiled Chl"
arcpy.env.overwriteOutput = True

#enable MGET modules and enable verbose logging
import GeoEco.Logging
from GeoEco.Logging import Logger
Logger.Initialize(activateArcGISLogging=True)
import GeoEco.SpatialAnalysis.Interpolation
from GeoEco.SpatialAnalysis.Interpolation import Interpolator

#set up variables
rasters = arcpy.ListRasters()
points = "Z:\\sharks\\Models\\Summer2013\\summer2013.gdb\\bll_wgs84"
fields = ["Chla_Daily"]
method = "Automatic"

#define function that extracts the gregorian date from the raster filename (which is in julian days)
def get_date(fname, nbeg, nend):
    yr = int(fname[nbeg:nbeg+4])
    j = int(fname[nbeg+4:nend])
    beg = datetime.datetime(yr,1, 1)+ timedelta(days=(j-1))
    return beg

#create feature layer for select layer by attributes analysis
arcpy.MakeFeatureLayer_management(points,"points_lyr")

#set up loop to step through all rasters
for i in rasters:
    currentRaster = get_date(fname=i, nbeg=16, nend=23)
    #set up cursor to step through all records in attribute table
    with arcpy.da.SearchCursor("points_lyr", ("Date",)) as cursor:
        for row in sorted(cursor):
            value = row[0]
            print value
            print currentRaster
            #if this record's date matches the date of the current raster, then select that record and perform interpolation
            if value == currentRaster:
                whereClause = """"DATE" = date 'currentRaster'"""
                print whereClause
                arcpy.SelectLayerByAttribute_management("points_lyr", "NEW_SELECTION",whereClause)
                Interpolator.InterpolateArcGISRasterValuesAtPoints([i], "points_lyr", fields, method,None, None,False,None,None,None,None)
            else:
                print "no"

arcpy.Delete_management("points_lyr")
del row, cursor


texts:
I have ILayer, how do i obtain INALayer,

if (TypeOf layer Is INALayer) is false?

I am using ArcGIS 9.3.1


texts: How to configure a write only layer on Geoserver?

texts:
This is probably a wrong way of doing it, but I was exploring this option only because I do not know how to implement the right solution.

We have a layer to which features are added using WFS-T. We have configured Geoserver to authenticate and authorize via LDAP.

While querying for features, we would like Geoserver to return features based on the user/role.

Since I do not know how to set feature based security (row level security), my thought was to see if we can make the layer write only and not allow any read operation.

The read will be done through a SQL  Parametric View layer will add a WHERE cause to filter by an unique value.

For doing that, I have this setting in layers.properties
workspace.layer.w=ROLE_USER
workspace.layer.r=SUPERUSER

However this doesn’t seem to work and I am not able to do any WFS-T on the layer although this user has the correct role.  Reading s

What would be the right strategy to implement this? Thanks in advance.


texts:
I've been developing a disaster risk management application that makes use of Java/Smart GWT, PostgreSQL/PostGIS, and Tomcat 7. I'm looking for a shared hosting environment for the beta testing period (about 6 months).

My current hosting provider does not provide PostGIS, or PostgreSQL greater than version 8.1. Any pointers or advice will be appreciated.


texts: Unable to delete feature class from file geodatabase

texts:
I am trying delete feature class which is corrupted from file geodatabase but its throwing error

Failed to delete selected object(s).
The item does not have a definition. The XML being loaded could not be parsed.[physical name:my feature class name]
Xml Parsing failed.
Reason:An invalid character was found in text content.


texts: How to add layers and update layer control dynamically : leaflet

texts:
Introduction

  I am working with the leaflet api, where the user can draw shapes to a
  map(image)...

  Initially the layer control(handling 1 layer) is added for base map
  using imageoverlay......

Problem

  Now, I have added a button of id newLyer to page where click event
  handles the creation of new layer.....i.e user can create new layer
  and update layer control(which is now handling 2 layers)....

  I have used several methods to create the layers and adding to control
  but failed....

Adding new layer to layerGroup

var layerGroup = new L.LayerGroup(),
                    imageOverlayUrl = 'aa.jpg',
                    // New imageoverlay added to the layergroup
                    imageOverlay = new L.ImageOverlay(imageOverlayUrl, bounds).addTo(layerGroup),
                    // New featuregroup added to the layergroup
                    featureGroup = new L.FeatureGroup().addTo(layerGroup);

LayerControl where i needed to add the control(if i am correct)

        var layerControl = new L.control.layers({
            'Main': layerGroup,
            //here i need to add new layer control
            }, null, { collapsed: false }).addTo(map);

OnClick function with so far static code, this will be executed on click

        $('#newLayer').click(function addNewLayer() {
            // Second layergroup not added to the map yet
            var layerGroupNew = new L.LayerGroup(),
                imageOverlayUrlNew = 'bb.png',
                // New imageoverlay added to the second layergroup
                imageOverlayNew = new L.imageOverlay(imageOverlayUrlNew, bounds).addTo(layerGroup2),
                // New featuregroup added to the second layergroup
                featureGroupNew = new L.FeatureGroup().addTo(layerGroupNew);
        });

In Short

  Initially, i have one layer with its control, now onclick function
  creates the new layer which will be added to the map but how i can add
  this layer into layerControl....


texts: Combining DEM Rasters from different sources in ArcGIS Pro

texts:
I am a GIS hobbyist trying to create a merged DEM of the St. Louis region in ArcGIS Pro 2.0. So far I have been successful in creating two mosaics, one of the Illinois side of the Mississippi River and one of the Missouri side. I created both by combining county DEM rasters from two different web sources. Here is an image of the two combined rasters:

Please note that originally these were two separate images that I combined after exporting from ArcGIS. This is not representative of where they are spatially. I would have them as two separate images if not for StackExchange's prohibition of more than two links and images in a question for new users.

The Issue

The problem is what happens when I put the Illinois and Missouri DEMs together in a mosaic:

As you can see, the Illinois side has higher elevation values than the Missouri side. When I compare the values of overlapping areas, the Illinois side is about 300 higher. The trouble is, the value increase in overlapping areas is not consistent. Sometimes it's about 250, other times 375.

I imagine that what I need to do is lower the elevation of the Illinois DEM or raise the elevation of Missouri, but I am not sure how. I have not been able to find anyone who has asked about an issue like this before. I would still like the DEM to be accurate after adjusting the elevation, but absolute accuracy is not essential. My goal is to create an elevation map of the St. Louis region for illustrative purposes.

Data Sources

Illinois

Specifically, I downloaded the derivatives for Jersey, Madison, St. Clair, Monroe, and Randolph Counties. The files were originally in .ovr format.

Missouri:

Missouri had two data sources, both from the same website. First I downloaded DEMS for each county from here. Those were geodatabases. I then had to download individual tile DEMS from here to cover gaps between the borders of each county. Those files differed in format; some were adf and some were geotifs.

Apologies for the awful link formatting. I am not allowed to post more than two links yet.

UPDATE

It turns out there was a unit issue, as explained by the top answer below. I just want to share a guide here for anyone looking at this in the future who has a similar issue with units.


texts: How to determine the visibility of a target at a fixed height above an observer on a given DEM?

texts:
I have a DEM, an observer position and an observer height above ground. Further I have a fixed height above the observer at which a flying target is expected. I now want to determine given the DEM, the observer and the object height all positions on the map where the observer can see the target at that particular fixed height.
I am not sure what kind of algorithm to look at exactly.
I have looked into viewshed algorithms but I can only find algorithms that determine the target height relative to the DEM at each location (not only relative to the observer height).
I also looked at algorithms that determine the visible horizon which is better.
But maybe someone with more experience has a better idea, what I should look into.


texts: Creating a "catchment area" geometry from place points in QGIS

texts:
I have data from Geonames.org on the location of population centres, provided in the format of Center Name, Latitude and Longitude.

Only one latitude and longitude point is provided for each population centre.

I would like to create a layer that turns these points into a shapefile of catchment areas around each Population centre, similar to State or Local Government Area set of boundaries.

This could be done e.g. by drawing a boundary between two points at the midpoint between those points.

Is there a way in QGIS to generate these shapefiles automatically?


texts: Measuring tool shows different measurements for same building in two different map documents?

texts:
So I have two map documents with two different Coordinate Systems set for the layers but using copies of the same shapefile.

Map 1: Coordinate System: NAD_1983_StatePlane_Massachusetts_Mainland_FIPS_2001_Feet
       Linear Unit: Foot_US
       The shapefile in this map has same coordinate system but is in Meters.

Map 2: Coordinate System: WGS_1984_Web_Mercator_Auxiliary_Sphere
       Linear Unit: Meters
       The shapefile's coordinate system in this map document is in the NAD_1983_StatePlane_Massachusetts_Mainland_FIPS_2001_Feet and its linear unit is Foot_US.

GIS Software is ArcGIS 10.2.2

I am measuring a distance of a building with the measure tool. In map 1 it is telling me the wall is 100 ft long. In Map 2 it is telling me that it is 135 ft long.

Which one (if either) is correct and, how do I know?


texts: Modified British System "South Italy Zone" grid - wide residual error georeferencing maps in ArcMap or QGIS

texts:
In accordance with the information of a previous thread on this platform, I am trying the georeference some war maps of Italy within the South Italy Zone grid of the Modified British System. I work with known points of a scanned map but, unfortunately, I get an average residual error of about 300-400 meters, both with QGIS and ArcMap, as well as using the 'coordinate translator' by Thierry Arsicaud’s website. I have never got a similar anomaly working with other map projections, and the results are always been acceptable.

However, They follow the projection parameters of the "South Italy Zone" grid:

Projection - Lambert Conical Orthomorphic
Ellipsoid: Bessel 1841
False Easting : 700000
False Northing : 600000
Central Meridian : 14.0°
Central Parallel : 39.5°
Scale Factor : 0.99906

New Qgis projection parameters:

Projection - Lambert Conical Orthomorphic  →    +proj=lcc
Ellipsoid: Bessel 1841                     →    +ellps=bessel
False Easting : 700000                     →    +x_0=700000
False Northing : 600000                    →    +y_0=600000
Central Meridian : 14.0°                   →    +lon_0=14
Central Parallel : 39.5°                   →    +lat_0=39.5 +lat_1=39.5
Scale Factor : 0.99906                     →    +k_0=0.99906
                       (other proj.4 terms)     +units=m +no_defs

the related Qgis proj.4 working definition:

+proj=lcc +lat_0=39.5 +lat_1=39.5 +lon_0=14 +k_0=0.99906 +x_0=700000 +y_0=600000 +ellps=bessel +units=m +no_defs

According to an informative article about the same subject, I could add the "Italy South" grid projection in ArcMap too. It has the additional complication to need further computation for the transformation between different geodetic systems, rather than QGis.

They follow the parameters which I used in ArcGis for the geodetic system:

Name: Bessel 1841
Angular Unit: Degrees
Prime Meridian: Greenwich
Datum: D_Bessel_1841
Spheroid: Bessel 1841
Semimajor Axis: 6377397.155
Semiminor Axis: 6356078.963
Inverse Flattening: 299.1521535

the parameters to set up the transformation from Bessel to WGS84:

dX = 682
dY = -203
dZ = 480

and from WGS84 to Bessel:

dX = -682
dY = 203
dZ = -480

Despite many attempts, I get the same wide residual error and I couldn't find a solution.

How can I calibrate the parameters in order to compensate the coordinate offset?


texts: Does a 'light' version of global OSM map exist?

texts: I need a global map for my project, but without POIs and other unimportant objects for me. Only big cities and country borders. My project will work in offline network therefore I need a local server with maps and I can't use a global map (40GB) because I haven't enough resources.


texts: Problems using the SAGA in processing framework with qgis 2.6.1 (ubuntu)

texts:
I've got problems using saga in the processing framework within qgis 2.6.1 and SAGA 2.1.4
To solve this problem I was obliged to correct the line 352 of the python file

'/usr/share/qgis/python/plugins/processing/algs/saga/SagaAlgorithm.py'

return 'io_gdal 0 -TRANSFORM 1 -INTERPOL 0 -GRIDS "' + destFilename + '" -FILES "' + source \


texts: Export selected fields in ArcGIS to Excel using Python

texts: I'm am creating a script that will allow a user to add XY coordinates to a point layer, then export (along with any addition user selected fields) to Excel. I am getting stuck on the iteration of the selected fields and how to map them to prepare to export to Excel. I have a parameter that derives the field names from the input point layer and the user can tick each field they wish to include in the Excel export.


texts:
On ArcMap 9.3.1, is there a way to ensure that a toolbar is always enabled?  In my case, I want the Standard toolbar (the one with new/open/save/print icons) to always be visible.  I'd like the toolbar to be enabled on ArcGIS startup, and to prevent the user from disabling it.

Is there a callback that fires when a toolbar is disabled where I can write code to prevent it from disabling, as well as API calls I can make during ArcMap startup to always enable the toolbar?


texts: I'm using the world borders common dataset from CartoDB, which I've then exported so I can grab the_geom column, which CartoDB gives me as a hex value. I'm taking this CSV and uploading it another database so I can have that hex value. But when I try to reimport that data with the hex data, either the file fails to upload or CartoDB is not able to interrupt the hex value. Any ideas why this would be happening?


texts: Rule-based labeling with optional second field

texts:
I need to label some polygons with 2 fields using only one label:

The first being a field in the attribute table named "Bio-Code" the second one being optional if a cell in the field "BBT-Code" is filled with numbers.

The first field and second field should be separated by a slash "/" (not the slash used as operator)

I have yet not been successful at creating an appropriate expression.
Could you point out where I have made mistakes?

Bio_Code
/
CASE
WHEN "BBT_CODE" IS '12500' THEN '12500'
WHEN "BBT_CODE" IS '12714' THEN '12714'
WHEN "BBT_CODE" IS '715212' THEN '715212'
WHEN "BBT_CODE" IS '715222' THEN '715222'
WHEN "BBT_CODE" IS NULL THEN 'None'
END


texts: Read JPEG2000 in R for subsequent random forest classification

texts:
I want to use random forests in R to classify an image. My image is in jp2 format. I read Open JPEG2000 (Sentinel 2) in R and updated my rdal:

rgdal: version: 1.2-5, (SVN revision 648)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15

Still, I see only JPEG not JPER2000 when I run gdalDrivers(). Translation into tiff is not an option for me as my image is ~800 Gb if uncompressed (more than my c drive) and gdal is not implementing random forests.

I also tried
gdal_chooseInstallation('JP2OpenJPEG')
library(rgdal)
myimage <- readGDAL('filename.jp2').

But even after exhaustive googling I can't figure out how gdal_chooseInstallation('JP2OpenJPEG') works. All I am getting is Error in .local(.Object, ...) :

How to I read jp2 image in R without translating it other formats so I can apply random forests?


texts: Where to get free seafloor elevation data

texts:
While I have found multiple sources for continental elevation layers, I have been unable to find sea-floor elevation layers.

Specifically, I'm looking for data for the sea-floor around southern South America.


texts: MapServer / OpenLayers: Is there a limit of 16 zoom levels?

texts:
I set up a MapServer (WMS) with a raster basemap like this (Mapfile):

MAP
  PROJECTION
   "init=epsg:900913"
  END
  WEB
    METADATA
      "wms_title" "mymap"
      "wms_srs" "epsg:900913"
    END
  END
  LAYER
    METADATA
      "wms_title" "myraster"
      "wms_srs" "epsg:900913"
      "wms_enable_request" "*"
    END
    PROJECTION
     "init=epsg:900913"
    END
    NAME         myraster
    DATA         "/path/to/myraster.tif"
    STATUS       ON
    TYPE         RASTER
    PROCESSING   "BANDS=1,2,3"
    OFFSITE      0 0 0
  END
END

Now I display this map with OpenLayers:

var mylayer = new OpenLayers.Layer.WMS( "My Map", "http://127.0.0.1/cgi-bin/mapserv?map=/path/to/mymapfile.map", {layers:'myraster'} );

That works fine so far. However, there seems to be a limit of 16 zoom levels and I would like to be able to zoom in deeper than that. I tried to add the option numZoomLevels:22 but that did not have any effect. My raster image should be large enough - opened in an image viewer the resolution is much higher than displayed with OpenLayers (I also tried using a larger raster image). Is there a way to enhance the maximum number of zoom levels?

EDIT - added OpenLayers code:

var wgs84 = new OpenLayers.Projection("EPSG:4326");
var sphericalMercator = new OpenLayers.Projection("EPSG:900913");

map = new OpenLayers.Map( 'map', { projection: sphericalMercator });

layer = new OpenLayers.Layer.WMS( "OpenLayers WMS",
        "http://127.0.0.1/cgi-bin/mapserv?map=/path/to/mapfile.map",
        { layers: 'myraster' } );

map.addLayer(layer);
map.setCenter(new OpenLayers.LonLat(7.7, 48).transform(wgs84,sphericalMercator), 13);


texts: How to convert x,y to longitude, latitude in Cesium?

texts:
Actually I have divided a line into 90 parts. I know the longitudes and latitudes of end points, and for rest of the intermediate points I know x and y. I have to find height for those 90 points and for that I need those points in proper format of longitude and latitude.

Otherwise, please suggest met another way to find n number of intermediate points between two given points. (latitudes and longitudes)


texts:
How would I change the following code (outJPEG section) if I'd like to have the original attribute table name (eg Adress) and the attribute field name (eg ZipCode) as file name? (Adress_ZipCode)

The code as I have it now would just name the output JPG after the field in the attribute table (ZipCode)

Maybe os.path.split is an option but I am not so familiar with python and cannot make it work..

# export as jpeg

outJPEG = os.path.join(outPath, attr + ".jpg")  #atrr and outPath are defined in another part of the code
arcpy.mapping.ExportToJPEG(mxd, outJPEG, "PAGE_LAYOUT", 200, 200, 300, False, "24-BIT_TRUE_COLOR", 90, False)


texts: Unique constraint on Geometry to resolve INSERT conflicts

texts:
I have a simple table with a geometry column like so:
CREATE TABLE geomtable
(
  gid serial,
  geom geometry,
  trip_count integer
);

I want a unique constraint on the geometry column, such that if the same geometry is being inserted again then, the ON CONFLICT clause should add the trip_count to the existing row.
So my insert query looks something like this:
INSERT INTO geomtable (geom, trip_count) values (<geom>, 123)
ON CONFLICT ON CONSTRAINT <equality_constraint>
DO UPDATE SET trip_count = geomtable.trip_count + EXCLUDED.trip_count;

I tried the approach to enforcing equality constraint as described here like so:
CREATE FUNCTION equality_constraint_func(
  id INT,
  gm GEOMETRY
)
RETURNS boolean AS
  $$
    SELECT NOT EXISTS (
      SELECT 1
      FROM   geomtable AS a
      WHERE  a.gid <> id
        AND  a.geom && gm
        AND  ST_Equals(a.geom, gm)
    );
  $$
  LANGUAGE sql
;

And adding a CHECK constraint like so :
ALTER TABLE geomtable
  ADD CONSTRAINT equality_constraint
  CHECK (equality_constraint_func(gid, geom))
;

However, this does not help me perform the DO UPDATE part of conflict resolution. It just raises an error on duplicate geometry. I of course can't create a primary key out of my geom column as I get the index row requires 42632 bytes, maximum size is 8191 error. I get a similar error on creating a unique index on the geom column. I tried creating a unique constraint using the above syntax by replacing CHECK with UNIQUE, but I get a syntax error. Is unique constraint using the function the way to go here, or do I need to consider some other approach like a trigger, or exclusion constraint?
My version information is as under:
PostgreSQL 12.4 (Ubuntu 12.4-1.pgdg20.04+1) on x86_64-pc-linux-gnu, compiled by gcc (Ubuntu 9.3.0-10ubuntu2) 9.3.0, 64-bit POSTGIS="3.0.1 ec2a9aa" [EXTENSION] PGSQL="120" GEOS="3.8.0-CAPI-1.13.1 " PROJ="6.3.1" LIBXML="2.9.10" LIBJSON="0.13.1" LIBPROTOBUF="1.3.3" WAGYU="0.4.3 (Internal)"


texts: Points within Polygons given certain conditions

texts:
I am trying to locate all the points within a set of polygons given a certain condition on the polygon in MongoDB using Python 3x.

This question differs from the other ones because i) it requires multiple points to be returned and has MULTIPLE polygons to search through, and ii) the polygons are cross-collection search.

Essentially, this is something along the lines that I need, where $geometry: is replaced by all the polygons it needs to search through in my other collection while applying the condition 'class' : 'B'.

db.locations.find({'geo_coordinate' :
                    { $geoWithin :
                    { $geometry :
                        db.polygons.geometry: {'class': 'B'}

                } } } } )

Locations structure (21k documents):

{
  "la": 43.57386,
  "lo": -79.5736,
  "address": "123 Stack Dr.",
  "geo_coordinate": {
    "type": "Point",
    "coordinates": [-79.5736, 43.57386]
  }
}

Polygon Structure (w/ coordinates removed | 16k documents):

{
  "type": "Feature",
  "properties": {  },
  "geometry": {
    "type": "MultiPolygon",
    "coordinates": [
    ]
  },
  "class": "A",

}


texts: Mapping multiple circle markers in a single layer with ipyleaflet and jupyter

texts:
How do I show multiple circle markers with ipyleaflet in jupyter?  I want to add them all as one layer (not individual layers) and I do not want to use the MarkerCluster class.  The documenation only references adding a single marker per layer - I tried: (using locations instead of location as referenced by the docs)
marker = CircleMarker(locations=[(172.7, -42.8),(172.7,-42.7)],
                          radius = 10,
                          color = "red",
                          fill_color = "red"
                         )

But that just got me one red circle in the middle of the map.  I can create a whole a list of circle makers and put them into a MarkerCluster but I don't really want to use this feature.  I just want small circles/points at the co-ordinates, in one layer.
df1 = map_df.set_geometry([Point(feature.coords[:][0][:2]) for feature in map_df['geometry']])
marker_list = []
for feature in df1['geometry']:
    marker = CircleMarker(location=(feature.coords[0][1],feature.coords[0][0]),
                         radius = 1,
                          color = "red",
                          fill_color = "red"
                         )
    marker_list.append(marker)

marker_cluster = MarkerCluster(
    markers=(marker_list),
    name = 'first point'
)

I also tried plotting the geodataframe but I could not see how to specify a cirle marker?  It uses the default blue pin which is not useful when I have a number of points to map.


texts: ArcPy Multiprocessing 000714 - Working in Debug Mode

texts:
In an attempt to render an interactive window via matplotlib I am using an out of process script in a toolbox.  If I debug the script, it is possible to step through and render the window as expected.  Running it through ArcMap (not in debug mode) causes a 000714 error.

Arc10.2
Running the script out of process
Tested with background processing on and off (Geoprocessing Options)
pythonExe is being properly found

Any thoughts?

import multiprocessing
import sys, os
import arcpy

import matplotlib.pyplot as plt

def show_form():
    fig = plt.figure()
    x = range(10)
    y = range(10)
    plt.plot(x,y)
    plt.show()

def main():
    pythonExe = os.path.join(sys.exec_prefix, 'python.exe')
    multiprocessing.set_executable(pythonExe)
    multiprocessing.freeze_support = True
    jobs = []
    pool = multiprocessing.Pool(1)
    jobs.append(pool.apply_async(show_form,[]))
    pool.close()
    pool.join()
    del pool
    del jobs
    arcpy.SetParameterAsText(0, True)

if __name__ == "__main__":
    main()

Update 1

As per Luke's comment, I tested cmd.exe /C test.py.  Other than the long import time for arcpy the script executed as expected with the Tk plot window appearing.

Update 2

By opening a command window and quickly snapping a screen capture, I can grab an error from the script before it dies.  This looks like an issue loading the libraries (32 vs. 64bit?)

Update 3

Arcmap is using the 64bit python by default when running out of process.  The next step is to figure out how to disable this, as it appears to be causing a Tk issue as @Luke suggested.


texts: How to mask raster image (DEM) on QGIS Grass

texts: I am doing research at graduate school. Using QGIS GRASS, I made a raster image (DEM) by spline interpolation from point data. Next, I decided to limit the range to within the vector (within the point data). When I masked it in the raster menu, it was not restricted in the vector, and the resulting raster image became a yellow image (with no value). What is the solution for this?


texts: PyQGIS Output leads to invisible layers that don't appear in the layer panel or the canvas but can be accessed in Vector tools

texts:
I'm working on a Processing Tool in PyQGIS and am struggling with the problem described in the title.
Layers are not visible in the "Layers Panel" or on the canvas but can be chosen and worked with in tools like Buffer or Dissolve.
I have 2 variations that produce the same result.
One looks like this:
        poolLayer = processing.run('native:dissolve', {
            'INPUT': outputs['ReprojectPoolkarte']['OUTPUT'],
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }, context=context, feedback=feedback)['OUTPUT']
        QgsProject.instance().addMapLayer(poolLayer)

The other case looks like this:
        alg_params = {
            'FIELD_LENGTH': 10,
            'FIELD_NAME': 'Flaeche_Pool',
            'FIELD_PRECISION': 2,
            'FIELD_TYPE': 0,
            'FORMULA': 'Flaeche_A/poolFl',
            'INPUT': outputs['AttributeNachFeldwert']['OUTPUT'],
            'NEW_FIELD': True,
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['Flaeche_PoolFieldCalculator'] = processing.run('qgis:fieldcalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['ErgebnisNeanteileFinal'] = outputs['Flaeche_PoolFieldCalculator']['OUTPUT']

        layer = context.takeResultLayer(outputs['AttributeNachFeldwert']['OUTPUT'])
        transform_context = QgsCoordinateTransformContext()
        save_options = QgsVectorFileWriter.SaveVectorOptions()
        save_options.driverName = "CSV"
        save_options.overrideGeometryType = QgsWkbTypes.Point
        save_options.fileEncoding = "System"
        save_options.actionOnExistingFile = QgsVectorFileWriter.CreateOrOverwriteFile
        QgsVectorFileWriter.writeAsVectorFormatV2(layer, r'c:\temp\output.csv', transform_context, save_options)

        return results

This actually works fine until I add the part of the code for exporting the layer as a CSV-File as well.
In both case the resulting layers are fine but I have to make them "visible" by manually running them through a buffer with 0m or something similar.


texts: Is points to route still available in QGIS 3.4?

texts:
I need to do the same thing as the post here as I need to create polylines based on given points with sequence attribute, and networks, but I cannot find the tool in toolbox:

my point data and network data

Is function points to route still available in QGIS 3.4?


texts: PyQGIS replace text label in Layout view by new text label from csv tabel

texts:
Currently, I'm stuck working on my first python plugin and hope to get some adivce here. The plugin should take a csv file (or .xlsx would be even better) which has two string columns. Column 1 has terms like "Project Name" or "Project Manager" and the second column has the actual names of the project and the manager. The terms of the first column already exist as text labels in the layout view of my QGIS Project. Thus, the plugin should take the table and replace the text labels in the layout view with the actual names in collum 2. I hope the purpose of the plugin is clear.

I used the "plugin builder" tool to get the basics and succesfully prepared the gui. After that it gets complicated

That's my code so far:

def select_table(self):
    """the name of the input table will be shown in the GUI"""
    filename = str(QFileDialog.getOpenFileName(self.dlg, "Select input Table",
                                           "", '*.csv, *.xlsx, * .xls, * .ods'))
    self.dlg.table_path.setText(filename)

    """This is supposed to add the table as a layer in QGIS"""
    uri = filename
    table = QgsVectorLayer(uri, "table", "ogr")
    QgsMapLayerRegistry.instance().addMapLayer(table)

def getTable(self):
    """gets table from the QGIS Project"""
    layer = None
    tablename = self.dlg.table_path.currentText()
    for lyr in QgsProject.instance().mapLayers.values():
        if lyr.name() == tablename:
            layer = lyr
            break
    return layer

def setVariable(self):
    """adds the table to a variable"""
    self.inTable = self.getTable()

def change_label(self):
    """gets the fields of the csv file"""
    fields = self.inTable.fields()

    """checks if the label texts are equal to those in the first column
    of the csv file. After replaces the label text with the value
    of the second column of the same object"""

    for text in QgsLayoutItemLabel.currentText():
        if text.name() == fields.value.firstCollum()
            text.name() = fields.value.secondCollum

I modified several tutorials to come so far but I guess there are some issues in every function I wrote. However especially tricky is the "change_label" function. It would be great if someone has some advice on that.


texts: ArcGIS - Merge error after Iterate Feature Selection in model builder

texts:
I would like to use the model builder to solve OD Cost Matrix in iterations because my origins/destinations FC need to be grouped by attributes. I use the Iterate Feature Selection tool to do so. After performing all OD Cost Matrix in loop, I want to merge all result lines into a single FC and export its table into an excel sheet.

I separate the Iteration part of the process into a subModel with a Collect Values as output:

Then I use the subModel output into another model to do the merge and the export to excel:

The subModel seems to run in iteration correctly, but when it's time to merge the output I got these errors:

      ERROR 000622: Failed to execute (Merge). Parameters are not valid.
      ERROR 000400: duplicate input not allowed

Can someone have an idea what is  incorrect in my process?


texts:
I have a set of (lat,lon) coordinates that compose a polygon.

Is there a straightforward method to get an analogous set of coordinates defining a smaller version of the polygon?

For example, given the following picture:

If I am provided with the set of coordinates defining the blue shape, how can I obtain the set of coordinates corresponding to the green one?


texts: Why are there strange horizontal lines in raster data?

texts:
I have a feature class (stored in a file geodatabase) containing information on land use for the Amazon rainforest. It is very detailed and looks like this

I need to convert this shapefile to a raster to combine it with other rasters using the Combine tool (I've encountered too many problems using spatial join with such detailed data). The conversion to raster yields the following image

Any idea why the raster image have horizontal lines on them? The horizontal lines seem to contain information (i.e. they are not null values)


texts: What does it mean "GML is an XML grammar"?

texts:
I'm trying to get my head around the theoretical side of XML and GML. The OGC web page for GML states that:

  The Geography Markup Language (GML) is an XML grammar for expressing
  geographical features.

But what does it mean by a "grammar"? I can't see that phrase (as a noun) used anywhere else. Is it a unique GML meaning?

I see the terms "language", "schema", "format" and "standard". Are they the same thing?

So can I say that GML is an XML schema? Or would that mean something different? If so, is there a single XSD file that defines the entire GML specification?


texts:
Has someone an idea how to get dynamic picklist widgets in the object window! Eg. You choose a country in a picklist so the next picklist only shows cities in that country! I know that you could program your own widget but i can not find any example.
Maybe there is even a simpler solution or workaround! I would be happy for any help!


texts: Multiple entries in attribute table for multiple buffers using ArcGIS for Desktop?

texts:
I am using ArcGIS 10.2.2 to draw 500ft buffers around a group of 133 bars.

The numerous buffer tools all produce a single entry in the attribute table for the collected group of buffers instead of one entry for each buffer when I use the dissolve setting.

Is there a way to produce a shapefile of dissolved buffers whose attribute table contains 133 rows, that it, one for each buffer, instead of the default that produces one row for all?


texts: Updating records based on value in another field using arcpy.da.UpdateCursor?

texts:
I am trying to write file names and paths to two separate fields in a feature class, and the value for each record is dependent on finding a match in a list. When I run my code it writes the same values for each record, so something is going wrong. I think its finidng the first match and writing that value over and over, but I want it to look at the value in the field my search cursor is set to, and then find the match in the list and write that out to the two fields. What am I missing?

import arcpy, os
photos_path = r'P:\Records\GIS\Projects\D04_OHS\OverheadSignStructurePics'
OHS_FC = r'D:\Records\GIS\Projects\D04_OHS\D04_OHS.gdb\D04_OHS'
ID_fld = 'OSMI_ID'
file_nm_fld = 'LINK_FILENAME'
path_fld ='LINK_PATH'
photos_list = []

for dirpath, dirnames, filenames in os.walk(photos_path):
    for filename in filenames:
        fullPath = dirpath+'\\'+filename
        if fullPath not in photos_list:
            photos_list.append(fullPath)

arcpy.AddField_management(OHS_FC,file_nm_fld,'TEXT',100)
arcpy.AddField_management(OHS_FC,path_fld,'TEXT',300)
srch_cursor = arcpy.da.SearchCursor(OHS_FC,ID_fld)
updt_cursor_fn = arcpy.da.UpdateCursor(OHS_FC,file_nm_fld)
updt_cursor_fp = arcpy.da.UpdateCursor(OHS_FC,path_fld)
for row in srch_cursor:
    val = str(row[0])
    for p in photos_list:
        f_name, f_path = os.path.split(p)
        if val in p:
            #print val
             for u_row1 in updt_cursor_fn:
                u_row1[0] = f_name
                updt_cursor_fn.updateRow(u_row1)
             for u_row2 in updt_cursor_fp:
                u_row2[0] = f_path
                updt_cursor_fp.updateRow(u_row2)


texts: How to handle multipoint symbology catalogue (as in mil2525c, tactical graphics appendix b)

texts:
Is there a general format (standards based) to store a vector based drawing description for multipoint symbology (means having several anchor points reflecting symbol point <-> geographic point association)?

OGC defined the Styled Layer Descriptor which uses symbolizers (point/line/polygon) but one problem is that they are better suited to be created for a single point representation (to place e.g. a raster based symbol at specific location, maybe with offset) or the have a stronger dependency to the geometry (following a line with an offset, placing somthing at the line, ...) like a line/polygon symbolizer.

If you look at the mil2525c standard there are several symbols where the graphical representation has nothing in common with the stored geometry except the anchor point association to "stretch" the symbol.

So I expect those problems using OGC conformant aproach using sld:

I don't see any default "multipoint" symbolizer, only a point symbolizer, a line symbolizer and a polygon symbolizer. Line and polygon symbolizers both don't fit because they only styles the line/polygon geometry ( => abstaction to mil2525c multipoint geometries and styling them using sld would be very complex... not to say impossible.)
Point symbolizer can call a mark factory (e.g. see GeoTools interface/implementation docs at http://docs.geotools.org/latest/javadocs/org/geotools/renderer/style/MarkFactory.html) to get a shape instead of a raster based symbol but there is no general approach to handle anchor points in a default manner.
If using the 2nd. approach using a point symbolizer: The problem of storing vector based data for symbology generation is not solved. ("Hand-made" generation by building the geometry seems to be an option, but is difficult to analyse or edit. What type is shape in the .net world or what would be the "universal" equivalent (System.Windows.Shapes.Path)?

Are there other standards for multipoint symbols I didn't noticed? How would you solve this problem?


texts: Identifying individual pine trees without using LiDAR data?

texts:
Are there any ways to identify individual pine trees without the use of LiDAR data? Crown identification seems irrelevant.

My goal is to do a rough treetop calculation to get an idea of how many trees are there and how they are layouted.

The segmentation function only classifies trees and the ground. Any other ways to get this calculation done?


texts: Leaflet Search for GeoJSON creating duplicate, non-functional buttons

texts:
I added a plugin for searching a GeoJSON file in leaflet, but now there are multiple search buttons and they do not recognize any input. I suspect that it is initializing for every single instance of a feature within the GeoJSON file. How do I change my script to only initialize once? (continued from this question: Activating Search Plugin with JS Leaflet

  // load GeoJSON from an external file
$.getJSON("Syriashape.json", function (data) {
// add GeoJSON layer and popups to the map once the file is loaded
var syriaLayer = L.geoJson(data, {
    style: function (feature) {
        return {color: 'red', weight: '1pt'};
    },
    onEachFeature: function (feature, layer) {
        layer.bindPopup(

          "Sheet Number: " + feature.properties.Sheet_Numb + "
" +

          "Catalog Record: " + "" + feature.properties.LOC_Catalo + "");

    //activate search button
    var searchControl = new L.Control.Search({layer: syriaLayer, propertyName: 'Sheet_Numb', circleLocation:false});

          searchControl.on('search_locationfound', function(e) {

        e.layer.setStyle({});
        if(e.layer._popup)
          e.layer.openPopup();

      }).on('search_collapsed', function(e) {

        featuresLayer.eachLayer(function(layer) { //restore feature color
          featuresLayer.resetStyle(layer);
        });
      });

      map.addControl( searchControl );  //inizialize search control

    }

}).addTo(map);
});


texts: Is there any Silverlight based WebGIS API like OpenLayers?


texts: Large map of Africa, tiled printing on A3

texts:
I would like to print a large map of Africa on several A3-pages. I want to hang the map on the wall. So the different pages should overlap a minimum.

The map will be used as a road map, so it should have the highest resolution possible, but consist of 5x7 A3-pages. (A3 is the largest format my printer can handle.) How do I print my map into tiled A3 pages with little overlap?

I am currently using QGIS and OSM-data, but other solutions are more than welcome.

Could anyone point me in the right direction?


texts: Reverse geocoding for multiple coordinates

texts:
My code is generating coordinates (Lat, Lon) of 1000 houses. Now I want to reverse these 1000 coordinates to get full address. The code is:

import pprint
from arcgis.geocoding import reverse_geocode
Houses = [0]*(1000)
Houses[0], Houses[1] = (143.5689855, -38.328956999999996), (143.5692555, -38.328993)
for i in range(2, 1000):
    latitude_diff = Houses[i-1][0] - Houses[i-2][0]
    longitude_diff= Houses[i-1][1] - Houses[i-2][1]
    temp = (Houses[i-1][0]+latitude_diff, Houses[i-1][1]+longitude_diff)

    Houses[i] = temp
pprint.pprint(Houses)

I tried reverse_code from ArcGIS but it only converts single coordinate like below:

results = reverse_geocode([143.5689855, -38.328956999999996])

How can I reverse the coordinates of 1000 houses into full addresses?


texts:
I've got geography in sqlserver that is invalid.  Running IsValidDetailed() I get the reason as

24412: Not valid because the interior of a polygon with rings (2) and (1) is not connected.

But I don't know what that means.  my polygon does have 3 holes in it, is it to do with ring orientation?


texts: Replicating Summary Statistics tool with case field using own calculation within ModelBuilder?

texts:
I thought this would be a simple problem to solve, by my googling hasn't provided the answers so far.

I have created a near table (Using the 'Generate Near Table' in Arc) comparing two sets of points, and I have included the option to get the 10 closest points for each of the input points. So for each of the input points I have 10 near features(this can be see in the IN_FID column).

I wish to calculate the average distance, plus the average angle (which I have converted to a bearing 0-360) of the 10 points for each unique value in IN_FID. I use the summary statistics tool with the case field set to IN_FID to get the average distance, however to get the average bearing I need to run a calculation that takes into account circular statistics.

My solution thus far within ModelBuilder, is to use the 'Row Selection' iterator with the IN_FID set as the 'Group by Field option', I then run the calculate field on the output for my average bearing. So for each unique value in the IN_FID column, I have one average bearing (from 10 points). However, this process is very slow, it takes about 3 seconds for each unique IN_FID, and with 10,000 points, is getting on for 8 hours of processing.

The summary statistics tool works in a matter of seconds, so I was wondering if there is a python script I can write and run in the calculate field tool code block, that means I don't have to use the 'Row Selection' iterator and speed up this process? (however, I'm not even sure this will speed it up!). My python isn't up to much at the moment, but my research has lead me to think I need to use search or update cursor, and some (nested?) for loops, but I don't have the knowledge to apply and write this yet, but I am learning using various examples I find online.

Does anyone have any python script examples that do a similar task that I could use as a basis for this task or, can point me in the direction of some useful tutorials?

Or do you have some alternative methods to overcoming this problem that I haven't thought of?


texts: Converting lines to polygon based on matching attribute

texts:
In QGIS, is there a Plugin or a geoalgorithm that allows you to convert lines to polygons but incorporates a particular attribute?
I have roadside kerbs which have a field for road name, and if possible would like to convert these to polygons using the road name so we have the road carriageway, so it would also need to close the polylines.
The conversion tools I've tried don't allow for attribute selection so it means it only creates polygons where there are joining lines, such as a whole block.


texts: LASlib internal error - reading las file with lidR

texts:
Should I assume that the following error occurs because the file is corrupt? Is there any way I can get more information out of this file?
bad_las <- readLAS("~/Downloads/3102_309374.las")
#> ERROR: reading 151140850 bytes of data into header.evlrs[0].data
#> ERROR: cannot open lasreaderlas with file name '/Users/lucasj/Downloads/3102_309374.las'
#> Error: LASlib internal error. See message above.
#> ERROR: reading 151140850 bytes of data into header.evlrs[0].data
#> ERROR: cannot open lasreaderlas with file name '/Users/lucasj/Downloads/3102_309374.las'
#> Error in C_reader(ifiles, ofile, select, filter, filter_wkt) :
#>   LASlib internal error. See message above.
packageVersion("lidR")
#> [1] ‘3.0.5’
packageVersion("rlas")
#> [1] ‘1.3.7’

UPDATE
lasinfo and las2txt produce the following error:
ERROR: reading 151140850 bytes of data into header.evlrs[0].data
ERROR: cannot open lasreaderlas with file name '3102_309374.las'
ERROR: could not open lasreader


texts: Concatenate from all the elements from a list of lists to have a string

texts:
I have data from a JSONs like as below
data={'type': 'Feature', 'geometry': {'type': 'LineString', 'coordinates':
[[410910.0885459831, 2479395.4132211693, 0.0], [410903.8088855865,
2479401.860380645, 0.0], [410845.2521416071, 2479436.215147523, 0.0],
[410844.8989142373, 2479436.568374893, 0.0], [410844.5338823536,
2479436.90979775, 0.0], [410844.1588620354, 2479437.2403241317, 0.0],
[410843.77385328256, 2479437.5581379607, 0.0]]}}

and I wish to make a concatenation on the elements like below in this line:
datas [i] ['geometry'] ['coordinates'] [j] [0]) +" "+ str (datas [i ] ['geometry'] ['coordinates'] [j] [1]of all existing sublists in coordinates.
Here is my code but unfortunately it does not work since it only gives me the last data [i] and not all the coordinates in the sublists
 path = os.path.join(dir, file)
 shapefile = ogr.Open(path)
 shape = shapefile.GetLayer()
 datas= []
 features=[]
 jsonstring=[]
 lines = []
 for i in range(len(shape)):
     features.append(shape.GetFeature(i))
     jsonstring.append(features[i].ExportToJson())   # myjson
     datas.append((json.loads(jsonstring[i])))

     for j in range(len(datas[i]['geometry']['coordinates'])-1):
         line = datas[i]['geometry']['type'] + "("+ str(datas[i]['geometry']['coordinates'][j][0])+" "+ str(datas[i]['geometry']['coordinates'][j][1])+")"
         lines.append(line)`

I'm using Python 3.x, and I would like to have in the end something like this:
LineString(410910.0885459831 2479395.4132211693, 410903.8088855865 2479401.860380645, ..... ,410843.77385328256, 2479437.558137960)

for each data[i]


texts: Where can I get a comprehensive list of Rules expression operators?

texts:
I'm selecting and styling my points layer on a contour map of Philadelphia using Rules Based Styling. (see: https://gis.stackexchange.com/questions/53346/can-i-use-rules-based-rendering-in-qgis-to-select-features).  The filtering process using [string_field = 'string'] turns out not to be accurate for one of my three rules.  I have 'Homocide - Justifiable'  and it returns points with (at least) 'Thefts'.

I'd like to try another expression.  In any case it would be nice to have a list of all the possible expressions.


texts: Arc Hydro Stream WSE to point WSE measurement tool

texts:
I wish to use stream wse to point measurement tool for finding out the flood inundation extent.

Can anyone provide an example of the input time series table to be used in the tool?


texts: Saving output in a gdb Feature dataset in ArcMap Modelbuilder

texts:
I'm using ArcMap 10.5.

I've got a model with a bunch of processes related to different topics such as ROADS, URBAN SQUARES, PARKINGS, etc. I'm currently working from a geodatabase that stores all the results in the same place and sometimes it's hard to find the outputs of a particular topic.

What I'd like to do is to store each output in is corresponding Feature dataset, meaning roads geoprocessing outputs will go to the ROADS Feature dataset and the same with the rest.

I thought that'd be something easy but turns out that when I set the save path to a Feature dataset, that process in particular returns this error:

This works:

C:\geodb.gdb\INDEX_FULLA

Whereas this does not work:

C:\geodb.gdb\RESULTS\INDEX_FULLA

This occurs only with raster outputs. I can't select the feature dataset when navigating and if I write it manually then the error is returned. I have checked the CRS and they are the same.

Also add that I'm trying to save in an empty Feature dataset, it didnd't even have any previous information.

Is there any way to be able to save raster outputs in Feature datasets?


texts:
I am using QGIS 2.6 (Brighton) and when doing a distance matrix, the output csv file appears blank. Is it possible that somehow the shapefile I am using is corrupted?Has anyone please got any hints or ideas to solve this issue?

The shapefile seems to be working well when using the attribute table! (i.e.: selecting points based on the existing fields, etc).


texts: I created a regular grid in QGIS and want to make an atlas that moves from grid to grid. Each grid has a unique id. Under atlas configuaration I chose my grid layer as coverage layer and page name is the id. As you can see in my overview map, my map extent is to high. The right and left side of the extend is perfect. I don't know how to adjust my map extent that it represents exactly my grids. Is that possible?


texts: Using SQL query in QGIS Database Manager?

texts:
I have below query which I am supposed to import into QGIS by amending certain aspects, to find out exactly what it does. But when I try to run it on the corresponding PostGIS layer, DB Manager warns me that there is not such "geom" attribute. I don't know what to do then. Is it a syntax issue?

UPDATE c_risque.ppr_inondation_surf
SET    geom = st_multi(st_simplify(st_multi(st_collectionextract(st_forcecollection(st_makevalid(st_snaptogrid(geom, 1))),3)),0))
WHERE  not st_equals(geom, st_multi(st_simplify(st_multi(st_collectionextract(st_forcecollection(st_makevalid(st_snaptogrid(geom, 1))),3)),0)));


texts: Can a raster map in GRASS be saved in JPEG/BMP format?

texts: I'm trying to save a raster map in GRASS in jpeg/bmp format without using the GUI for d.out.file/d.save. Is there a script to do it? I tried r.out.gdal but it produced a map that is black and white.


texts: Plotting alignment issues with plotting raster list of raster objects and multiple polygons together

texts:
I am very new to R and am trying to plot one object from a raster list with multiple polygons (boundary shapefiles) together in one plot. When I resize the plot window, however, the alignment of the layers change. I am assuming that their alignment needs to be fixed somehow, but I don't know how to do this.

I've included all code for importing raster data (NSIDC sea ice concentration rasters) and shapefiles (one is a management boundary and the other is a digitization of land masses). I want to remove (mask) all land pixels and pixels outside of the management boundary so that I can focus on only the ice/water pixels within the management boundary.

# set up raster -------------------------------------------
pixel <- 25000 #pixel dimension in meters for both x and y
xMin <- -3837500 #From NSIDC: ulxmap -3837500
xMax <- xMin + (pixel*304)
yMax <- 5837500 #From NSIDC: ulymap 5837500
yMin <- yMax - (pixel*448)
r <- raster(nrow=448, ncol=304, xmn=xMin, xmx=xMax, ymn=yMin, ymx=yMax)
projection(r) <- '+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273     +b=6356889.449 +units=m +no_defs'

# bin converting fx
bin_raster <- function(ymd){
filename <- paste("data/","nt_", ymd,"_v1.1_n.bin", sep="")
con <- file(filename, 'rb')
x <- readBin(con, "raw", 300)
x <- readBin(con,"int", size=1, signed=FALSE, 150000)
close(con)
x[x==251] <- 1
x[x>251] <- 0 # set all non-ice pixels to 0% ice cover
x[x<251] <- x[x<251]/250 #scale ice cover to 1
rr <- setValues(r, x) # place result in raster

#change file names to just the date
geoTiffFile<- paste('geoTiffs/', ymd,'.tif',sep='')
rr <- as.list(rr)
}

# create raster list

# extract file names for loop
files <- list.files("data", pattern=".bin$")

# string manipulation
dates <- str_sub(files, 4, -12)

# create an empty list
raster_list <- rep(NA, length(dates))

# loop over all the files and fill the list
for(i in seq_along(dates)){
raster_list[i] <-   bin_raster(dates[[i]])
}

# rename rasters with dates
list_names <- str_sub(dates, 0, -5)
names(raster_list) <- list_names

# Clip to the management area shapefile -----------------------------
# Import shapefile
boundary <- readOGR("filepath/Management_area.shp")
plot(boundary)
# Change projection of shapefile to match raster list (i.e. polar stereographic)
proj4string(boundary)
boundary_proj <- spTransform(boundary, crs('+proj=stere +lat_0=90 +lat_ts=70    +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs'))
proj4string(boundary_proj)
plot(boundary_proj)

# Crop and mask list to the shapefile
raster_list <- lapply(X=raster_list, FUN=crop, y=boundary_proj)

# mask pixels outside of management area
for(i in 1:length(raster_list)){
raster_list[[i]] <- mask(raster_list[[i]], boundary_proj)
}

# look at resulting raster
plot(raster_list[[4]])
plot(boundary_proj, add = TRUE)

# Clip to land ---------------------------------------------
# Import shapefile and check projection
land_mask <- readOGR("filepath/boundary_minusland.shp")
plot(land_mask)
plot(boundary_proj, add = TRUE)
proj4string(land_mask) # at this point the plotting looks good

# Crop and mask list to the shapefile
raster_list <- lapply(X=raster_list, FUN=crop, y=land_mask)

# mask pixels outside of DS management area
for(i in 1:length(raster_list)){
raster_list[[i]] <- mask(raster_list[[i]], land_mask)
}

# look at resulting raster
plot(raster_list$'19930209')
plot(boundary_proj, add = TRUE)
plot(max_ice_extent, add = TRUE)
#This final plot looks fine until you start playing around with resizing the plot window,
#then the alignment changes so that the land shapefile (land_mask) does not align with the raster pixels


texts: List all tables with ogr_fdw ODBC connection

texts:
I have been using ogr_fdw to connect PostgreSQL/PostGIS databases to FireBird databases via ODBC connection on Windows machines. The connections have been working great, and are immensely useful for transferring data between the databases  - I am using the foreign data wrapper to keep data synced as we transition from FireBird to PostgreSQL.

FireBird is not spatially enabled. On one of the databases, I had created my own "GEOMETRY_COLUMNS" table in FireBird and used SQL to create views and tables with WKT and WKB geometries. This was partly an experiment, but at the same time shoe-horned some spatial functionality into the FireBird database that we needed.

On all of the other databases, I can use "import foreign schema" to connect to all tables, including non-spatial. However, on this one database I am limited to only the tables in the geometry_columns table. Based on this post Can't connect to foreign SQL Server table from PosgreSQL using OGR_FDW, it seems that I need to set variables to list_all_tables and use_geometry columns. However, the post was for MSSQL and I am using ODBC with FireBird. I suspect these settings tie into GDAL: https://gdal.org/drivers/vector/pg.html. I've been trying to create a foreign server that will ignore the GEOMETRY_COLUMNS table and show all tables, not just spatial ones, but have not been successful.

This is my typical create server SQL:

CREATE SERVER fb_db_svr FOREIGN DATA WRAPPER ogr_fdw OPTIONS (datasource 'ODBC:FB_DB')

CREATE SCHEMA IF NOT EXISTS firebird;

I then import schema for individual tables as follows:

IMPORT FOREIGN SCHEMA ogr_all LIMIT TO ("TBL_NAME") FROM SERVER fb_db_svr INTO firebird;

I tried taking a guess at a different foreign server SQL statement based on this https://pgxn.org/dist/ogr_fdw/

CREATE SERVER fb_db_svr FOREIGN DATA WRAPPER ogr_fdw
OPTIONS (datasource 'ODBC:ODBC:FB_DB', config_options 'LIST_ALL_TABLES=YES USE_GEOMETRY_COLUMNS=NO');

I also tried setting a few environment variables on the server:

Still no luck. Is this possible, or is my only option to delete the GEOMETRY_COLUMNS table in the FireBird database?


texts: Is there any way to debug Geoserver WPS deployment?

texts:
There are many reasons for Geoserver WPS processes to fail when initializing, but the most possible is missing dependences. Anyway, it is not always clear which one, since Geoserver provides absolutely no output when a process fails to init for some reason. Setting log4j.category.org.geoserver=TRACE and log4j.rootLogger=TRACE, stdout in resources\log4j.properties does nothing in my Geoserver 2.11.0 standalone.

Only thing I need is a message that shows what is wrong with a process class. Any directions please?


texts:
I am trying to pan to a selected feature. When I run my code it does not show any errors, but when I open the mxd it does not pan to the extend of the selected feature. Here is my current code:

import arcpy, os
ws = arcpy.env.workspace = r"C:\Users\vlaw\Desktop\test\1313\figures"

for mxd in arcpy.ListFiles("*.mxd"):
    current_mxd = arcpy.mapping.MapDocument(os.path.join(ws, mxd))
    for lyr in arcpy.mapping.ListLayers(current_mxd):
        if "site" in lyr.name:

            lyr.replaceDataSource(r'C:\Users\vlaw\Desktop\test\1313\project_data\1313.gdb','FILEGDB_WORKSPACE','Site_Feature_Polygon')

    for df in arcpy.mapping.ListDataFrames(current_mxd):

        for lyr in arcpy.mapping.ListLayers(current_mxd,"",df):
            if lyr.name == "site_boundary":

               query = "'[Type]' = 'Site Boundary'"
               extent = arcpy.SelectLayerByAttribute_management(lyr,"NEW_SELECTION", query)

               df.panToExtent(lyr.getSelectedExtent(extent))
               print lyr.name
    current_mxd.save()


texts: CreateFeatureclass_management Not Displaying in Map - Null Geometry

texts:
I have two feature classes that share a common ID field. I wrote a script to create a new feature class that contains the records in the first FC that are not in the second FC.

It works, in that it creates a feature class with the records that are different, however the new feature class will not display. I have tried running Check Geometry on it, and the "Problem" says "null geometry" for every record in the new FC. It has a spatial reference, but the extent is blank.

I can get the newFC to display by joining it to the first FC and exporting it, but that seems like a workaround that I could avoid.

Anyone have any ideas?

import arcpy

def createMatchingFieldList(fc1, fc2):
    lst1 = [fld.name for fld in arcpy.ListFields(fc1)]
    lst2 = [fld.name for fld in arcpy.ListFields(fc2)]
    return list(set(lst1) & set(lst2))

def createWhereClause(fc, fld_name, value):
    if len(arcpy.ListFields(fc, fld_name)) == 1:
        fld = arcpy.ListFields(fc, fld_name)[0]
        if fld.type == "String":
            where = "{0} = '{1}'".format(arcpy.AddFieldDelimiters(fc, fld_name), value)
        else:
            where = "{0} = {1}".format(arcpy.AddFieldDelimiters(fc, fld_name), value)
    return where

def getPrimaryFieldValues(fc, field):
    return [r[0] for r in arcpy.da.SearchCursor(fc, [field])]

def getSelectCursor(fc, flds, whereClause):
    return arcpy.da.SearchCursor(fc, flds, whereClause)

def diff(a, b):
    return list(set(a) - set(b))

#Get Parameters
#feature class 1
source1 = arcpy.GetParameterAsText(0)
#feature class 2
source2 = arcpy.GetParameterAsText(1)
#new feature class location
my_path = arcpy.GetParameterAsText(2)
#new feature class
new_fc = arcpy.GetParameterAsText(3)
#Name of Key Field in feature class 1
fieldName1 = arcpy.GetParameterAsText(4)
#Name of Key Field in feature class 2
fieldName2 = arcpy.GetParameterAsText(5)
#Spatial Reference
spatial_reference = arcpy.Describe(source1).spatialReference

#my attempt to get the new_fc to display
#define coordinate system
#arcpy.env.outputCoordinateSystem = arcpy.SpatialReference("NAD 1983 UTM Zone 15N")

# create empty feature class
arcpy.CreateFeatureclass_management(my_path, new_fc,"", source1,"","", spatial_reference)

# create a list of field names which are in both featureclasses
flds = createMatchingFieldList(source1, source2)

source1Values = getPrimaryFieldValues(source1, fieldName1)
source2Values = getPrimaryFieldValues(source2, fieldName2)

additions = diff(source1Values, source2Values)

with arcpy.da.InsertCursor(new_fc, flds) as insertCursor:
    for a in additions:
        where = createWhereClause(source1, fieldName1, a)
        insertRows = getSelectCursor(source1, flds, where)
        for r in insertRows:
            insertCursor.insertRow(r)


texts: Finding the shortest path between points to LineString PostGIS

texts:
I have a MultiPolygon. Something like this:

All one colored circles make up one multipolygon. I would like to find the shortest line to connect all one colored polygons.

The first idea was to find all possible lines between one group of points. Then sort these lines by distance and point_id and loop over all points and keep the lowest distance line for each point and then use linemerge to make it one big line.

Expected result would be something like this:

Don't know if this idea would work and don't have any clue how to implement it in postgis. Any ideas?


texts: Evaluate if at least one column contains the value "Yes"

texts: In QGIS, I'm trying to input either 'Yes' or 'No' for a new field based upon the input of values of another 9 fields. I've attached a screenshot of what I'm trying to accomplish. Essentially if there is a 'Yes' in any of the 9 previous fields, then the new field will input 'Yes'. If however all 9 fields are NULL, then the new field will input 'No'.


texts: Intersecting categories within shapefile using ArcGIS for Desktop?

texts:
I have a shapefile composed of fire scar polygons identifiable as unique categories based on year of fire.

I am interested in measuring the area of overlapping fire scars from different years.

Is it possible to perform an intersect of categories within a shapefiles or will I have to split each category into a separate shapefile first?


texts: PostgreSQL Trigger to Change geometry of Point on Different Polygon

texts:
Issue :
I have table in Schema GIS
1. Asset ( Id, User, Begin Date, Status, Master , Geometry)
2. Room ( Id, User, Begin Date, Status, Master , Geometry)
And I have table in Schema Public

1. Asset ( Id, Code, Name, Room, WKT)
2. Room ( Id, Code, Name, WKT )
I have a problem to create trigger in PostgreSQL to change room when asset relocate to some different position of room (change geometry of point on different polygon).
how to adjust for that issue using this query
=================================================================

CREATE OR REPLACE FUNCTION get_address () RETURNS TRIGGER AS
$BODY$
DECLARE
BEGIN
raise notice 'point trigger starts now: %', now();
IF TG_OP = 'INSERT'
THEN
    IF
       (SELECT COUNT(*)
        FROM
            (SELECT t.id
            FROM point as t,
                polygon as s
            WHERE st_Within(NEW.geom, s.geom)) AS foo) > 0
    THEN
        SELECT polygon.address
    INTO NEW.address
    FROM polygon
    WHERE ST_Intersects(NEW.geom, polygon.geom);
    raise notice 'point trigger ends now: %', now();
        RETURN NEW;
    ELSE
            RAISE notice 'no intersection';
            RAISE notice 'point trigger ends now: %', now();
        RETURN NEW;
    END IF;
ELSIF TG_OP = 'UPDATE'
THEN
IF
(ST_Equals (NEW.geom, OLD.geom)=FALSE)
THEN
    IF
       (SELECT COUNT(*)
        FROM
            (SELECT t.id
             FROM point as t,
                polygon as s
        WHERE st_Within (NEW.geom, s.geom)
            AND (t.id  OLD.id)) AS foo) > 0
    THEN
SELECT polygon.address
INTO NEW.address
FROM polygon
WHERE ST_Intersects(NEW.geom, polygon.geom);
    RAISE notice 'Interesction fund';
    RETURN NEW;

    ELSE
    SELECT NULL
        INTO new.address;
  RETURN NEW;
  raise notice 'point trigger ends now: %', now();
  END IF;
ELSE
Raise Notice 'Update of attribute data';
raise notice 'point trigger ends now: %', now();
Return NEW;
END IF;
END IF;
END;
$BODY$
LANGUAGE plpgsql;

CREATE TRIGGER update_point_address
BEFORE INSERT OR UPDATE ON point
FOR EACH ROW
EXECUTE PROCEDURE get_address ()

query from: @neogis https://gist.github.com/neogis-de/a1d08c38d8b9c5d316c7


texts: what is the difference between '$x_at' and 'x', (Labeling) in qgis?

texts:
In properties > Label > placement > data defined > coordinates:

'x'... when I move the label manually,
but what is '$x_at' for? Do they have a huge difference?


texts:
I have a python script that runs a set of queries in parallel using multiprocessing.pool and map_async.

The script from @Aaron works perfectly.

Now I want to integrate the same functions into a PyWPS service, but I do not succeed. The processing does not start and I get no error in return.

Here is my PyWPS script:

import os, sys
from pywps.Process.Process import WPSProcess
import multiprocessing
from multiprocessing import Pool
from pg import DB
from functools import partial....
db = DB(dbname='****', host='****', port=****, user='****', passwd='****')

class Process(WPSProcess):
     def __init__(self):
         # init process
         WPSProcess.__init__(self,...)
         self.input1 = ...
         self.output1 = ...

     def execute(self):
        tup = db.query("""
                SELECT distinct code_cc from comm_pop limit 200
                       """).getresult()

        code_list = list(zip(*tup)[0])

        db.query("""DROP TABLE IF EXISTS tpop.tpop_sample""")
        db.query("""CREATE TABLE IF NOT EXISTS tpop.tpop_sample(...)""")

        pool = multiprocessing.Pool(processes=7)
        func = partial(self.tpop, self.input1.getValue())
        pool.map_async(func, code_list)
        pool.close()
        pool.join()

        out = db.query("""SELECT * from tpop.tpop_sample""").getresult()
        self.output1.setValue(out)

     def tpop(self, tabcouv, code):
        db.query("""
                INSERT INTO tpop.tpop_sample...do something"""% (tabcouv, code, code))

I have no error message in the pywps log file, I have only an info message:
INFO: Status [processstarted][0.0]: Process pop started

I have already tested without the multiprocessing and I got results.
For information I'am using PyWPS 3.2.5 with Python 2.7


texts: Mosaic to Raster - successful but the layer does not show in layers or documents

texts: I have used the Mosaic to New Raster tool and it works however the raster tif does not show in the layers or my documents, why is this happening and how can I solve this issue?


texts:
I'm writing a bachelor thesis about air pollution from point source and I want, to put it to GIS, but I don't know really how...

I'm working with Gaussian plume model, so I have pollution concentration in distance from the source. What else do I need to make GIS pollution model, if I have concentration in distance from source?

I want to make something like this http://ready.arl.noaa.gov/gaussian-bin/usrdisp2.pl in QGIS

I have installed QGIS 2.4, this is my first experience with GIS.


texts: Why does my MapServer WMS not zoom correctly in OpenLayers 3?

texts:
Below code is showing the addition of mapserver file in openlayers 3, it is successfully showing me image in the browser but when i zoom that image it does not zoom. It comes back to originial view.

  var perimetro = new ol.layer.Image({
    source: new ol.source.ImageWMS({
    url: 'http://localhost:8081/cgi-bin/mapserv.exe?',
    params: {
      LAYERS  : 'landuse10' ,
      CRS     : 'EPSG:3857',
      map : 'C:/ms4w/apps/osm/basemaps/osm-google.map',
      mode : 'map',
      FORMAT : 'png'},
    serverType  : 'mapserver',
    extent: [8046840.710000,3906011.335000,8217827.220000,4076997.845000]
    })
 });
 var map = new ol.Map({
   target: 'map',
   layers: [perimetro],
   projection : new ol.proj.Projection({
    code: 'EPSG:3857',
    units: 'm' }),
    view: new ol.View({
       center: [8046840.71, 3907521.8],
       zoom: 16 })
 });

This is my map file layer code:

LAYER
  TYPE POLYGON
  STATUS ON
  PROJECTION
    "init=epsg:3857"
  END
  NAME "landuse10"
  GROUP "default"
  CONNECTIONTYPE POSTGIS
  CONNECTION "host=localhost dbname=osm user=postgres password=postgres123 port=5432"
  DATA "geometry from (select geometry,osm_id,type,name as name from osm_new_landusages_gen1 where type in ('forest','pedestrian','cemetery','industrial','commercial','brownfield','residential','school','college','university','military','park','golf_course','hospital','parking','stadium','sports_center','pitch') order by area desc) as foo using unique osm_id using srid=3857"
  PROCESSING "LABEL_NO_CLIP=ON"
  PROCESSING "CLOSE_CONNECTION=DEFER"
  MAXSCALEDENOM 650016
  MINSCALEDENOM 325008
  CLASSITEM "type"
  LABELITEM "name"

I modified the code in this as per told in previous answer but it is still giving me error of "Unexpected identifier".

<title>MapServer Simple Viewer</title>
<link rel="stylesheet" href="http://openlayers.org/en/v3.0.0/css/ol.css">


texts: Finding shortest route between multiple stops using ArcGIS for Android?

texts:
ArcGIS is already equipped with API service which calculates the shortest path between a start and an end point based on length, time etc.
The following is how the service is called to returns in a RouteResult object, which in turn is used to draw a route on the map (I omit that part):

try {
            NAFeaturesAsFeature rfaf = new NAFeaturesAsFeature();
            // Create the stop points (start at our location, go to pressed location)
            ArrayList<StopGraphic> stopGraphics = new ArrayList<StopGraphic>();
            for(Point pt : directionStops){
                stopGraphics.add(new StopGraphic(pt));
            }
            Graphic[] stopGraphicsPrim = new Graphic[stopGraphics.size()];
            stopGraphicsPrim = stopGraphics.toArray(stopGraphicsPrim);
            rfaf.setFeatures(stopGraphicsPrim);
            rfaf.setCompressedRequest(true);
            A.getRouteParameter().setStops(rfaf);
            //A.getRouteParameter().setImpedanceAttributeName("length");
            // Set the routing service output SR to our map service's SR
            A.getRouteParameter().setOutSpatialReference(wm);
            mResults = A.getRouteTask().solve(A.getRouteParameter());
        } catch (Exception e) {
            e.printStackTrace();
            dialogRouting.dismiss();
            mHandler.post(new Runnable() {
                @Override
                public void run() {
                    Toast.makeText(getActivity(),"directions cannot be found ",Toast.LENGTH_SHORT).show();
                    hideSoftKeyboard(getActivity());
                }
            });
        }

I now have multiple destination points. Say, someone needs to visit 5 coordinates from its own location, but, he/she needs to follow a path which ends up in a shortest travel distance, having visited all the 5 destinations.

I could not find a service from ArcGIS server which provides me such a general path.
If I just provide RouteParameter object with multiple stops with method: A.getRouteParameter().setStops(rfaf); the service will return a path in the given coordinate order. i.e. if my first stop coordinate is the farthest to the start coordinate (own current coordinate in my case), the the route first points to there instead of finding the closest one from the all given stops.

Now, are there any ways to find this shortest general path using any service (maybe apart from ArcGIS), or do I have to make multiple service calls between two points to eventually calculate the shortest general path?
i.e. for N multiple stops I have to make Nx(N+1)/2 point to point service calls with the method: mResults = A.getRouteTask().solve(A.getRouteParameter());, as mentioned in the above code snippet.


texts: Adding intermediate contour values automatically in QGIS?

texts:
We have 100's of elevation contour shapefiles. In each file, only the 10th contour line has values (200, 450, 700, etc..). We have to attribute values for all the intermediate contour lines (2nd line to 9th, 11th to 19th, etc..). While searching for a solution, happened to see this (link below). ArcGIS CalcContour script(as mentioned in the link) will only work in ArcGIS 9.3. So I couldn't even test the script. Please check the "Limitations script CalcContour" note in this link

https://translate.google.co.in/translate?hl=en&sl=pt&u=http://www.processamentodigital.com.br/2014/04/27/arcgis-9-3-manual-tecnico-para-georreferenciamento-de-um-arquivo-dgn-do-ibge/&prev=search

The script will work like (from what I can see from the YouTube video - link below) we have to snap a point on a contour line and have to draw a line. Finish the editing and enter the increment values. The intermediate lines between the starting and ending point will get the respective values. I've attached a screenshot of that video

https://www.youtube.com/watch?v=Uzk5V0MXhUg

Is it possible to do the same or any plugin available for this kind of work in QGIS ?

This is not duplicate of How to automatically plot intermediate contours with QGIS, GRASS GIS or SAGA GIS?

Our intermediate contour doesn't have any values. And we don't have to interpolate any new contour. We have to just assign the respective values for the contours which are plotted. We know the values for, eg: 200m & 450m contour which will be 1st and 10th line.


texts: Formula/ calculation for ESRI map coordinate conversion (meter to decimal degrees)

texts:
I'm new to ESRI mapping. I have a set of points in ESRI, with x, y, z coordinates, which I believe is in meters.

I am now writing an Android program, and I do not want to depend on ESRI, thus I would like to get the coordinates and pin them on Google map.

I have tried to use ESRI and export the coordinates to decimal degrees and pin them on Google map. But the position is not very accurate (varies from original position about 500m).

So I was thinking if there any hard coding calculations that I could do to convert meters to decimal degrees? I'm looking for accuracy at around 1m.


texts:
I have a geotiff image that I am clipping to a geojson polygon using rasterio.  When I load input the image in QGIS, it has an alpha band (4 bands total).  rasterio says it only has 3 bands, which leads me to believe that the rasterio sees the alpha band as a dataset mask.  I am using rasterio.mask.mask to clip the image to the polygon.  The output of rasterio.mask.mask when written to a geotiff is a 3 band image without the alpha band.  Is there a way to include the alpha band/dataset mask in the clipped image?

Here is my code:

src = rasterio.open('<input image path>')
output_meta = deepcopy(src.meta)
output_meta["driver"] = "GTiff"

masked_image, masked_image_transform = rasterio.mask.mask(src, [shape(<geojson geometry>)], all_touched=True, crop=True)

output_meta.update({
    "height": masked_image.shape[1],
    "width": masked_image.shape[2],
    "transform": masked_image_transform
})

with rasterio.open('<clipped image path>', 'w', **output_meta) as dest:
    dest.write(masked_image)


texts: Put together Shapefiles with data for individual years into Multiple Year Space-Time Cube

texts:
I have 16 shapefiles (listed above).  Each shapefile contains county-level data on mortality rates for years 1999 to 2015, respectively.  I organized the mortality rate data in Excel.  I joined the Excel data to an empty US_county_polygon data map by joining FIPS (commonly used to identify counties in census data).

I am creating a space time cube in ArcMap. Space time cubes require a SINGLE input file with data for all years. I have MULTIPLE input files.

I cannot match my multiple files to the same US_County_polygon reference shapefile.  Perhaps, because I join them using FIPS.  Each time I join a new single-year shapefile using FIPS to the master US_County_polygon shapefile, it writes over the prior file.

How can I amalgamate the data to create a space time cube?


texts: Include map overview details in a legend?

texts:
I'm creating a map to show BRUV sites in different marine zones.

There are 6 zones (created an overview for each - please see attached image). The zones are individually displayed in smaller map insets below the legend with site locations and BRUV status.

Is it possible to include the overview information for marine zone in the map legend, for example... I would like the legend to display Zone and 1-6 as map items.

The maps were created manually as I was unsure how to utilise atlas (that's another question but I will not go into that)

UPDATE/EDIT!

Ok, today I created what I desired by hand, manually drawing text boxes, choosing colour, measuring symbol space...  This was time consuming. I don't mind taking the time to do it, but, it's a pain in the bum when I need to repeat the process for another project + I would like to  experiment with colours, text or placement and so forth. It would be splendid if there was an automatic way.


texts:
Under "Settings -> Options -> System" the SVG-paths are set. There are two default svg-paths:

C:/PROGRA-1/QGIS3/apps/qgis-Itr/./svg/

and

C:/Users/jan/AppData/Roaming/QGIS/QGIS3\profiles\default/svg/

I would like to permanently remove both paths.

I can remove them by clicking on the "-" symbol but when I restart QGIS, both of them reappear. Is there a way to remove both paths permanently?


texts: Using Lizmap plugin with QGIS Server to manage user roles

texts: I want to manage user roles in server side but I found out that QGIS Server does not have this option so Lizmap may be a solution. Can I use Lizmap plugin with QGIS but not to use the Lizmap web client because I want a personalized interface?