texts: Finding number of points within each respective polygon in table

texts:
I have:

a shapefile of circles at 15 meters radius and
point data with GPS coordinates.

How can I count the number of points within each circle and have a raw number that can show such? I've tried spatial joins, but I can't seem to find the proper way to determine the number of points.


texts: Clip RGB composite raster by shapefile: the process takes a very long time

texts: RGB-composite raster (6,27 MB, derived from Landsat-8 band images) is clipped by a shapefile (Toolbox - Data management tools - Raster - Raster processing). It takes too long time (foreground processing shows 4% and it goes on at least 1 hour). ArcGIS 10.5.1, RAM 4 GB, Windows 7 x64. How to solve this problem?


texts: How to convert Shapefile to TopoJSON with grouping neighbor shapes by some criteria

texts:
I have a shapefile with countries, which I convert to GeoJSON using ogr2ogr and then converting that json to another TopoJSON, using topojson command-line tool.

Everything works like charm, but now I want to learn to concatenate/merge/join/cluster (just don't know what is the most appropriate term in GIS context) different neighbor shapes to a bigger ones - say, all African countries to one big Africa shape. Let's assume that all Hdata needed (like continent name) is present in the original shapefile.

So, what is the simplest way to merge neighbour shapes while converting shapefile to topojson format?

Sorry if this question is about something trivial, it's just that this is a new topic to me, though I'm very glad I should deal with such issues.


texts: Converting (geodetic pose + cartesian offset) to geodetic position?

texts:
I am in search of an algorithm to solve the first/direct geodetic problem precisely in 3D. The input is:

a geodetic (not geocentric) position (WGS84, height is above the WGS84 reference ellipsoid)
heading (relative to true north)
pitch and roll (relative to gravity, which isn't necessarily pointing towards the center of the earth)
a 3d-offset from the pose P (defined by 1,2,3) to a point A, given in cartesian coordinates (meters).

The output should be the geodetic position (again, WGS84) of point A.

I have read Algorithm for offsetting a latitude/longitude by some amount of meters and this is how I've computed the solution in the past. It is obviously computationally efficient, but not precise enough anymore.

Currently, my idea is to

convert 1. to ECEF, populating the 4th column of a 4x4 matrix with the solution.
incorporate the angles/orientation into that matrix
multiply the matrix with the 3d vector to A to compute the position of A in ECEF.
convert this back to Geodetic WGS84 (which isn't quite straightforward).

Of course, I still would like things to be fast, but precision has a definite priority here.

My questions:

What do you think of this approach? Does it work, is it precise? (I fear that compensating for pitch/roll not being relative to the center of the earth would probably imply a LOT of additional complexity).
Is PROJ.4's pj_transform() capable of computing the conversions in steps 1 and 4? If not, what C/C++ library can you recommend? This must be a common routine in some library?!

Note to @whuber who referred to 30448! I think I'm slowly starting to resolve misunderstandings, getting an idea of what you do. Questions about your answer to that question:

When you say "geocentric", you mean cartesian coordinates from earth center (=ECEF), right? (I first interpreted "geocentric coordinates" just like "geodetic coordinates", except the latitude being determined using earth-center instead of ellipsoid-normal. Now I know that "geocentric latitude" != "geocentic coordinates". Whew.)
In the first image of 30448, you show "the geocentric frame at the North Pole". The origin of that coordinate system really is at the center of the earth, right? But if your first image started at earth's center, simply rotating wouldn't get you to the origin of the local coordinate system. So why does it start at (0,0,6370) in the first picture? You could have also chosen (6370,0,0), but then rotations would be different. Is it that each initial offset corresponds to an Euler convention (Z-Y'-Z'')?
If so, it then seems to me that my drafted approach is very much similar to your solution at 30448. Are there other approaches that are more precise? Given I compute using doubles, what accuracy do you expect?
Slightly OT: At the end, you add a local Z displacement motivated by differences against MSL - but I thought the WGS84 datum you referenced earlier is not a tidal/MSL datum, but an ellipsoidal. Assuming you sticked to WGS84, you'd add the local Z offset based on the distance to the ellipsoid, right?

@whuber, @martin f provided quick replies and patience - I read the whole day, trying to understand more.


texts: Deleting an invalid inmemory layer from the Table of Contents

texts:
I am writing an ESRI addin (in c#) where in_memory layers are used. When the user opens the map document and activates the tool, the tool checks whether the in_memory layer is valid or not, and if it is not valid (has a red exclamation mark) it should delete the layer and recreate it for the current session.

I am using the following to delete the layer

mxDoc.FocusMap.DeleteLayer(dsAOI);

but once it deletes it i then get a message box saying:

  Object not set to an instance of an object

I have tried to disconnect the layer before I delete it using IDataLayer2.Disconnect

   IFeatureLayer toDelete = dsAOI as IFeatureLayer;
   IDataLayer2 dl2 = toDelete as IDataLayer2;
   dl2.Disconnect();

but that results in the error message:

  Error HRESULT E_FAIL has been returned from a call to a COM
  component

Has anyone got ideas on what I might be doing wrong, or a better way to get round this.


texts: ArcObjects .NET 9.3 upgrade to 10 breaks IRasterDataset.Copy(...) code

texts:
I have to upgrade an application from ArcObjects 9.3 (C# .NET) to 10. Unfortunately, there seems to be some incompatibility regarding the IRasterDataset.Copy method. The following code worked fine in 9.3:

my_raster_dataset.Copy("my_copied_raster", my_sde_workspace);

my_raster_dataset is a raster dataset called "INPUT_RASTER" loaded from a file geodatabase. my_sde_workspace is an ArcSDE 9.3 geodatabase running on SQL Server 2005.

Using ArcObjects 10, I get the following error message:
System.Runtime.InteropServices.COMException: The table was not found. [VAT_INPUT_RASTER]
   at ESRI.ArcGIS.Geodatabase.IRasterDataset.Copy(String copyName, IWorkspace copyWorkspace)

Any ideas or suggestions? Do I have to upgrade ArcSDE to 10 to make this work?

Thanks in advance.

EDIT - Found the answer

Found the problem. I was checking out esriLicenseProductCodeEngine instead of esriLicenseProductCodeEngineGeoDB. Strange, that the code worked with ArcObjects 9.3 since 2010 and suddenly stopped working with ArcObjects 10.


texts: Floating point: understanding their inaccuracy

texts:
I am trying to understand why a floating point is inaccurate. Please bear in mind that I am not a mathematician.

What is exactly the difference between floating point and double integers? and why use one rather than the other?

Thanks for any "easy to understand" replies,

Robert


texts: How can I recalculate an attribute based on overlap with another layer?

texts:
I have 2 layers: The first one is a point Layer (converted from SRTM) with height information. The second one is a polygon layer (shows the forest).

I want to reduce the height value - let's say about 15m if the point value is in a forest area. How can I do that?


texts: Remove tables - not just layers - from the map document

texts:
How do you remove tables/table views from the table of contents (TOC) in the map document?

I've written code that cleans things up after each tool run, but I cannot figure out how to get rid of the tables. The most I've been able to do is delete them from the geodatabase, but then they just remain in the TOC with broken connections:

# ---------------------------------------------------------------
# Clean up the gdb & map docoument from previous script executions
# ---------------------------------------------------------------
goodFCs = ["NYC_Subway_Stations", "NYC_2010_CensusBlocks_BasicDemographics", "NYC_Borough_Boundaries"]

mxd = arcpy.mapping.MapDocument("CURRENT")
for df in arcpy.mapping.ListDataFrames(mxd):
    for lyr in arcpy.mapping.ListLayers(mxd, "", df):
        if lyr.name not in goodFCs:
            arcpy.mapping.RemoveLayer(df, lyr)
    for tbl in arcpy.mapping.ListTableViews(mxd):
        arcpy.Delete_management(tbl)
mxd.save()

allFcs = arcpy.ListFeatureClasses()
for fc in allFcs:
    if fc not in goodFCs:
        arcpy.Delete_management(fc)
    else:
        pass

Image 1: Once it's run a second (or third or fifth, etc.) time, my script will successfully remove the map layers highlighted in yellow, but not the tables inside the red box:

Image 2: After the second script run, yellow map layers will have been successfully replaced with new yellow map layers, but the old tables will remain with broken links. I want to get rid of them completely:

I'm aware that I don't have any code written to remove the tables, but that's because I don't know what code to write. Tables aren't inside of the dataframe so I can't use arcpy.mapping.RemoveLayer(), which takes a dataframe as its first required parameter...


texts:
I know this question has been addressed in the past but I am having trouble implementing what has been suggested here on my set up and I was hoping someone could help me hack through all of the information.

I am running GeoServer on port 8080 (version 2.5) on Windows server 2012. I would like to have an OpenLayers 3 application on port 80. I can draw the layers without issue but getting information using this example just isn't working for me.

I am new to all of this so configuring anything server side is a challenge for me to say the least. that being said I have tried the following things based on what I have read to no avail.

Presenting GeoServer on port 80 on a Windows IIS Server

Application Request Routing module and the URL rewrite module on IIS

I have not attempted an open layers proxy since all of those posts seem to come from a few years ago.

I apologize if this is a repeat but to my untrained mind I see several possible solutions and none of them are firing for me.


texts: Open mail client when clicking feature in Attribute Table

texts:
One of the attributes/features of my layer contains an email address. I'd like to be able to open the default mail client by clicking on the mail address from the Attribute Table.

After some initial investigation I've found, what seems to be, two different possibilities.

1) I could make an action 'Layer Properties/Actions', or

2) perhaps I could change the fields 'Edit widget' in Layer 'Properties/Fields' to a 'Web view', and that would take care of it(?)

On option 1) I'm a bit unsure how to create an action programmatically using PyQGIS.
I've found that a QgsVectorLayer has a function actions() that returns a QgsAttributeAction object. The QgsAttributeAction has an addAction function, so I should be able to use it like this...

layer = iface.activeLayer()
aa = layer.actions()
aa.addAction(...)

but I'm not a 100% on the arguments. From the documentation:

addAction(QgsAction::ActionType type, QString name, QString action, bool capture=false)

Anyone done this before that could give a few hints?


texts:
If I set up a Web Feature Service that provides point (location) data, can my WFS also provide a bunch of additional data for events that occurred at those points?  Is there a better way to do that?  Would a GIS tool such as an ESRI product even know how to deal with the event data?  Should the WFS only include the point data?

What am I trying to do?  I have a database full of data that has samples/observations that occurred at monitoring locations over time.  Our database users want to be able to pull the data into their GIS tools over the web (not through database connectivity) and then do all sorts of cool stuff with the data in their GIS tools.  Is WFS the rights type of service for what I'm trying to do?


texts: How to suppress noises of this edge-detected image?

texts:
I've recently developed a new approach to detecting image edges. Everything seems to be normal except the noises! As you can see in the following image, there are many noises I'm not able to handle them properly. I used some bluring filters on the output, but the results were not as acceptable as those I expected. Moreover, since my digital number results generally range from 0.2e-8 to 0.5, thresholding methods, especially those of Python which only work on 8-bit images, have not been effective at all.

I would like to know your suggestions and comments on this denoising problem.


texts: How to join an external table to a set of polygons, and symbolize on the joined values, in a web map?

texts:
Say I have a table showing a number of statistics - for example, various gun crime by US State measures. I also have a dataset of state boundaries.
How can I join the table to the features in a web map, and symbolize the polygons using attributes from the table? The joined table must be specified within the web application (rather than hard-coded via ArcMap) since the statistics will be derived on-the-fly and will be specific to each user.
I need the ability to change the symbology on-the-fly, for example by choosing a different crime statistic to map, in which case the map should update immediately.
ArcGIS Server
Using ArcGIS Server 10.1 the REST API supports table joins. Are there any sample scripts which step through how to set this up?
Other environments?
Can the same thing be achieved using FOSS software?
Thanks. (This is a modified version of my earlier question since that hasn't yet received any responses)


texts: How to save attribute query on qgis 2.14? I have more than 20 queries. The recent queries save only 20 queries. If I can save them, it will be more effiecient.


texts: How do I extract building elevation values from Multipatch GDB

texts:
I have a GDB Multipatch Buildings layer that I open with ArcScene or ArcGlobe.
I want to extract building elevation values, but this is not stored in the attribute table. I can run a "Add Z Information", but I was wondering if someone can guide me on how to extract it with FME.


texts: How to control state of buttons in Python Add-In?

texts:
I have bunch of buttons to perform separate tasks within a Python Add-In. Buttons are greyed out and required to be enabled after initial tools run. In the following example, tool1 is turned off after execution using its ID and button1 is enabled. Classes are arranged alphabetically at the time of Add-In construction. This configuration of classes works sometimes and fails at others with message: "NameError: global name 'button1' is not defined". This is confusing as to why it runs sometimes and fails at others. Can anyone please help with explanation as to why is it the case and how to better control the state of buttons without error?

class Calc(object):
    """Implementation for test_addin.button1 (Button)"""
    def __init__(self):
        self.enabled = False
        self.checked = False
    def onClick(self):
        """do something"""

class process(object):
    """Implementation for test_addin.tool1 (Tool)"""
    def __init__(self):
        self.enabled = True
        self.cursor = 3

    def onMouseDownMap(self, x, y, button, shift):
        """do something"""

        tool1.enabled = False
        button1.enabled = True


texts:
ArcGIS Desktop 10.6 won't open. It just shows the splash, but then doesn't do anything. No error, no nothing. I turned my laptop off and back on; I renamed the normal.mxt file; I deleted said file; I deleted the folder it's in. None of these actions worked.

What else can I try?


texts: Is anyone aware of an open source exhaustive list of geocoded (WGS84 lat long) zip codes for the Republic of Ireland.


texts: No Data beyond 1:2290 with Openlayers Plugin

texts:
No Data beyond scale 1:2,290 with Openlayers Plugin. I have not installed the plugin that you should deinstall to fix the problem: Unable to zoom beyond scale of 1:2,256 in QGIS?

How can I fix this?


texts: Some problems with understanding zoom, resolution on map and device

texts:
I have some questions over map resolution and zoom over Mapbox GL in maps and device
I have found a similar reference to understand
https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Resolution_and_Scale
Within it there is an equation
scale = 1 : (screen_dpi * 1/0.0254 in/m * resolution)
My queries are

IF using the above equation, the result of scale would be   inches/pixel?
ie from the page at zoom level 18, I would be getting 1/2116 inches / pixel for 90 dpi
Am I right?
As I am using Mapbox (512x512), it would mean going down 1 level ?
I am using Ipad as display that would mean 264 dpi
REference: https://www.apple.com/sg/ipad-pro/specs/
Lastly I am using Mapbox studio to style a line(vector) from my database, I am trying to determine in my style how thick the line (in pixel) should be to be as close as possible to real life markings. Will the style I made in PC be the same as in IOS (as in same real width) ? I am displaying an offline map from Mapbox in Ipad

Thanks


texts: projectAs on original geometry vs part geometry

texts:
I need to get XY and lengths of a part of a multipart feature in different CRS.  I can use the projectAs just fine when the it's the original feature geometry, but when I try to use it on the feature part, it doesn't do anything.

with arcpy.da.SearchCursor(fc,['NAME','SHAPE@']) as sc:
for row in sc:
    geometry = row[1]
    print "GeoType {0}".format(geometry.type)
    print "Orig Geo: {0}".format(geometry.firstPoint)
    print "Orig GeoPrj: {0}".format(geometry.projectAs(spatialRefSPNM27).firstPoint)
    geometryPart = geometry.getPart(0)
    print "Part Type {0}".format(arcpy.Polyline(geometryPart).type)
    print "Part: {0}".format(arcpy.Polyline(geometryPart).firstPoint)
    print "PartPRJ: {0}".format(arcpy.Polyline(geometryPart).projectAs(spatialRefSPNM27).firstPoint)

When I run that test I get this results.  Notice how projectAs works on the 3rd line but not on the 6th.

GeoType polyline
Orig Geo: -103.950469357 32.198115364 NaN NaN
Orig GeoPrj: 618430.281549183 436005.055399118 NaN NaN
Part Type polyline
Part: -103.950469357 32.198115364 NaN NaN
PartPRJ: -103.950469357 32.198115364 NaN NaN


texts: What system tables get INSERTS when registering object with Oracle EGDB (ST_Geometry)?

texts:
In ArcGIS Desktop, I can register ST_Geometry objects with the geodatabase (Enterprise 10.7.1; Oracle 18c):
Register with Geodatabase

Register feature classes, tables, views, and raster layers that were
created in the database using ArcGIS clients with the geodatabase.
Once registered, information about the items—such as table and column
names, spatial extent, and geometry type—is stored in the
geodatabase's system tables, allowing these registered items to
participate in geodatabase functionality.

It would be helpful to know what happens when I register an object with the geodatabase — for troubleshooting purposes.

For example: Occasionally, GDB objects become orphaned; the database object has been deleted, but the geodatabase system records still persist (an artifact of the object shows up in ArcCatalog that can't be deleted via ArcCatalog).

It's my understanding that the action of registering an object with the geodatabase INSERTS records into system tables.
Question:
What system tables get the INSERTS?

TABLE_REGISTRY
ST_GEOMETRY_COLUMNS
Other?

Note: my feature classes have domains, but not relationship classes.


texts: Displaying markers instantly when loading in cluster?

texts:
I am loading a big json for about 10000 pins with some text, about 3-4mb.
And it takes some time to load all the pins and display them.

Is there a way to diplay the pin right after loading?

So to see the progress of pins loading?

Here is my jquery json loading

$.ajax({
         type: 'POST',
         url: 'api.php',
         data: {},
         success: function(data) {

             data.map(function(item) {

                 var marker = L.marker(new L.LatLng(item.lat, item.lng), {
                         icon})
                     .bindPopup(content)
                     .on('click', function(e) {
                         console.log(marker);
                     });
                 clusterMarkers.addLayer(marker).addTo(map);
                             })
                      },
         //contentType: "application/json",
         dataType: 'json'
     });

Update: the answer and example is here https://github.com/Leaflet/Leaflet.markercluster#chunked-addlayers-options


texts: Python Toolbox parameters set up error: Value is Required

texts:
My Python Toolbox needs to complete two different kinds of functions, Urban Growth Model and Energy Sector Model, for the Model Type parameter. When I choose Urban Growth Model, the tool window will enable the Plat Features in GDB, and disable the last four parameters. When I choose Energy Sector Model, the tool window will enable the last four, but not the Plat Features in GDB. All of them are the required parameters for the corresponding model type.

The question is I can't let the Plat Features in GDB blank when I run the Energy Sector Model and I can't let the last four parameters blank when I run the Urban Growth Model. I want to find a way to leave the parameters empty when I choose another model type to run.

Any suggestions? I am using ArcGIS Desktop 10.3.1 Advance License.


texts: Query for Selecting only Selected Fc's from a GDB

texts:
I need to version n number of fc's which start with sde."ENT-QA". I need to select only those feature classes and run versioning for it. Can anyone help me in writing the query for selecting only fc's which start with sde."ENT-QA". Remeber ENT-QA is mentioned has quotes. Thanks in advance.

import arcpy
arcpy.env.workspace= r"xyz.sde"
fdList = arcpy.ListFeatureClasses("","","")

for fc in fdList:
    print fc
    arcpy.RegisterAsVersioned_management(fc, "NO_EDITS_TO_BASE")
    print "fc registered"


texts: Identifying any gaps "broken span" in network of ducts using ArcGIS for Desktop?

texts:
I have a huge network of ducts and I need to identify if there are some gaps in the network, for example if some lines are not connected to each other in ArcGIS so I would like to select them.

I was using Selection by Location  by target the same layer using the method "touch the boundary of the source layer feature" and the result is not accurate.


texts: Issue with ModifyFeatures using FeatureStore

texts:
We are trying to modify features using SimpleFeatureStore and Transaction as shown below:-

Transaction transaction = new DefaultTransaction("Example1");
SimpleFeatureStore store = (SimpleFeatureStore) sfSource;
store.setTransaction( transaction );

Filter flt = ff.within(ff.property("the_geom"),ff.literal(polygon));

SimpleFeatureType featureType = store.getSchema();
try {
   store.modifyFeatures( "0", object, filter );
   transaction.commit();
}
catch( Exception eek){
   transaction.rollback();
}

and This results us in an error.

ERROR 6: SetFeature() must be used on a feature with a FID.
java.lang.RuntimeException: OGR Error: General Error at
org.gdal.ogr.ogrJNI.Layer_SetFeature(Native Method) at
org.gdal.ogr.Layer.SetFeature(Layer.java:145) at
org.geotools.data.ogr.jni.JniOGR.LayerSetFeature(JniOGR.java:287) at
org.geotools.data.ogr.OGRDirectFeatureWriter.write(OGRDirect
FeatureWriter.java:137) at
org.geotools.data.InProcessLockingManager$1.write(InProcessL
ockingManager.java:337)

This process is done on OGRDataSource.

After this we are trying to modify values using featureIDs but that results me in same error.

Transaction transaction = new DefaultTransaction("Transaction");
    Maplayer MpLayer = this._MapListLayers.get(cmbLayers.getSelectedItem().toString());
    FilterFactory2 ff = CommonFactoryFinder.getFilterFactory2();
    Filter flt = ff.id(ff.featureId("long cl0805143437.44177"), ff.featureId("long cl0805143437.44178"),
            ff.featureId("long cl0805143437.44179"), ff.featureId("long cl0805143437.44180"),
            ff.featureId("long cl0805143437.44181"), ff.featureId("long cl0805143437.44182"),
            ff.featureId("long cl0805143437.44183"), ff.featureId("long cl0805143437.44184"),
            ff.featureId("long cl0805143437.44185"), ff.featureId("long cl0805143437.44186"),
            ff.featureId("long cl0805143437.44187"), ff.featureId("long cl0805143437.44188"),
            ff.featureId("long cl0805143437.44189"), ff.featureId("long cl0805143437.44190"));
    try {
        MpLayer._SimpleStore.setTransaction(transaction);
        MpLayer._SimpleStore.modifyFeatures("Equipment", new String("Equipment"), flt);
        transaction.commit();
    } catch (Exception e) {
        e.printStackTrace();
        transaction.rollback();
    }

Please help us to modify OgrDatasource.
We have done both the operations on CSVDataSource it works perfectly.
Issue is only in OGRDataStore.


texts:
How could I extract all peaks from a DEM (large mountainous area) using GRASS?

Peak means something like: a pixel which is sourrounded by other pixels with smaller elevation numbers in a area arround the "peak-pixel" of about 100m.


texts: Obtaining elevation data for each grid cell in raster image using QGIS

texts:
I am obtaining a raster image from the internet for elevation data. There are two-three files available for a single location :

National Elevation Dataset 3 Meter, 1 map,
National Elevation Dataset 10 Meter, 1 map,
National Elevation Dataset 30 meter, 1 map.

What I want to do is for each raster image obtained, I want to make a grid of specific cell separation (preferably in terms of km/miles on actual raster image) and for each cell's (in grid) centre point obtain the latitude longitude information and associate it with mean elevation level of that grid cell. Ultimately I want a tabular format for coordinates along with their mean elevation value.

I am working on QGIS.


texts: Oracle SDO_GEOM.SDO_MBR Wrong Bounding Box Coordinates

texts:
I am using Oracle Spatial to match geo-specific images to the appropriate country they represent.

I have a table with a row for each country. Each row consists of a SDO_GEOMETRY column for 3 items: SHAPE (overall boundaries), CENTROID (centerpoint of the country) and MBR (Minimum Bounding Rectangle). The Centroid and MBR fields are auto-computed using SDO_GEOM.CENTROID and SDO_GOEM.SDO_MBR respectively.

Mostly there are no issues, however... For some countries that span the 180deg Longitude line, the MBR created is wrong. I believe this is caused by Oracle handling the (co)Ordinates as X,Y values.

For instance, using Russia as an example. The LowerLeft value of the MBR is returned as (41.19268056, -180) while the UpperRight value is returned as (81.858710028, 180). When mapping this bounding box, it covers a swatch that spans 360deg (and encompasses many countries along the way).

Another example is Kiribati. The LL returned is (-11.46111418, -174.543405728) and the UR is (4.7230899, 176.850922). Mapping this shows the bounding box as NOT crossing the 180Lat line, but rather spanning -174 to +176, when it should be progressing from 176...178...180...-178...-176...-174.

Any thoughts or ideas as to an solution?

===========================

EDIT

SQL used to populate the MBR row:

UPDATE schema.table SET mbr=sdo_geom.sdo_mbr(geometry);


texts: SQL / PostGIS: Expanding polygons contained inside another polygon until one vertice touches bounding box

texts:
We have a problem that we are not quite sure how to approach, we have a number of polygons contained inside a larger bounding polygon, We need to expand these internal polygons as much as possible but only until a single edge touches the containing polygon.

I hope a diagram will help with my explanation:

The polygons in red can be assumed as already created, the dotted polygons in green indicate the final polygon that we want to create (enlarged just enough until one side touches the walls of the bounding polygon).

This will be run on big data (hundreds of millions of polygons with hundreds of millions of inner polygons) as such efficiency is one problem we are facing. In code this is reasonably trivial, however its not feasible for the dataset we have.

We have a large cluster setup that can process Geospatial SQL efficiently which we intend to use for this task but other than arbitrarily passing a value to an expand function we aren't quite sure how to do this intelligently over massive amounts of data.

I believe we need to scale the polygon using ST_Expand and stop when ST_Intersects on the main polygon is true - I just have no idea how to do that in a single SQL call.


texts:
I am trying to randomly select some points from a shape file using ogr2ogr and an SQL query via sqlite but it seems to baulk at the sqlite version:

ogr2ogr -sql "SELECT * FROM probe_2004to2012 WHERE D2004 > 0 ORDER BY RANDOM() LIMIT 10" -dialect sqlite ./Probing/probe_2004_r10.shp ./Probing/probe_2004to2012.shp

ERROR 6: The SQLite version is to old to support the SQLite SQL dialect

Given that I tried updating sqlite3 in my /usr/bin/ from the jan 2013 version and it is now the latest (Oct 2013) version and ogr2ogr still throws up this error I suspect it may be something else that's wrong.
I know I can select random points in QGIS using the "Research Tools" toolbox but I would really like to be able to script this as I have a lot of random selecting to do.

-update:
its seems ogr2ogr supports

-dialect sqlite3

as a valid flag and that seems to get a simplified sql query to parse but the random() selection won't get through.

After the discussion below it seems that ogr2ogr is not parsing the -dialect flag "properly". The sqlite3 flag is ignored and ogr defaults to its built in sql, which is why the simple queries work. Why the sqlite flag is coming up ERROR 6 I don't know. Yet.


texts: Generate a Grid / Points around a single point in QGIS

texts:
I have Baseline points. Around these points I would like to generate a grid that is a at a specified distance. Lets say that this grid will be 10km x 10km. It should have a limited amount of rows and columns for-instance 3 x 3. Please see what I mean in the attached picture.

I can generate a Fishnet and then move the grids by hand but that will be labour intensive due to the amount of Baseline points I have.


texts: QField - Layers and basemap are toggled off by defeault when opening a project

texts:
I just created a data collection app using QGIS and QField. Everything looks fine on QGIS and I exported using QSync plugin. Whenever I open the project on QField, all layers are toggled off. I have to go through every single layer to check the box 'show on map'. This gets very annoying when working with multiple layers.
Does anyone have an idea on how to resolve this issue? I would like to have all the layers visible whenever I open the project.


texts: What steps would I need to take to create an offline 7.5 x 7.5 min USGS Map Viewer similar to this but entirely dependent of an external hard drive for the large data/PDFs? It unfortunately, has to work off an external hard drive because the intent is to give our Search & Rescue folks the ability to print maps on the fly of the desired area when they are deployed across the country for natural disasters. Additionally, it is assumed that connectivity is unlikely for the first several days during a Search & Rescue deployment.


texts: Is there a website that has collection of GeoJSONs?

texts: I specifically look for GeoJSON for all dry lands in the globe.


texts:
I need to develop a web-based map using existing shapefiles.

Using the web-map, we must be able to select a region and produce contours for that selected region.

Can any one suggest the best open source solutions for this?


texts: IDL 7.1.2; How to plot negative correlation in Taylor diagram

texts:
I am working on a model's run performance evaluation, three runs of GFDL-CM2.0 and five runs of GFDL-CM2.1 are used for temperature and precipitation for the years 1979-2000. The objective is to validate the runs with respect to observed data sets e.g CMAP for precipitation and NCEP/NCAR for temperature. In this regard I want to plot a Taylor diagram. I am using "cgtaylordiagram", but this code generates only positive correlations on plot and my data have some negative correlations as well with respect to the observed dataset. Further more I need some guidance regarding the RMS corporation in Taylor diagram.
The following code is being used and shows only positive correlations on plot:

  !P.CHARSIZE = [0.5]
     labels = ['2.0R1', '2.0R2', '2.0R3', '2.1R1', '2.1R2', '2.1R3', '2.1R4', '2.1R5']
     stddev = [0.23, 0.15, 0.39, 0.27, 0.28, 0.12, 0.18, 0.08]
     correlation = [-0.96, -0.86, 1.00, 0.95, 0.33, 1.00, -0.11, 0.47]
     ref_std = 0.22
     stddev_max = 0.5
     OVERPLOT=overplot
     cgTaylorDiagram, stddev, correlation, REF_STDDEV=ref_std, $
       STDDEV_MAX=stddev_max, LABELS=labels, OUTPUT='taylor56.png'
   END

Is there alternative code or method that can plot negative correlation as well in the diagram?
Note: I did not find any solution for e


texts:
I have a maps where i want to show 6 different sites around europe. Instead of creating 6 different maps, i want to use the atlas generator.

I'm using QGIS 2.18.3 (Las Palmas).

I have made a coverage.shp including 6 polygons to show the extent of my different maps. Also, i set visibility presets (always showing 2 Layers: one is my backround google satellite image and one is always my vector shape).
I have consulted a previous question:
(QGIS Atlas generation for multiple map extents?).

The Problem I have is that i want each map to have a different scale. I added an attribute to my coverage layer with the scale of my choice (eg. 30000, 16000...etc.). In the Print composer I went to the scale of the Map, data defined override > Edit, and set the scale of my attribute table. This is however not working.
where am I going wrong?


texts: Fails ST_intersects(geom, select statement) for drive time

texts:
I am trying to select houses in a drive time isochrone.
I have the code working fine to make the isochrone.
I can't place the results in a seperate table
I need the isochrone in a where clause

Is there anyway to achieve this?

Here is where i get to:

SELECT "uk_properties"."id" AS "id",
  "uk_properties"."geom" AS "geom",
  "uk_properties"."address" AS "address",
  "uk_properties"."type" AS "type"
FROM "public"."uk_properties" AS "uk_properties"
where ST_Contains(
      SELECT cast(ST_convexhull(ST_collect(tableauways.the_geom))
        as geometry(Polygon,4326)) from tableauways join (
SELECT * from pgr_drivingDistance(
        'SELECT gid as id
        , source
        , target
        , length_m/(26.82*10) AS cost from tableauways'
        , (select source from tableauways
                    order by tableauways.the_geom <->
                        (
                        select geom from asda_items
                        where storedesc = 'ISLE OF DOGS'
                        limit 1)
            limit 1),
                10,FALSE)
) as route on tableauways.target = route.node

, "uk_properties"."geom")

The first half of the st_contains (the isochrone generating bit) works perfectly on its own and returns a single geom.

Any ideas why I get:
ERROR:  syntax error at or near "SELECT"
LINE 14:   SELECT cast(ST_convexhull(ST_collect(tableauways.the_geom)...


texts: Running ArcGIS Model tool from Python script?

texts:
I've broken a complex task down to a very simple exercise and still can't get it to work... I need to set up a Python Script to Run a Model.

model: "SelectTest"
toolbox: "INTEGRITY_TOOLS.tbx".
location: c:\Work\INTEGRITY\INTEGRITY_TOOLS.tbx

import arcpy
arcpy.ImportToolbox (r"c:\Work\INTEGRITY\INTEGRITY_TOOLS.tbx", "SelectTest")
arcpy.SelectTest()


texts: QGIS start-up error message: "Couldn't load plugin MetaSearch due to an error when calling its classFactory() method"

texts:
I'm an Ubuntu 16.04 (64-bit) user and I just recently upgraded my QGIS instalation (from 2.14, "Essen", to 2.16, "Nødebo").

I'm not a programmer and I know little about the details of how python works, much less about how QGIS and python relate, so when I get this message at start-up, I don't really know what to do about it:

Couldn't load plugin MetaSearch due to an error when calling its classFactory() method

ImportError: No module named requests
Traceback (most recent call last):
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 333, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "/home/jmb/.qgis2/python/plugins/MetaSearch/__init__.py", line 29, in classFactory
    from MetaSearch.plugin import MetaSearchPlugin
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 607, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/jmb/.qgis2/python/plugins/MetaSearch/plugin.py", line 31, in
    from MetaSearch.dialogs.maindialog import MetaSearchDialog
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 607, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/jmb/.qgis2/python/plugins/MetaSearch/dialogs/maindialog.py", line 44, in
    from owslib.csw import CatalogueServiceWeb
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 607, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/owslib/csw.py", line 27, in
    from owslib.util import OrderedDict
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 607, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/owslib/util.py", line 35, in
    import requests
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 607, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: No module named requests

Python version: 2.7.11+ (default, Apr 17 2016, 14:00:29) [GCC 5.3.1 20160413]
QGIS version: 2.16.1-Nødebo Nødebo, 8545b3b

Python Path:
/usr/share/qgis/python
/home/jmb/.qgis2/python
/home/jmb/.qgis2/python/plugins
/usr/share/qgis/python/plugins
/usr/lib/python2.7
/usr/lib/python2.7/plat-x86_64-linux-gnu
/usr/lib/python2.7/lib-tk
/usr/lib/python2.7/lib-old
/usr/lib/python2.7/lib-dynload
/usr/local/lib/python2.7/dist-packages
/usr/lib/python2.7/dist-packages
/usr/lib/python2.7/dist-packages/PILcompat
/usr/lib/python2.7/dist-packages/gtk-2.0
/home/jmb/.qgis2//python

  Note that I can ignore this message and work in QGIS without further issues.

Anyways, can anyone tell me how to fix this?

In case it matters, here are the commands used for upgrading QGIS:

sudo apt-get update
sudo apt install qgis python-qgis qgis-plugin-grass

And these are the lines in my relevant lines from the /etc/apt/sources.list
file:

deb http://qgis.org/debian xenial main
deb-src http://qgis.org/debian xenial main

  Note2: can this be originated by the non ASCII character in the QGIS version "Nødebo"?


texts: Getting "Points of LinearRing do not form a closed linestring"


texts: Changing Virtual field properties in QGis 3.0

texts:
I have some existing virtual field with set Length and Precision

An example:

The formula was first set to rounding it by 2 digits
With the current settings it displays 10.00 instead of 10 because of the precision.

How do I change those two properties in the existing virtual field without deleting it and setting up a new one?

QGis displays numbers with points instead of commas.
For calculating with excel I do need commas.
A solution to this would also help a lot.


texts:
I'm looking for some insights into the accuracy (x,y,z) of oblique imagery from Pictometry in terms of providing accurate measurements for property evaluation. Is the high cost of this product justified by its potential to provide not only high visual/pixel resolution, but also highly accurate measurements? Essentially, I'm curious about alternatives that are more cost effective and provide equal or better accuracy in the 2-D plane.

I find a lot of sales/marketing info from Pictometry, but not much quantitative metadata in terms of consistency for what DTMs they use to rectify the images. There seems to be a lot of momentum behind this company and I want to know if it's justified.


texts:
I develop a map with live data. The data is stored in a PostGIS (PostgreSQL) database. Geoserver is used to provide the WMS.

My problem is that I need to redraw my layer, let's say, every 5 seconds to see changes of the data. The most common way seems to be something like this:

window.setInterval(function(){
    var params = source.getParams();
    params.t = new Date().getMilliseconds();
    source.updateParams(params);
    source.refresh();
    window.setTimeout(function(){
    map.updateSize();},200);
}, 5000);

Unfortunatelly it does not work for me. New data is shown only if I zoom or move the map. I cannot see GET requests from the browser.

How can I update my map?


texts: QGIS Heatmap by value (NOT by point-densty)

texts:
I need a heat map generated by points value irrespective of how close or far apart the points may be in QGIS 2.6.1   So far I'm using the raster/heatmap/heatmap function. I need the same colour to run through all the points of the same value.

From the heatmap dialog, i'm using "use weight from field" which points to the right attribute, but I'm not getting the right result. As you can see this image shows that the heatmap that is generated has the same colour that runs through different values so it can't be right.

How do i get a heat map that is NOT controlled by point-density


texts: Populate a layout text element with the description from a coded value domain

texts:
ArcGIS 10.1+ lets you easily put dynamically populated text boxes in your print layouts.  Doing so puts some text in the text element that looks something like this:

<dyn type="page" property="attribute" field="sampleTable.sampleAttribute"domainlookup="true"/>

It works fine for me when I have it just pulls from the index feature class.  My problem is that I have data in other feature classes in the same database that need to be joined to the index layer so I can pull relevant attributes out of there and display them in other text elements.  As soon as I join one of those feature classes to the index feature class the domain lookup no longer works regardless of what it is set to.  That is to say I only get the coded values in my text elements.  ESRI has logged this bug (#NIM104451), but I need a workaround sooner than when (if) they get around to fixing it.

My first thought was to use some python to take care of this, but it's getting more unwieldy the more I think about it.  My general thought process is this:

Make sure the text elements I want to use share the name of the domain they will pull from
List layout elements
List domains in the database
For every layout element, step through the domain list until the domain equals the element name
Take the current element text (the coded value) and replace it with the equivalent description from the domain

I got this far in my code before thinking there has to be a more efficient way:

import arcpy
mxd = samplemxd

domains = arcpy.da.ListDomains("Database Connections\\database.sde")

Elements = arcpy.listlayoutelements(mxd)

For elem in Elements:
    for dmn in domains:
        if elem.name == dmn.name:
            coded_values = domain.codedValues
            for val, desc in coded_values.iteritems():
                if val == !layer.attribute!
                    elem.text = desc

Any ideas?  I'm open to non-python ideas as well, in case there's something database side I can do.


texts: Where can I find a shapefile for the coastline of North America?

texts:
I am looking for a shapefile representing the coastlines of North America, such as something like this:
http://www.diva-gis.org/gdata

Does anyone have any suggestions?


texts: Iterating over years for features in public feature collection using GEE?

texts:
I want to calculate the surfaces deforested every year for every protected area in Democratic Republic of the Congo. I am building from the script provided as answer by Rodrigo Principe in reply to this similar question. But this example used a private shapefile loaded in assets. When I try to reproduce this solution with a public dataset (in-land protected areas in DRC, so my question somehow completes the official GEE tutorial on quantigfying deforestation in DRC), it doesn't work : output "areas" remain desesperatelly empty.
The code bellow is the same as proposed by Rodrigo Principe, except the lines marked with an"// EDIT" comment. Does anyone knows how to fix it?

// Load and filter the Hansen data
// EDIT : The source in following line has been updated to load v1_3
var gfc2014 = ee.Image('UMD/hansen/global_forest_change_2015_v1_3')
              // END EDIT
              .select(['treecover2000','loss','gain','lossyear']);

// list for filter iteration
var years = ee.List.sequence(1, 14)

// turn your scale into a var in case you want to change it
var scale = 100

//add country districts as a feature collection
// EDIT: loads public source+filters inland protected areas in DRC
var distr = ee.FeatureCollection('WCMC/WDPA/current/polygons');
// Filter to keep only inland PA in 2 countries
var distr = distr.filter(ee.Filter.and(
    ee.Filter.eq("ISO3", "COD"),
    ee.Filter.eq("MARINE", "0")));
// END EDIT

//look at tree cover, find the area
var treeCover = gfc2014.select(['treecover2000']);
var areaCover = treeCover.multiply(ee.Image.pixelArea())
                .divide(10000).select([0],["areacover"])

// total loss area
var loss = gfc2014.select(['loss']);
var areaLoss = loss.gt(0).multiply(ee.Image.pixelArea())
               .divide(10000).select([0],["arealoss"]);

// total gain area
var gain = gfc2014.select(['gain'])
var areaGain = gain.gt(0).multiply(ee.Image.pixelArea())
               .divide(10000).select([0],["areagain"]);

// final image
var total = gfc2014.addBands(areaCover)
            .addBands(areaLoss)
            .addBands(areaGain)

Map.addLayer(total,{},"total")

// Map cover area per feature
var districtSums = areaCover.reduceRegions({
  collection: distr,
  reducer: ee.Reducer.sum(),
  scale: scale,
});

var addVar = function(feature) {

  // function to iterate over the sequence of years
  var addVarYear = function(year, feat) {
    // cast var
    year = ee.Number(year).toInt()
    feat = ee.Feature(feat)

    // actual year to write as property
    var actual_year = ee.Number(2000).add(year)

    // filter year:
    // 1st: get mask
    var filtered = total.select("lossyear").eq(year)
    // 2nd: apply mask
    filtered = total.updateMask(filtered)

    // reduce variables over the feature
    var reduc = filtered.reduceRegion({
      geometry: feature.geometry(),
      reducer: ee.Reducer.sum(),
      scale: scale,
      maxPixels: 1e9 // EDIT : added to allow large areas
    })

    // get results
    var loss = ee.Number(reduc.get("arealoss"))
    var gain = ee.Number(reduc.get("areagain"))

    // set names
    var nameloss = ee.String("loss_").cat(actual_year)
    var namegain = ee.String("gain_").cat(actual_year)

    // alternative 1: set property only if change greater than 0
    var cond = loss.gt(0).or(gain.gt(0))
    return ee.Algorithms.If(cond,
                            feat.set(nameloss, loss, namegain, gain),
                            feat)

    // alternative 2: always set property
    // set properties to the feature
    // return feat.set(nameloss, loss, namegain, gain)
  }

  // iterate over the sequence
  var newfeat = ee.Feature(years.iterate(addVarYear, feature))

  // return feature with new properties
  return newfeat
}

// Map over the FeatureCollection
var areas = districtSums.map(addVar);

Map.addLayer(areas, {}, "areas")

print(areas)


texts: Geoserver dxf plugin - How to export attributes labels?

texts:
I'm trying to export from Geoserver a set of layers with WFS request and dxf plugin for dxf conversion.

There is an option withattributes:true that I have added for each vector layer a LAYERNAME_attributes, but I don't know how to display labels in my dxf drawing.

I want to see street names, numbers.... and so on

Here is my example

I have attached an image showing the file loaded in DraftSight a CAD software


texts: I'm trying to use GPS clocks to synchronize some measurements and I'm a little confused about the leap second correction between GPS time and UTC. I was wondering - does the GPGGA NMEA header have the 18s leap second offset between GPS time and UTC corrected already? Further, are there GPS NMEA headers that don't have that correction factored in? I'm using two different GPS systems to produce time stamps and they have an 18s offset between each other, indicating that one is correcting for it and the other is not.


texts: Open mdb database (created in ArcGIS10) in QGIS 2.2/2.4

texts:
I have a bunch of orthophotomaps (about 250 tiff files) organized in 3 mdb databases created in ArcGIS10.
I'm now trying to open the mdb files one by one in Qgis 2.2 and 2.4 (Layer->Add vector layer->Browse->Open) but I cannot browse within the database to select the table to import. I can just select the mdb file and I keep getting the same error message "Invalid Data Source: E:\Orthophotos\Block 3\block3.mdb is not a valid or recognized data source".

I know several questions have been asked already on this matter but none of the solution suggested helps me to solve this problem.

Does anyone know how to handle the issue? If QGIS is not able to open such files, I'm available to start from scratch (i.e. from the tiff files) but I don't know how to do it.


texts: Processing a very large hosted KML file to identify what makes it so big

texts:
I've run into a problem where I have two very large layers that I'm tasked to overlay on a map, showing cellular coverage across a the whole South Africa.

Layer 1 was given to me in various formats and it was really easy to work via the Javascript API after getting it into a fusion table.

Layer 2 is my struggle currently, as it was given to me as a hosted KML that is updated automatically every 24 hours. The only reason it's a problem is because of its size, being 90MB!

I've explored several options including:

Downloading and simplifying with ogr2ogr, but it drops the size with about 20MB's, which is still nowhere near ideal.
QGIS - trying to work with the layer hangs the system
pyKML and fastKML - both seem to be capable of reading from it but I'm not sure if this is the best option. I like working in Python and not afraid of a learning curve in parsing the correct data, but I haven't been able to find a clear example of this use case.
Google Earth - Seems to be too big even for GE as it also doesn't show up or hangs the program.

I'd like to know if there's a way to simplify such a big file. Right now I can open it in QGIS and can see Layer 2 there, but like I said, trying to work with it is memory intensive.

Can it be that the KML contains data that wasn't optimized for such a task? I'd like to be able to look at the KML data in a text editor, but even that hangs my Atom editor and others I've tried. I want to see what I'm working with in order to try and remove unnecessary stuff in there.

The url is here: http://siteadmin.wbs.co.za/kml/coverage.kml

I just need a few suggestions on how to handle this.


texts:
I have a raster with null data. I want to create a mask for where ever I have values in my raster. (not null cells)
When I use translate(convert format) in qgis (Raster-->conversion--> translate)to transfer the no data in the original file to 0, instead of changing nodata to 0, it is changing them to -999!

And when I use raster calculator and divide the raster by itself to create mask and then using translate to convert null data to 0,  the tool convert the null values to 3.4028234663852886e+38. How I can fix this problem?


texts: Pyproj inverse transform returns NaN with EPSG that is not 4326

texts:
I am using pyproj inverse transform to add azimuth and distance "info" to an ordered geodataframe and keep getting unexpected NaN results when I use a local UTM-based EPSG.
It works fine in WGS84 (my earlier question here), but projected systems are better than geographic for getting reliable azimuth and distance metrics, which is why I want it to work in other EPSGs.
Below are two examples for two different EPSG codes, and both return NaN for the azimuth and distance fields. Both examples use the same syntax/approach, the point is that the issue is repeatable in two different EPSG codes.
Running the following:
Windows 10
conda 4.8.2
Python 3.8.3
shapely 1.7.0 py38hbf43935_3 conda-forge
pyproj 2.6.1.post1 py38h1dd9442_0 conda-forge
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib.ticker import ScalarFormatter
import contextily as ctx
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from shapely.geometry import LineString
import pyproj
from pyproj import CRS

Example 1
myid = [1, 1, 1]
myorder = [1, 2, 3]
x = [550338.0319, 550428.0048, 550523.9951, 550589.9544]
y = [3795929.972, 3795798.055, 3795659.962, 3795528.029]
myepsg = 32611

df = pd.DataFrame(list(zip(myid, myorder, y, x)), columns =['myid', 'myorder', 'y', 'x'])
gdf_pt = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['x'], df['y']))
gdf_pt = gdf_pt.set_crs(epsg=myepsg)

ax = gdf_pt.plot();
ax.set_aspect('equal')
ax.set_xticklabels(ax.get_xticklabels(), rotation=90);
ax.xaxis.set_major_formatter(ScalarFormatter())
ax.ticklabel_format(style='plain', axis='both', useOffset=False)

geod  = CRS.from_epsg(myepsg).get_geod()
for i, r in gdf_pt.iloc[1:].iterrows():
    myinfo = geod.inv(gdf_pt.x[i], gdf_pt.y[i], gdf_pt.x[i-1], gdf_pt.y[i-1])
    gdf_pt.loc[i, 'az_fwd'] = myinfo[0]
    gdf_pt.loc[i, 'az_back'] = myinfo[1]
    gdf_pt.loc[i, 'dist'] = myinfo[2]
    gdf_pt.loc[i, 'bearing'] = max(myinfo[1], myinfo[0])

display(gdf_pt)

Example 2
myid = [1, 1, 1]
myorder = [1, 2, 3]
lat = [5174925.07851924, 5174890.26832387, 5174855.45812849]
long = [1521631.6994673, 1521667.11033893, 1521672.52121056]
# typo above, it says lat/long but it really is UTM-y & UTM-x
myepsg = 2193

df = pd.DataFrame(list(zip(myid, myorder, lat, long)), columns =['myid', 'myorder', 'lat', 'long'])
gdf_pt = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['long'], df['lat']))
gdf_pt = gdf_pt.set_crs(epsg=myepsg)

geod  = CRS.from_epsg(myepsg).get_geod()

ax = gdf_pt.plot();
ax.set_aspect('equal')
ax.set_xticklabels(ax.get_xticklabels(), rotation=90);
ax.xaxis.set_major_formatter(ScalarFormatter())
ax.ticklabel_format(style='plain', axis='both', useOffset=False)

for i, r in gdf_pt.iloc[1:].iterrows():
    myinfo = geod.inv(gdf_pt.long[i], gdf_pt.lat[i], gdf_pt.long[i-1], gdf_pt.lat[i-1])
    gdf_pt.loc[i, 'az_fwd'] = myinfo[0]
    gdf_pt.loc[i, 'az_back'] = myinfo[1]
    gdf_pt.loc[i, 'dist'] = myinfo[2]
    gdf_pt.loc[i, 'bearing'] = max(myinfo[1], myinfo[0])

display(gdf_pt)


texts:
Hopefully this isn't so far afield as to get banned - apologies in advance to those who think it doesn't belong.

This article discusses a region of Antarctica which is described as "West Antarctica". However, Antarctica is a roughly circular mass centered on the South Pole. Hence any boundary region is north of the center of mass. I suspect there's a convention for naming a region with a compass direction, perhaps relating to the prime meridian or some such, but I can't find information on this. It's just a matter of curiosity, but I would like to know how the region in question comes to be "west".


texts:
I am wondering if it is possible to create scaled labels so that at a certain zoom level only one label appears but when you zoom in closer more labels appear for different fields.

I have this VBScript expression to allow for multiple fields to appear in labels (This is what I would like to see at a zoomed in level e.g at 1:10,000 (my layer consists of large polygons fyi))

Function FindLabel ( [Name], [Kvalue], [Pvalue], [LimeRequirement], [PHvalue] )
  FindLabel = "Plot Name: " & [Name]&vbnewline & "K Value: " & [Kvalue]&vbnewline & "P Value: " & [Pvalue]&vbnewline & "Lime Requirement: " & [LimeRequirement]&vbnewline & "pH : " & [PHvalue]
End Function

Output of expression:

Plot Name: Balally Park

K Value: 1.76

P Value: 5.97

Lime Requirement: 4.37

pH : 5

At a zoomed out level (from 1:500,000 to 1:10,000) I would just like Plot name and K value.

Is there a way of doing this without creating copies of the layer and setting scales on the layer itself?


texts: Edit Field Name In Attribute Table (QGIS 3)

texts:
This seems like such a simple task, yet after researching how to execute it in QGIS I've been unsuccessful.  I'm merely trying to rename columns in my attribute table.  My existing field is currently named "employment", and I would like to rename it "BG_Desc".

I've read the "table manager" plug-in will do this, but it appear the plug-in is no longer available for version 3.0.

I've used the "replace" tool under "String" in the field calculator, but it does not work when I enter the following:

replac("employment",employment","BG_Desc")

Suggesstions?


texts:
I am using GeoServer+OpenLayers and I have a base layer (Google map) and some additional layers (OpenLayers.Layer.WMS); users can show or hide these layers.
I would like when user show a layer then the map should zoom to the area where layer has data.

I know I have to use map.zoomToExtent(bounds) but I don't know how get this information from the layer.

Do you have any suggestion?


texts: ArcGIS find conflict between layer and geometry

texts:
I have a problem I could use some help with. In our application we have multiple FeatureLayers and on some of these layers a user can draw new objects.

When a user draws a new object, lets call it a tent (any type of Geometry), I need to check if the position is valid. The object is valid when it does not consume the same position as a building. The buildings are stored on a separate featureLayer that is stored on a MapServer (and loaded in the browser). The drawing, the tent, is only stored locally (for now).

So to get to the point.. what is think I need is a service that finds the distance between all objects on a layer(the buildings) and the object that is drawn at-the-moment..

I did have a look at the SummarizeNearby services but at this moment we don’t have the license for it. Then I had a look at the GeometryService and I would like to know if this is the correct way to do it..

I think it is possible to use the difference or intersect method of the GeometryService to find the difference between all the graphic/geometries on the buildings-layer and the drawn-object. If the difference is the same size as the graphics in the buildings-layer, then there is no conflict. But if the result is different then the buildings-layer, then the buildings-layer did contain the object that is just drawn, and there-for the drawn object is not valid.

Please let me now if this is the way to go, or if there are other ways to calculate/find objects that intersect/overlap.


texts: Openlayers zoomToExtent from array of coordinates

texts:
I have a dropdown list. On Selection, it will request extent of selected property to database (PostgreSQL). I am using PHP for server-side processing.

At client-side on console log gives extent output result (variable res):

[74.3965356273929, 75.2484319931094, 35.6450698628058, 34.776235256942]

Also I am sure that coordinates are float and not string using parseFloat() function. later I am  passing that variable to OpenLayers.Bounds.

    var extent = new OpenLayers.Bounds(res);
    map.zoomToExtent(extent);

OpenLayers.Bounds gives following output:

OpenLayers.Bounds.OpenLayers.Class.initialize {left: 74.396535627393, bottom: 75.248431993109, right: 35.645069862806, top: 34.776235256942, left: null…}

Which is I guess correct. But Openlayers is not zooming to the provided extent and it is zooming to maximum zoom level at current map view only.


texts: Performing raster calculations with more than two bands in PostGIS?

texts:
We need to perform raster calculations using more than two bands in PostGIS.  At the moment, ST_MapAlgebra has a two-band form but doesn't support more than that.

I'm open to any solution others would like to offer, but am currently favoring writing a custom plpython function that might be called as follows

SELECT
    ST_FromGDALRaster(
        some_custom_function(
            ST_AsGDALRaster(
                some_raster_column,
                'GTiff'
            )
        )
    ) new_raster_column
FROM some_table

where the custom function accepts and returns a bytea PostgreSQL datatype.  Inside the function I am able to import gdal and use following

vsipath = '/vsimem/buffer'
gdal.FileFromMemBuffer(vsipath, gdal_rast)
rast = gdal.Open(vsipath)

gdal_rast above is the parameter being passed into the custom function resulting from the ST_AsGDALRaster function in the select clause.  From there, I can turn the bands into numpy arrays and process as desired.

My problem comes in returning the processed raster back from the custom function to the database.  I have tried the MEM driver to create an in-memory raster, but when I return that raster to the ST_FromGDALRaster function from the custom function I receive this error:

ERROR: RASTER_fromGDALRaster: Could not open bytea with GDAL.
Check that the bytea is of a GDAL supported format

Can the plpython function approach be salvaged?  Or is there a better solution that avoids plpython altogether?


texts:
Check this straightforward code that exports the subset of a shapefile (path_of_shp) using where_clause:

where_clause = '"new" = 1'
arcpy.MakeFeatureLayer_management(path_of_shp, "feature")
arcpy.SelectLayerByAttribute_management("feature", "NEW_SELECTION", where_clause)
arcpy.Copy_management("feature", "C:/GIS/output.shp")

This gives an error:

  arcgisscripting.ExecuteError: Failed to execute. Parameters are not
  valid. ERROR 000840: The value is not a Data Element. Failed to
  execute (Copy).


texts:
I got this error, RuntimeError: wrapped C/C++ object of type QLineEdit has been deleted and
this one RuntimeError: wrapped C/C++ object of type LeastCostPathAlgorithm has been deleted.

I think that is the reason that I cannot execute many tools such as Fill Sinks (Wang and Liu) but I am not sure. For example when I use the Fill Sink tool the parameter for Fill Sink is executed but the flow direction parameter and watershed basins are not.

I was trying to understand what is going wrong. So when I went in settings-options-scripts and chose Script Folder and I click OK I got those errors.

I think that some operations cannot be executed because there is problem with python but I am not sure if I am correct.

Do you know something about it?
Is there a relationship in terms of the errors when I run the tools that I mention?
Do you know how can I fix it?

I use the last version of QGIS 3.10.1 'A Coruna', actually I downloaded the OSGeo4W Network.


texts: Exporting points according to a route I made

texts:
So I work on a ship and get given a shape file with points on it and I am required to make a vessel route out of this.
The problem is that once I have created a route between all the points and export them as a shape or csv, they are not arranged according to the route, but the way they were sent. I then have to manually arrange them in Excel.
For example: I am sent the points in order of A C D B and the route goes A B C D, however when exported it is still A C D B.
I assume I have to link the point and line feature somehow, but not sure how.
Software on the vessel is QGIS.


texts:
I have been searching, but have not found a way to do this...

I have several 'generic' layers that I want to show in my Control for users to toggle on and off... That works fine.

However, I would like to let the user select additional layers from a list, and have them toggle on and off as well. BUT, I do NOT want those layers added to the control. Instead, I want to add them to an external list (shown outside the map), and have a checkbox for each that the user can use to show/hide that layer on the map.

I want to do this via Angular/JavaScript.

Potentially, the user could have 20+ additional layers selected, so adding them to the control is not good, as the list would be very long.

Any ideas?

And no, I do not have any code to share, as I cannot find any way to add a layer to the map, without adding it to the control....

--Sam


texts: Set web container with timezone properly configured on GeoServer

texts:
I am learning how to implement time series with raster on GeoServer following the image mosaic documentation .

I am struggling a little bit to understand what this sentence means:

  To support time-series layers, GeoServer needs to be run in a web
  container that has the timezone properly configured. To set the time
  zone to be Coordinated Universal Time (UTC), add this switch when
  launching the Java process:

  -Duser.timezone=GMT

  If using a shapefile as the mosaic index store (see next section),
  another java process option is needed to enable support for timestamps
  in shapefile stores:

  -Dorg.geotools.shapefile.datetime=true

How could I set those?
Maybe when I start GeoServer?
I am running GeoServer on a docker container.


texts:
I'm using Win10 (with the bash shell) and one of my python scripts requires gdal/osgeo.

I ended up trying multiple install methods (gdal was successful with all of them, but I couldn't use osgeo), but the last one was Anaconda.

But when I try and run the script, I get:

from osgeo import ogr
ImportError: No module named osgeo

And gdal info yields:

Usage: gdalinfo [--help-general] [-json] [-mm] [-stats] [-hist] [-nogcp] [-nomd]
                [-norat] [-noct] [-nofl] [-checksum] [-proj4]
                [-listmdd] [-mdd domain|`all`]*
                [-sd subdataset] [-oo NAME=VALUE]* datasetname

FAILURE: No datasource specified.

I did some digging and found that:

a) my GDAL folder doesn't have Python in it (I don't know if this is necessary, but some other questions I looked at seemed to specifically mention C:/GDAL/Python

b) in my C:/Python27/Lib/site-packages folder, I found an osgeo folder.

At risk of being terse: how do I make it work? (If this is a duplicate or there's a similar post that I missed, please link it)


texts:
Using ArcGIS 10.3.

The .kmz files show up as saved to the appropriate location in the shared drive of the computer. The drive appears to be "talking" to the GIS catalog fine; I can create a folder within the catalog and it will appear in the computer drive. I cannot, however, populate it with files that will then show up in the catalog.

It has been refreshed, and restarted, to no avail. The files work fine in Google Earth, and I have had no issue converting between the two in the past.


texts: Rotating map canvas for QGIS 2.8.3 atlas generator?

texts:
I have a atlas coverage layer that follows a water network for use with the QGIS atlas generator. I am having issues getting the map canvas rotation to update dynamically with each coverage feature.

What I've done:

set up a print composer with the map added
On the generate atlas tab, I checked the box to generate an atlas and chose my coverage layer
On the map object properties, I checked Parameters controlled by atlas, with 10% margin
Under map rotation, I clicked the expression drop down, chose edit and entered "Angle" (the name of my field with an angle of rotation).

After generating the atlas and clicking through coverage features, the scale updates, but not the rotation (stays 0). Is there maybe something wrong with my expression?


texts:
I would like to transform SDO_GEOMETRY from WGS84 projection to my custom projection in Oracle Database.

Details for custom projection is as follows:

+proj=lcc +lat_1=40.5 +lat_2=37.5 +lat_0=25 +lon_0=36 +x_0=1003827.11 +y_0=-1183453.08 +datum=WGS84 +units=m +no_defs

PROJCS["Customized Lambert",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",40.5],
    PARAMETER["standard_parallel_2",37.5],
    PARAMETER["latitude_of_origin",25],
    PARAMETER["central_meridian",36],
    PARAMETER["false_easting",1003827.11],
    PARAMETER["false_northing",-1183453.08],
    UNIT["Meter",1]]'

Running SDO_CS.VALIDATE_WKT for custom projection returns TRUE.
But,

SELECT SDO_CS.TRANSFORM(SDO_GEOMETRY('POINT(29 41)',8307),1000003006) FROM DUAL;

returns with ORA-13282 & ORA-06512 errors.

What is the problem here ?


texts: Web API for reverse geocoding geographic features

texts: Is there a free web API that can reverse geocode geographic feature names? I'm currently using the Google API, which will return political and street address info, but I would like to return the names of mountains, islands and especially bodies of water for a given lat/lng pair.


texts:
I am loading an nc4 raster downloaded from a public dataset to QGIS and it first displays an error that the data does not have CRS. Then, it charges the map but the reference point 0 is shifted the map starts in Europe (I see most of the maps start in America on right side). Therefore, when I add an openstreet map layer, both maps don't match. I have tried changing the CRS but when I do that the map appears cropped.
I attach an image of the problem


texts: How to correctly create a polygon from an extent?

texts:
I'm tying to create a new polygon based on some calculated extent. I'm doing this way, as suggested by documentation:

let newPolygon = new Polygon().fromExtent(newExtent);
let newFeature = new Feature({
    geometry: newPolygon,
});

This is giving me an error in the new Polygon() line:

  Cannot read property 'length' of undefined.

newExtent is an array and is in accordance with the documentation as well.


texts: Zooming in after setting coordinate system in ArcMap?

texts:
I'm using ArcGIS 10.3 and after setting up the coordinate system (WGS 1984) I added some georreferenced *.JPG files

However when I try to zoom in/out to digitize, the zoom goes crazy: it doesn't display the image correctly and I'm not able to controll the zoom anymore, I just can visualize the overview by 'Zooming to Layer'. .

I tried to convert the image into a TIFF file but I'm getting the same bug.

Any hint on what should I do to solve this problem?


texts: R, given a set of center points how to define the cell area?

texts:
I have the coordinates for the police stations in a city, they are grouped close together at the center and sparse towards the suburbs.

Since I do not have the shapefile for the jurisdiction of each station, I want to do a rough estimate of the area by defining non overlapping variable size and dimension polygons with the police station as the center.

Is there a standard function in R or any other open source toolkit which will do this?

If a standard function is not available what algorithms is best suited for this task?


texts: QGIS Raster to Vector Conversion - 3.4.6-Madeira

texts:
I have some contour lines in raster format, I want to convert them vector. I tried with QGIS (version 3.4.6) I got the vector form.

But when I tried to select the contour lines to put the actual values those lines represent, I noticed that these lines are formed by rectangles as shown below.

Now the problem is I can't select the entire contour and put the required values. How can I get a line instead of these?


texts: Parameter dependencies based on two datasets in Python toolbox

texts:
I am setting up the parameters in a Python toolbox. I have a parameter that I would like to list fields based on a user input feature class. I am wondering if it is possible to list fields based on two feature classes.

def getParameterInfo(self):
    input_data = arcpy.Parameter(
        displayName = "Input data to be fuzzed to a grid",
        name = "input_data",
        datatype = "GPFeatureLayer",
        parameterType = "Required",
        direction = "Input")

    existing_grid = arcpy.Parameter(
        displayName = "Existing Grid",
        name = "site_desc",
        datatype = "GPFeatureLayer",
        parameterType = "Optional",
        direction = "Input")

    output_fields = arcpy.Parameter(
        displayName = "Fields to be Included in Output Fuzzed Data",
        name = "output_fields",
        datatype = "Field",
        parameterType = "Required",
        direction = "Input",
        multiValue = True)
    output_fields.parameterDependencies = [input_data.name]

I would like for the user to be able to also choose fields from the existing_grid dataset, in addition to the input_data dataset.

How do I include the fields from both in the field list?

I tried adding it separately and including it in the brackets, but this just caused the tool parameter to be blank.


texts: Field properties in QGIS are lost when polygon layer shape file is reloaded

texts: I have been creating Shapefile polygon layers.  The Attributes Form is used to create a Value Map in the 'Widget Type' section and everything works fine.  The Value Map option creates a combo box pull down of options.  The layer is saved off as a shapefile, but when it's reloaded, the 'Value Map' settings are lost along with any default settings for other fields.  Is there a way to preserve these settings?


texts: MultiLine field-name label referencing in PyQGIS?

texts:
A labeling routing (which is updated often) is being implemented via pyqgis and a QComboBox.

The code below worked great, that is until needs for narrowed column headers spawned the use of pythons '\n' (and subsequent problems thereafter).

This string formatting switch has caused any multiline field reference to return blank labels, while neighboring fields with single lines are displayed correctly.

A remedy I initially sought was to use fieldname (as strings) for the label settings, but this did not work. Then it was supposed that a index based context would work better (avoiding syntax), but such did not solve the problem

I am suspicious that I might be missing a line in my publishing of the QgsPalLayerSettings() class to the layer.  Because when I do it manually through the QGIS labeling menu, multiline field names work perfectly fine.

Related to a previous question:
Connecting vector label size with map-units?

field = str(self.csvGridFieldNames[index])
activeLayerName = self.iface.activeLayer().name()
csvGrid = ftools_utils.getMapLayerByName(unicode(activeLayerName))
# fieldIdx = csvGrid.fieldNameIndex(field)
textSize = 'coconut hoves'
valLabel = QgsPalLayerSettings()
valLabel.readFromLayer(csvGrid)
valLabel.enabled = True
valLabel.placement = QgsPalLayerSettings.OverPoint
valLabel.fieldName = field
# valLabel.fieldIndex = fieldIdx
valLabel.fontSizeInMapUnits = True
valLabel.setDataDefinedProperty(QgsPalLayerSettings.Size,True,True,'%f' %(textSize),'')
valLabel.writeToLayer(csvGrid)
valLabel.writeToLayer(csvGrid)
iface.legendInterface().refreshLayerSymbology(csvGrid)
iface.mapCanvas().refresh()

For sake of missing information, the following is a list of a few fieldNames (as contained in the list 'self.csvGridFieldNames');

While 'Status' and 'Kx' publish labels correctly, both 'Top\nElevation','Bottom\nElevation' do not.

  ['Status','Top\nElevation','Bottom\nElevation','Kx']

Per advice, the following also does not work:

(And 'print repr(field)' returns expected 'Top\\nElevation')

    rawField = str(self.csvGridFieldNames[index])
    if '\n' in rawField:
        beginingOfString,endOfString = rawField.split("\n")
        field = beginingOfString+'\\n'+endOfString
    else: field = rawField
    valLable.wrapChar = '\\n'  #<<=== with/without this line, still blank labels

Link to a ultra-beta QgsPointLayer Labeler pulled from a larger (yet unpublished) plugin.
Extract contents to Qgis Python Plugins Dir as normal.

ReadMe instructions are a must for functionality.

QgsPointArray Labeler


texts: How do I delete a point (not on a polygon)?

texts:
I want to delete points that is not part of a polygon. I can add points using Capture Point, but I can't delete them. The following solution does not work for points:

Select your layer
Enabled Editing
Choose the Node Tool
Click a point on your polygon
Use Command-click & drag to select multiple points with the marquee tool. (might be shift-click on Windows)
Use Delete key to delete the selected points (Fn+Delete on MacBook Pro)


texts: Working with Raster data/DEM (mosaic to new raster, viewshed) in ArcGIS Desktop?

texts:
I'm trying to work with raster data from this source: http://www.soest.hawaii.edu/coasts/data/maui/index.html

I downloaded the DEMs from different islands and imported them into ArcGIS. But I can't seem to do anything with them. First, I tried to merge the different DEMs into one new raster. I also can't calculate a viewshed from a known point. I constantly get errors from the like of "dem not valid".

How do I convert them into formats I can actually work with?

I'm not too familiar with raster datasets.


texts: Leaflet plugins - how to use two plugins together

texts:
I want to use two leaflet plugins: ExtraMarkers and IconLabel.
How to combine them together? The result should be ExtraMarker with label.

When I create ExtraMarkers:

        var icon = L.ExtraMarkers.icon({
            shadowSize: [0,0],
            prefix: 'fa',
            extraClasses: 'fa-fw',
            shape: 'circle',
            icon: 'fa-'+Icon,
            iconColor: IconColor,
            markerColor: MarkerColor
        });

or IconLabel

        var iconLabel = L.Icon.Label.extend({
            options: {
                iconUrl: 'img/icon.png',
                shadowUrl: null,
                iconSize: new L.Point(24, 24),
                iconAnchor: new L.Point(0, 1),
                labelAnchor: new L.Point(26, 0),
                wrapperAnchor: new L.Point(12, 13),
                labelClassName: 'icon-label'
            }
        });
        var icon = new iconLabel({labelText: title})


texts: If two columns match, return only closest absolute value from other two columns

texts:
I have a shapefile with points that are joined to the closest twenty points from another shapefile. The joined shapefile has columns from the original shapefile: name (ID), elevation (Z), name from joined shapefile (name) and elevation from the joined shapefile (Z2).

So twenty rows have the same ID and name but differing Z and Z2. I would like to use the field calculator to return the row that has the smallest absolute difference between Z and Z2 with the matching ID and name.

Using a case statement in the qgis field calculator I can get all the results to return the difference in a new column:

CASE WHEN "Name" = "ID" THEN abs("Z" - "Z2") END

but I would like only the smallest value to return.


texts: Set buffer radius in metric units in qGIS

texts:
Possible Duplicate:
  What is up with QGIS buffer tool units?

A similar question to this one has already been asked How to change the unit of measure from degrees to meters in QGIS?, but none of the steps suggested has worked for me.

I have a data set with over 37000 points. I would like buffer the points with 5km, 10km, and 15km radii using f-tools, dissolving the buffers to see how the total surface under the buffers evolves.

By default, my qGIS has the CRS set at WGS84 (settings-> project properties). My shape file with the points was also initially set at WGS84.

As I mentioned, I would like to draw the buffers with the radii set in metric units, however, the units are always in degrees.

I have tried adjusting my qGIS project settings both before and after opening the shape file, modifying the shape file settings (always changing both to EPSG:3414 - SVY21 / Singapore TM as this is data for Western Malaysia) directly, etc.

Somehow, I always get buffers in degrees rather than metric units.


texts: duplicate polygons once for each point inside

texts:
I have a set of feature classes: points representing addresses, and polygons representing buildings. Sometimes there are more than one addresses per building.

I use Maplex to generate and place labels. Instead of labeling the points with the address, I would like to take advantage of the Maplex "building labels" features to place the address labels nicely within the polygon boundaries.

I need to create a polygon for each of the points, using the geometry of the polygon that the point is within, with the attributes of the point.

Thus far I have able to use the spatial join tool to create a set of polygons with the attributes of one point within them, or a set of all the points with the attributes of the polygon they are within.... neither solution is quite what I'm looking for.


texts: Can I draw GRAPHICS (not vector) like in ArcMap and convert it to shape?

texts:
In ArcMap one can draw graphics (lines, rectangles, circles, etc) that are NOT vectors but can be then converted to shape files. Is same possible in QGIS?

I have only found to make directly vectors in edit mode.


texts: Redrawing WMS layer on OpenLayers without flickering from GeoServer?

texts:
I have a WMSlayer that is updated every 10sec with new data from SQL.
Whenever I redraw the layer using redraw(true) , the map screen flickers.

Is there a way to do it without flickering ?


texts: How do I import data from a delimited file into a Feature Class?

texts:
I have a feature class created in a file geodatabase. How do I import data from a comma-delimited file into the feature class?

Update #1 - The destination is a feature class, not a table.


texts: Google map tile size and number for an area

texts:
Is there any way to estimate the number of tiles in an area for Google Maps for each zoom level? I need something like this app:

http://tools.geofabrik.de/calc/#type=geofabrik_standard&bbox=-14.885271,8.083153,7.614729,27.393223

but for Google Maps. I also think it would be possible to calculate using this equation for whole earth:

sigma(i^4)

which i goes from 1 to 19 (zoom level) but it calculates for the whole world and gives tile numbers, I can later multiply it into average tile size (I have no idea about this value too).


texts: Difference between radiometric correction at surface and at sensor in ENVI

texts: I want to ask you about the difference between the radiometric correction at surface and at sensor in ENVI software, I do the two in flash and radiometric correction and the results was same image. Can you help me to understand?


texts: How to use GRASS module to delineate watershed area？

texts:
Update#1：I changed to work under qgis 1.8, I only input the DEM file and run the r.watershed module, and I got this：

I have a DEM file of one specified region, and I take the following steps to get files that might be used：

Fill Sinks(this one I use SAGA module instead of GRASS)
r.fill.dir

After step 2, I want to get flow accumulation file, but I can not find this kind of module in GRASS, so I pass it and directly try the r.watershed module.

There are several input options in this module, like this：

I tried any of files I have as input files, but I can't get the output files through this module, and always got this error message：

The log below told me that the data source is invalid, and I go to check the sextante message log：

ERROR|Wed Jan 22 2014 12:38:01|Could not load model watersheds.model|ALGORITHM:saga:catchmentarea(parallel)

What should I do to successfully use the r.watershed module？

Thanks for giving me any suggestion!


texts: Determine sewer tap elevations along a sewer line shapefile

texts:
I have a sewer collection line shapefile with upstream and downstream invert elevations (Up_el and Dwn_el). I have created a point shapefile of service connections based upon the endpoints of a private lateral line shapefile, using the data management tool "feature vertice to point" (I used endpoint).

I suppose that I could split each sewer line at the service connection point (sometimes called a "tap") and calculate the y value using the individual line slope and the distance the tap is down from the upstream manhole, but I am hoping there is a more eloquent tool or code to populate a tap elevation field for this point shapefile attribute table.


texts: Deleting features with Null value in attribute field

texts:
I have polygonized a single band raster and created an area field. I want to keep the 10 largest polygons by area and delete the rest. The following script computes the area for all polygons and sorts the area field and keeps 10 largest values and assign NULL to rest of the features.

img = r"C:\Users\test\input_img.tif"

# mapping between gdal type and ogr field type
type_mapping = {gdal.GDT_Byte: ogr.OFTInteger,
                gdal.GDT_UInt16: ogr.OFTInteger,
                gdal.GDT_Int16: ogr.OFTInteger,
                gdal.GDT_UInt32: ogr.OFTInteger,
                gdal.GDT_Int32: ogr.OFTInteger,
                gdal.GDT_Float32: ogr.OFTReal,
                gdal.GDT_Float64: ogr.OFTReal,
                gdal.GDT_CInt16: ogr.OFTInteger,
                gdal.GDT_CInt32: ogr.OFTInteger,
                gdal.GDT_CFloat32: ogr.OFTReal,
                gdal.GDT_CFloat64: ogr.OFTReal}

ds = gdal.Open(img)
prj = ds.GetProjection()
srcband = ds.GetRasterBand(1)

# Create shapefile layer
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource(r"C:\Users\test")
srs = osr.SpatialReference(wkt=prj)
dst_layer = dst_ds.CreateLayer("output_shp", srs=srs)

# Create area field
col = ogr.FieldDefn('AREA', ogr.OFTReal)
dst_layer.CreateField(col)

# Compute area and sort values
features = []

def get_area(feature):
        area= feature.GetGeometryRef().GetArea()
        return area

for feat in dst_layer:
    features.append(feat)

# Sort 10 largest values
sorted_areaVals = sorted(features, key= get_area, reverse=True)[:10]

# Write sorted values to the filed
for sv in sorted_areaVals:

    sv.SetField("AREA", get_area(sv))
    dst_layer.SetFeature(sv)

Area field gets the 10 largest values while rest of features get NULL.
How can I  delete the features with NULL in area field?


texts: Joining lots of small polygons to form larger polygon using PostGIS?

texts:
I have the following layer using SRID 27700 in postgis:

It's every administrative region in the UK, and (as you can see from the colour grouping) each of them has a text field specifying the county they lie in.

What I'd like to do is to make larger county polygons from the smaller ones in a given county, so EG in the picture above all the teal colour polygons would form one large polygon from the single outer ring that contains all the polys in that colour, like wise all purple, brown, pink, grey etc should all form one polygon.

I've already tried the following:

insert into parishesmerged (geometry)
select astext(multi(ST_Union(the_geom))) as the_geom from parishes
group by county_name

But it keeps generating broken geometries which I then have big problems processing further.

I'm trying to make a simpler county level map with the major output areas in.

Any solutions don't have to be in Postgis either, I have the full OS4Geo stack installed, the latest version of QGis and more utils than I can shake a stick at.

The only things I don't have are the big boys like ArcGis (Although I may have an Old Mapinfo lying around somewhere)

For the record, the dataset I'm trying to create is to accompany a GIS book I'm on writing aimed at .NET programmers wishing to write GIS applications using .NET

After trying the suggestions below, the one that worked the best was 'Paul Ramseys' solution.

I now have a nice simplified counties & boroughs file that's just simple enough for my book, but complex enough to allow me to demonstrate some interesting geo-spatial SQL.

Even though Paul's solution ultimately was the one that worked for me, I also drew on the other answers for things like simplifying the polygon map and reducing the complexity further.

On thing I did observe while doing this however, while ST_Collect is indeed faster than ST_Union, run for run it was also the one mostly responsible for broken geometries.  My guess is the speed increase is obtained at the expense of less accuracy in the core function.


texts:
Is there a way to connect directly form QGIS to a Google spreadsheet?

The purpose for pulling data from a Google spreadsheet is diverse. This could be to pull coordinate values from a questionnaire or content for atlas making.


texts: Projecting features to Enterprise Geodatabase in Oracle gives ERROR 000210: Cannot create output

texts:
I am trying to project features from a geodatabase and write into an enterprise geodatabase connection. I have permissions to create data and the name I am using also works, but I keep getting error:

  ERROR 000210: Cannot create output

It is being loaded in sde:oracle11g database. I am facing issues with only point feature classes

My code:

in_dataset = r'R:\50_Scratch\4AnnaArora\MigrationProject\test\test_data.gdb\Test'
out_dataset = r'R:\50_Scratch\4AnnaArora\MigrationProject\test.sde\WS_Test'
out_coor_system= 4326
sr = arcpy.Describe(in_dataset).spatialReference
arcpy.Project_management(in_dataset, out_dataset,out_coor_system, sr)


texts: Compiling QGIS 2.0.1 on Centos (SIP problems)

texts:
Trying to compile QGIS 2.0.1 on CentOS 6.4 (64 bit):

/usr/share/sip/PyQt4/QtGui/qsyntaxhighlighter.sip: In function ‘PyObject* meth_QgsExpressionHighlighter_setCurrentBlockUserData(PyObject*, PyObject*, PyObject*)’:
/usr/share/sip/PyQt4/QtGui/qsyntaxhighlighter.sip:78: error: ‘qtgui_wrap_ancestors’ was not declared in this scope
make[2]: *** [python/CMakeFiles/python_module_qgis_gui.dir/gui/sipguipart3.cpp.o] Error 1
make[1]: *** [python/CMakeFiles/python_module_qgis_gui.dir/all] Error 2
make: *** [all] Error 2

This is with SIP 4.15.2, Python 2.6.6, and  PyQt-x11-gpl-4.10.3

Any help please?


texts: What raster-based GIS systems really work?

texts:
Raster GISes like GRASS, ArcGIS/Spatial Analyst, and Idrisi can perform a rich set of data processing and analytical procedures loosely known as "map algebra."  In today's computing environment it is becoming commonplace to maintain rasters of 100,000,000 cells or more in many different formats and to demand relatively complex calculations like viewsheds, watersheds, and terrain identification, as well as image processing capabilities.

It seems that many open source, free, and inexpensive solutions are out there.  But which ones really hold up in practice?  That is, which ones can handle large grids efficiently, can easily get data in and out, are reasonably bug free, and offer a full complement of analytical procedures?  What are the pitfalls or hidden limitations that you don't find out until you have invested lots of time in learning these systems?  (This last question is one that is not easily answered with Web searching and where I am hopeful respondents can offer valuable advice.)

I am especially interested in solutions that can both integrate well with and compete with the popular (but expensive) commercial systems (which means Windows compatibility is important).


texts: Offline Map Service with Mapnik and Leaflet

texts:
I'm relatively new to GIS, so I hope my asking doesn't appear to be ignorant in any way.

I currently work on an offline map service solution for usage in areas without internet connection. For that purpose I want to download the data from a provider (i.e. geofabrik.de), load it into a PostGIS enabled database by using osm2pgsql, process the data with Mapnik and store the rendered tiles on the local machine. The map data doesn't need to be up-to-date but should be able to get updated once or twice a year.

From there, Leaflet will be used to display the rendered tiles and also satellite tiles if possible, as a slippy map. The map user should also be able to edit data on the map and to save them as vector layer to the database.

I already did something similar by using Geoserver, GeoWebCache and OpenLayers but I'd like to have at least two approaches for comparison.

My questions are:

Where is the difference in using Mapnik or Geoserver? This may seem to be a redundant question but I couldn't answer it - so far.
Where do I get the satellite tiles? I read something about OSM providing servers for them but only give out higher capacity downloads to chosen projects.
Are there any performance issues I should look out for? The provided systems are old Panasonic Toughbooks with limited CPU and RAM. I already tried the Geoserver way on them. It worked, since there is no real-time rendering. But maybe I'm overlooking something.
Any flaws with the approach?

Thanks for investing time in this!


texts: How to create a rule on QGIS to show only features that have last week as date?

texts:
Basically I have a PostGIS layer with a date column.

I only need to show features that have the last week as a value in the date column. How to achieve this?

Unfortunately now() - 1 week doesn't seem to work


texts: Attribute Joining with .shp and .dbf showing NULL values in ArcMap?

texts:
I have a shapefile that only has 2 fields, one of which is the same as one of the fields in the .dbf file. I tried joining these 2 by their common field but after I do so, the 10 other fields in the .dbf get assigned NULL values.

I am using ArcMap 10.5.

Some additional context since I've already done some Google searching:

This was initially an Excel spreadsheet but I changed it to a .dbf
They are both strings/text (though they are different lengths which I've seen could cause a problem but there doesn't seem to be an easy solution around this)


texts:
I am creating/drawing a new polyline successfully; the only problem is it's going under my other graphic layers. Is there any ability to set the z-index of this to bring it above the over layers? I read in 4.x there is this ability with 'collections' but I find the documentation very vague.

  var graphicsLayer = new GraphicsLayer();
  map.add(graphicsLayer);

  var simpleLineSymbol = {
      type: "simple-line",
      color: [13, 121, 190, 1],
      style: "short-dash",
      width: 5
  };
 ..............

    var polyline = {
        type: "polyline",
        paths: [
            [lineY, lineX],
            [coordinatesBx, coordinatesBy]
        ]
    };
    var polylineGraphic = new Graphic({
        geometry: polyline,
        symbol: simpleLineSymbol
    })


texts: Exporting layers from ArcGIS into PostGIS while preserving attributes such as color

texts:
I have managed to import/convert a KML into ArcGIS 10.0 (ArcMap) using the built-in conversion tools and now have a series of layers. Next, I export these layers to a PostGIS database using the following command:

ogr2ogr -f "PostgreSQL" PG:'host=localhost port=12345 dbname=mydb user=removed password=removed' Roads.gdb -progress --config PG_USE_COPY YES

Note that when doing so, four separate layers are stored as a single layer in PostGIS - I'm unsure whether or not that is relevant. Anyhow, next, I use GeoServer with the PostGIS database as a store to serve the layer using WMS, and have a Leaflet client display the tiles overlaid on various other tiles.

Somewhere in this process all color attributes from the data are apparently removed (The data are a set of lines over roads, colored based upon speed limit - each color is a separate ArcGIS layer).

How can I preserve the attributes of my data? I have access to ArcGIS Desktop, PostGIS,  OGR/GDAL, and GeoServer.


texts: Geometry generator line from centroid to edge of irregular polygon

texts:
I have a point layer styled with the following geometry generator formula for an irregular polygon. It is the union of (a) a circular buffer with (b) a tree canopy spread buffered by 1m (where that extends beyond (a)).

 union(buffer($geometry,"tpzbsrad_m"),buffer(make_polygon(make_line(
 translate($geometry,0,max("N",0.2)),
 translate($geometry,0.6*max("E",0.2),0.6*max("N",0.2)),
 translate($geometry,max("E",0.2),0),
 translate($geometry,0.6*max("E",0.2),-0.6*max("S",0.2)),
 translate($geometry,0,-max("S",0.2)),
 translate($geometry,-0.6*max("W",0.2),-0.6*max("S",0.2)),
 translate($geometry,-max("W",0.2),0),
 translate($geometry,-0.6*max("W",0.2),0.6*max("N",0.2)),
 translate($geometry,0,max("N",0.2)))),1))

I would like to use the geometry generator to draw a 45° line from the tree point ($geometry) to the edge of this irregular polygon.

For most cases the following suffices - a line to the edge of the circular buffer.

make_line(
 $geometry, make_point($x + "tpzbsrad_m"*cos(radians(45)),
                       $y - ("tpzbsrad_m"*sin(radians(45)))
))

But where the tree canopy extends beyond the buffer in that area it doesn't match up. See image below:

What expression can I use to get the line to be long enough to always touch the edge of the irregular polygon at 45°? shortest_line() results in angles all over the place, naturally.

(Looking for an expression that ideally works in QGIS 2.18 and 3.x)


texts: Converting area from degrees to square kilometres using OpenLayers geometry?

texts:
I get a multypolygon from my web service and want to calculate the area of that polygon. I was trying to use this piece of code:

var poly = new OpenLayers.Geometry.MultiPolygon(poligon);
var area = geom.getArea();

But that returns the area in square degrees.  The projection is EPSG:3243.
I then tried to use

OpenLayers. Layer. SphericalMercator

But that didn't help me a lot.

Is there any other solution or some other hint?


texts: I added the shapefile from thematic layers but I am not able to view the imported shapefile in the project window on the image.


texts: ArcGIS: Get layer by name using Javascript API and modify it

texts:
I'm very new to ArcGIS and I'm developing a custom widget in the ArcGIS WebApp Builder that will allow me to change the style attribute of a specific layer in my map.

To manually change the style attribute I would open the map, hover over a layer and click 'Change style':

Then I would select the style attribute from the list that appears:

At the moment my JavaScript widget allows the user to pick from a list of (manually entered) layers and style attributes:

The function that is called when the user clicks 'Update' resembles the following:

_updateAttribute: function() {
    var lyr = document.getElementById("layer");
    var att = document.getElementById("attribute");
    var selectedLayer = lyr.options[lyr.selectedIndex].value;
    var selectedAttribute = att.options[att.selectedIndex].value;
    // Pseudo-code;
    // map.layer(selectedLayer).Attribute = selectedAttribute;
}

Is there a way to carry out the operation represented by my pseudo-code?