texts: How to create a grid of pixel coordinates from the corners of a MODIS tile?

texts:
I have an HDF4 file whose StructMetadata.0 contains the following attributes:

UpperLeftPointMtrs = (-20015109.354000,1111950.519667)
LowerRightMtrs     = (-18903158.834333,0.000000)

These are X and Y distances in meters of the MODIS Tile for L3 Gridded product (Sinusoidal Projection). I want to extract/create the coordinates of all the pixels (240 x 240) in this tile given the pixel resolution is 5km. How can I acheive this in Python?


texts: pgRouting: unable to execute routing_core.sql

texts:
I'm unable to execute the routing_core.sql, receiving following error.

I'm using pgAdmin III, Windows 7, PostgreSQL 9.1, PostGIS 2.0.

ERROR:  could not load library "C:/Program Files/PostgreSQL/9.1/lib/librouting.dll": %1 is not a valid Win32 application.

*** Error ***

ERROR: could not load library "C:/Program Files/PostgreSQL/9.1/lib/librouting.dll": %1 is not a valid Win32 application.
SQL state: 58P01

Regards,
Nas


texts: How to find polygons inside polygon with PostGIS 2.0

texts:
I have a relation representing all the counties in my country and the other one containing all the districts - e. g. each and every county is composed of several districts.
How do I select them with spatial query?
I tried with

select nazorp from kraje as k, orp_wgs as o WHERE ST_Within(o.geom,k.geom) AND
k.nazev = 'Liberecký'

but had wrong results returned (well, they are probably not wrong, they are just not what I expected them to be). ST_Within only returns districts that don't share a boundary with the county, but I need to get all the districts within the county. Is that possible?
I haven't found any built-in function suitable for my needs yet.


texts: Creating dynamically styled symbols (fill, stroke, stroke-width ..) with GeoServer?

texts:
Case scenario:

I need to create symbology that takes fill, stroke and stroke width with information from the database. Problem arises with custom defined symbols like the one on the image. There are lot of them.

One of the solution that came to my mind was to create font and style that font. But that approach is not possible because I can only set fill to the font, not stroke and stroke width. Next idea was to use SVG and somehow try to dynamically set fill, stroke and stroke width. But I had no success. I tried the approach at How to create svg symbols that have modifiable fill color, stroke color and stroke width? but with no success.

Is there a way to do this using GeoServer?


texts: QGIS Exponential Notation in Raster Calculator

texts:
I am looking to execute an equation in the QGIS raster calculator using multiple instances of exponential notation. The current format I am using is as follows:

.00005*exp (x,y)

Where x is the input raster and y is the exponent I am attempting to raise the raster to. I am getting the error of invalid expression, any guidance would be very helpful!


texts: Is it possible to get a dynamic legend in QGIS 3.6?

texts:
I need to send a GetPrint request for a map with a dynamic legend that only shows layers/items/symbols relevant for the current map frame.

In QGIS 2.14 one can select "Filter Legend by Map Content" in the composer and export a filtered legend, but if you use a GetPrint request no legend is shown. I have managed to work my way around this by adding together GetLegendGraphics requests for each layer.

Now I'm trying to upgrade from QGIS 2.14 to QGIS 3.6, but I'm not sure this is possible. In QGIS 3.6.3 "Filter Legend by Map Content" does not seem to work in the composer and GetLegendGraphic requests are apparently ignoring the BBOX parameter.

Is this a bug? Is there any way to get around this, or do I need to keep on using QGIS 2.14?


texts: Is Lock Graphics functionality from ArcGIS 9 ArcScript available for ArcGIS 10.1?

texts:
I have created grids and have some layout graphics and I dont want to accidentally move these elements.

I found a old script here which does what I want, and I understand that the Production Mapping extension does this but is more advanced.

I cannot recompile from ArcGIS 9 to 10 (not smart enough), but do you know if there is an finished DLL file like the one in the link but for 10.1?

I also want to know if it is possible to lock Map frame element and/or save layout settings in ArcMap 10.1


texts: Getting OpenLayers to scale past nativeZoom like Leaflet?

texts:
Is there any way to get OpenLayers 3 to scale past the native zoom as is possible in Leaflet?

For reference, in Leaflet simply set maxZoom higher than maxNativeZoom and it stretches or perhaps scales the tile for a closer view.

In this example tiles go to 16. I'd like to scale 16 up to 18 for a closer view:

map = new ol.Map({
    target: 'map_div',
    renderer: 'canvas',
    view: new ol.View({
        center: ol.proj.transform([X_LO,X_LA], 'EPSG:4326', 'EPSG:3857'),
        zoom: 16,
        minZoom: 4,
        maxZoom: 16
        })
    });

L1=new ol.layer.Tile({
    source: new ol.source.XYZ({
        url:'http://blahblah.com/map_tiles.php?Z={z}&X={x}&Y={y}'
        }),
    maxResolution:152.87,// 9
    minResolution:1.94// 16
    });
map.addLayer(L1);

Solution from @erilem

map = new ol.Map({
    target: 'map_div',
    renderer: 'canvas',
    view: new ol.View({
        center: ol.proj.transform([X_LO,X_LA], 'EPSG:4326', 'EPSG:3857'),
        zoom: 16,
        minZoom: 4,
        maxZoom: 18
        })
    });

L1=new ol.layer.Tile({
    source: new ol.source.XYZ({
        url:'http://blahblah.com/map_tiles.php?Z={z}&X={x}&Y={y}',
        maxZoom: 16
        })
    });
map.addLayer(L1);


texts: Is GeoServer getting coordinate output in GeoJSON WFS wrong?

texts:
When I do a WFS request to GeoServer and specify the SRS like this:

http://myserver.com/geoserver/wfs?
    service=WFS&
    version=2.0.0&
    request=GetFeature&
    outputFormat=application/json&
    srsName=EPSG:4326&
    typeName=MYWORKSPACE:mylayer&

I get exactly what I expect and it works fine.

When I remove the srsName=EPSG:4326 parameter I get GeoJSON that still identifies itself with:

"crs": {
    "type": "EPSG",
    "properties": {
        "code": "4326"
    }
}

but the order of the coordinates is reversed, so it doesn't work.

Are both of these valid EPSG:4326 output, or is one of them wrong?

The layer I'm working with has these settings:

Native SRS: EPSG:4283
Declared SRS: EPSG:4326
SRS handling: Reproject native to declared


texts: Convert point layer to text in v8 design file

texts:
I have a shape file point layer with field  "name" having attribute values e.g. 'city name'. I want these attribute values as text in v8 design file using FME.
Any help please.


texts: How to split filenames for their first word in Python?

texts:
I have 2 shapefiles, each one placed in a different folder:

"C:\...\Folder O\Germany (O).shp"
"C:\...\Folder P\Germany Berlin (P).shp"

The following script finds the shapefiles and contains a couple of if conditions. The first if condition strips the "(O).shp" and "(P).shp" from the shapefiles. This works perfectly.

I want to set up a second if condition whereby the first name of both shapefiles matches, so in this case it would be "Germany". I thought the string.split method would be ideal but I'm clearly not implementing it correctly.

Could someone kindly point out my mistake please?

def run():
    for o in glob.glob("C:\Users\me\Desktop\Test\Folder O" + "\*.shp"):
        for p in glob.glob("C:\Users\me\Desktop\Test\Folder P" + "\*.shp"):
            if os.path.basename(o).rstrip("(O).shp") == os.path.basename(p).rstrip("(P).shp"):
                print p
                # works perfectly up to here

                if os.path.basename(o).split(o[[1]]) == os.path.basename(p).split(p[1]]):
                    print p


texts:
First attempt, I want to create a table at CartoDB from another Postgis database:

snow:Docker-Postgis alasarr$ ogr2ogr --config CARTODB_API_KEY daf3960b41733ef71e03ea77019642761d547f5b -f CartoDB "CartoDB:alasarr" PG:"host=192.168.59.103 user=postgres dbname=eiel_huesca port=5433" "geometries_eiel.alumbrado_ok"

// output
CartoDB driver does not support data source creation.

Second attempt, I try that ogr2ogr create a table at CartoDB from a Shapefile:

snow:Docker-Postgis alasarr$ ogr2ogr --config CARTODB_API_KEY XXXX  -f CartoDB "CartoDB:alasarr" ~/dev/naturalearth/ne_110m_admin_0_countries/ne_110m_admin_0_countries.shp

// output:
CartoDB driver does not support data source creation.

Third attempt:

I import the shapefile at CartoDB.
Truncate the table. Truncate table ne_110m_admin_0_countries;
I try to import using ogr2ogr.

Here the command:

ogr2ogr --config CARTODB_API_KEY XXXX  -f CartoDB "CartoDB:alasarr" -append -t_srs EPSG:4326  ~/dev/naturalearth/ne_110m_admin_0_countries/ne_110m_admin_0_countries.shp

//output
ERROR 1: HTTP error code : 400
ERROR 1: Error returned by server : relation "ne_110m_admin_0_countries" already exists
ERROR 1: Terminating translation prematurely after failed
translation of layer ne_110m_admin_0_countries (use -skipfailures to skip errors)

I'm using Gdal 1.11.0 that in theory supports CartoDB format.

Any ideas? Should I try with a newer version of GDAL.


texts:
I have a problem with my first QGIS 3 plugin. I want to calculate lenght each segment in my polyline and next creating buffer equal 1/2 length of line.

The first 13 lines code have to do calculating the length of the selected layer and it works. But I have a problem with the buffer. I tried many ways that I found on the internet and it still doesn't work. Now I have no errors, but the plugin does not return anything.

Can you solve this problem?

if result:
    for layer in QgsProject.instance().mapLayers().values():
        features = layer.getFeatures()
        coluna = QgsField('Length', QVariant.Double)
        res = layer.dataProvider().addAttributes([coluna])
        leng = layer.dataProvider().fieldNameIndex
        layer.updateFields()
        fieldIndex = layer.dataProvider().fieldNameIndex( "Length" )
        layer.startEditing()
        for f in features:
            geom = f.geometry()
            leng = geom.length()
            layer.changeAttributeValue(f.id(),fieldIndex, leng)
        result = processing.run("native:buffer",
            {'INPUT': layer,
            'DISTANCE' : 100,
            'SEGMENTS' : 10,
            'DISSOLVE': True,
            'END_CAP_STYLE' : 1,
            'JOIN_STYLE': 0,
            'MITER_LIMIT': 10,
            'OUTPUT': 'memory:buffer'})
        layer.commitChanges()
        layer.updateFields()


texts: Cannot set the natural color for Sentinel-2 raster data

texts:
I'm trying to create a new raster data with nature color (B4-B3-B2) for the left raster data in QGIS 3, however, I failed. This should be similar to the right one.

Both of them were downloaded from https://scihub.copernicus.eu/dhus/#/home.

The left raster data:

  S2B_MSIL2A_20200227T030659_N0214_R075_T48PXS_20200227T071451

The right raster data:

  S2B_MSIL2A_20200227T030659_N0214_R075_T48PYS_20200227T071451

How can I solve this problem?

Edit:
Thank Jake for helping me.
Actually I just followed this step.
First I imported B4, B3, B2 raster bands

Then I build a new virtual raster

Next I assigned the bands to the channel. I don't know why there is no band 4 in even thought I've added before.

However it's still not working

Edit 2: It's worked!
However I don't know why the band order is 3-2-4 for R-G-B instead of 4-3-2. Could anybody explain it ? (Sorry that I'm a new member)


texts: Determining which shapefile lies in another using ArcGIS Desktop?

texts:
I am trying to find out which region (shapefile 2) each polygon in shapefile 1 is in.  In some cases it may be in more than one region.  In this case I am trying to determine in which ecosystem (Terrestrial Ecosystems of the World) each of the World's protected areas are in.

I have tried most combination of things in ArcMap 10.3 without success (I can get some to work but not others).


texts: Offset shapefile to start from prime meridian

texts:
I have Natural Earth shapefile (for example this file: link) covering whole globe in -180:180 range.
I want to offset this file, so that it starts from Greenwich prime meridian instead -180.

How can I do this with gdal/ogr?


texts: Is it possible to add custom SQL-based layers in GeoServer?

texts:
I would like to add a layer to GeoServer that builds its features from an SQL query executed against a PostGIS database (something like a PostGIS layer in QGIS, but customized with additional SQL query).

Is it possible to have not a whole PostGIS geometries table, but some 'SELECT' from it to be published by GeoServer?


texts: Has anyone found that using CalculateDefaultGridIndex in arcGIS9.3 does not return the ';' delimited string suggested in the documentation?


texts: How to create a shapefile having Voronoi diagrams from the shapefile of an existing region?

texts:
I have a shapefile of a region and I want to generate random sub regions of that shapefile using Voronoi class of scipy.spatial package of Python. But I am a newbie to GIS. I found out about Voronoi from StackOverflow comments. I don't know how to do that for a city with contoured boundaries.

Any suggestions, code snippets or suggestions will be a great help to me.


texts: Setting placement to a group of selected labels

texts:
At work with the Ontario Geological Survey where we use ArcMap (10.6), I can select a number of annotation text and change the properties of the selected text.

To date, at home using QGIS 3.18 (Linux), I am unable to figure out if this is even possible, let alone how to do it.

Sorry for not being clear enough, here's some additional info.

Easting and Northing labels at the end of each line (nodes) whereby the Easting labels at the top of the map are centered on the node 1mm above the neatline, the Easting labels are centered 1mm below the neatline, the Northing labels on the right side of the map are 1mm right of the neatline and the Northing labels on the left side of the map are 1mm left of the neatline.

In ArcMap, I can select the left Northings and apply right alignment, select the right Northings and apply left alignment, etc. It appears that alignment works differently in QGIS but I've figured out that Placement>Offset from point using the top/bottom/right/left quadrants with a 1mm offset does exactly what I want it to do but I don't know how to select the labels on the left and apply their specific Placement settings that are different from each of the groups of labels.


texts:
ArcGIS Pro 2.4.0 - all extensions

I have a dataset showing Water Mains (lines) for an area. I have been tasked with finding water mains that go from a given diameter to a smaller diameter, and back to a larger diameter.

For example: An 8" water main that links into a 6" water main, which then links back into an 8".

Ideally I would want to select just the 6" water main in the example but I'm unsure how to go about this. Another problem I think I might run into is that the water main lines are made up of many small segments. For example a 40 foot line is actually 8 different 5" segments.

I'm unsure whether I can use a tool for this, definition query, or write a Python script which I'm not very good with.


texts: How to create a NoData raster mask in ArcGIS

texts:
I need to create a mask of NoData pixels from an elevation raster dataset. ArcGIS directly recognizes the NoData values as such, so I don't see an option to address them e.g. in the Raster Calculator to create a simple raster Mask where 1 stands for the data gaps in my original elevation raster.

Could anyone help me with a default method or at least a flexible workaround for this?


texts: Huge file size after averaging two rasters

texts:
After I average two raster files (singleband gray, initial size ~700MB each) through QgsRasterCalculator() with formula "Raster1@1/2 + Raster2@1/2" the size of the resulting file is 2.9GB.

Is there any possibility to keep the filesize somehow lower by modyfing something? Switching to some external toolbox?


texts: qgis2web Polygon Layer Overwhelming Popup Selection of Point Layer

texts:
I have a map of a point shapefile layer over a poylgon shapefile layer of legislative districts. When I export with qgis2web, the popup will only show the fields from the polygon layer. Here are examples of what happens with OL3:

With Leaflet, I have to change the polygon background to transparent or it covers over the dot layer. The issue is still the same with the popups not ever displaying the point layer information when the polygon layer is engaged.

Is there anything I can do to allow for the point layer to be clicked to a popup in this situation?

If it is relevant I am using Windows 10, QGIS 2.12.3 and qgis2web 1.19.0


texts: 'Add XY Data' and 'Make XY Event Layer' are plotting points at an inaccurate latitude

texts: I have a .csv file with 10 wells and their respective longitude (X) and latitude (Y). Every time I try 'Add XY Data' or 'Make XY Event Layer', my wells are being plotted at some random place in the world and not on our job site. So to see what was going on, I added two new columns to the attribute table and used field calculator to calculate the longitude and latitude of wherever these wells are plotted. They are plotting at the correct longitude, but the latitude for each well is plotting at 0.000369. I can't understand what is going on. I have done these exact steps hundreds of times and have never had this issue.


texts: Connect not working for custom QGIS form

texts:
I have a custom form for editting feature attributes. However, I can't seem to correctly disconnect and reconnect certain bindings. The init code is below and should work for any simple QT widget with buttons.

It doesn't throw an error. The OK button is disconnected from the original connect that QGIS assigns. But it does not reconnect to the validate function.

Can anyone point me in the direction of the mistake here?

from PyQt4.QtCore import *
from PyQt4.QtGui import *

# QGIS passes in the dialog, he layer_id and the feature_id
def form_open(myDialog,layer_id,feature_id):
    # QGIS passes in the dialog, he layer_id and the feature_id
    # I send all this to be handled by its own class. So I can use
    #'self' and avoid having to use globals...
    return MyClass(myDialog,layer_id,feature_id)

class MyClass(QDialog):
    def __init__(self,myDialog,layer_id,feature_id):
        self.ui = myDialog
        self.buttonBox = self.ui.findChild(QDialogButtonBox,"buttonBox")

        # The problem is here...
        self.buttonBox.accepted.disconnect(self.ui.accept)
        self.buttonBox.accepted.connect(self.validate)

    def validate(self):
        msgBox = QMessageBox()
        msgBox.setText("I validate here...")
        msgBox.exec_()


texts:
I've been hunting around ArcGIS 9.3 and having no success at finding a tool/function that can give me the median (not the mean/avg) tabular value of a spatial database.  For example, if I have a data layer of mountain peaks, I would like to be able to determine the elevation of #35 of 70 mountain peaks that exist within a county.

Conceptually it should be pretty easy.  Count the total records, sort by elevation ascending, get the value from the elevation field for record # (count/2).  Would I need to script this out in Python, or am I missing something obvious that would allow me get this more quickly?  I find it hard to believe no one else has had a need for this function and created a solution.

To add some complexity to the matter, what if I wanted the median elev value for peaks in every county within the state in a single table?  Basically the same output format as Frequency or Summary Statistics tools when using the Case Field option.

I was able to adapt the script linked by @VietThanh Le to meet my needs. I have posted the updated script and ArcToolbox at: http://resources.arcgis.com/gallery/file/geoprocessing/details?entryID=A0CD1751-1422-2418-882E-001EE0DC0D35


texts: Neighbourhood statistics in Google Earth Engine after reclassifying (remap)

texts:
I want to reclassify a categorical raster in Google Earth Engine, and then calculate the proportion of non-0 pixels in a neighbourhood. I've remapped to values of 1 or 0.

I'm having trouble with the reduceNeighbourhood() function. With the code below, I expected this to sum the values of all the pixels in the neighbourhood (i.e., sum the 1s), but it doesn't seem to do this? My plan was then to divide this by the number of pixels in the neighbourhood using either a constant (i.e., calculated from my ee.Kernel.square() radius) or from ee.Reducer.count().

// import image (image with multiple bands)
var image = ee.Image('projects/mapbiomas-workspace/public/collection3_1/mapbiomas_collection31_integration_v1');

// xMin, yMin, xMax, yMax.
var testArea = ee.Geometry.Rectangle([-60,-20,-40,0])

// Get information about the bands as a list.
var bandNames = image.bandNames();
print('Band names: ', bandNames);

// select single band
var v1985 = image
  .select("classification_1985")
  .clip(testArea);

print(v1985)

// Map.addLayer(v1985)

// do what I want with one image as a trial.
var rc_1985 = v1985.remap([2,3], [1,1],0);

// Map.addLayer(rc_1985)

var k1_1985 = rc_1985.reduceNeighborhood({
  reducer: ee.Reducer.sum(),
  kernel: ee.Kernel.square(15, 'pixels')
  });

print(k1_1985)
Map.addLayer(k1_1985)


texts:
My company recently installed MapMarker v28.0, and the product description says it comes with something called MapMarker Streets.

My company has an older web application, and the streets we use are kept in a file called streets.tab. It's outdated, and a reason why we updated to the latest version of MapMarker was to get the updated street information.

I can't find any information on how to use MapMaker Streets or access it as data.


texts: High PostgreSQL/PostGIS memory and disk usage when saving QGIS project

texts:
Why would simply pressing the "save" button in QGIS cause my local PostgreSQL instance to cause memory and disk usage to spike to "Not Responding" levels?  Seems to me like the PostgreSQL shouldn't even care about saving the project.  Are there any suggestions for improving the PostgreSQL DB performance?
UPDATE: Through some debugging in QGIS, I think I at least now understand why I am getting so much throughput to the DB on save.  QGIS is throwing dozens of errors like below (one for just about every layer in the project) in the PostGIS tab of the Log Messages Panel:
2020-09-02T09:19:51     WARNING    Error executing query:
SELECT id,styleName,description FROM layer_styles
WHERE f_table_catalog='<dbname>'
AND f_table_schema='<shema_name>'
AND f_table_name='<layer_name>'
AND f_geometry_column='geom' ORDER BY useasdefault DESC, update_time DESC
2020-09-02T09:19:51     WARNING    No styles available on DB

I am not intentionally saving any layer symbology in the DB.  Any thoughts on what is causing this query to happen, and how I can remove it?
Background:
When I save my existing QGIS project (qgz) that has about 100 layers, all sourced from PostgreSQL/PostGIS, I frequently get a "Not Responding" message from QGIS.  While trying to debug the problem, I started watching the system resources and noticed that during these not responding times, the memory and disk i/o go crazy for the PostgreSQL instance, spiking the disk i/o to 100% (as high as 350 MB/s) and consuming 20% of the CPU.
I am running QGIS 3.12, PostgreSQL 9.5, and PostGIS 2.3 on a Windows 10 laptop with an Intel i5 CPU with 2 cores/4 processors and 12GB of RAM (physical memory) and a 1TB SSD hard drive.


texts: Why unable to create interactive FeatureSet features for ModelBuilder tool?

texts:
I am attempting to set up an "interactive clip" tool that will allow the user to draw clip features on the fly and use those as the clip features in the Clip tool in the analysis toolbox. I have added the Clip tool to ModelBuilder and exposed the clip features as a model parameter. I also changed the type to a FeatureSet, which I understand is supposed to allow the user to interactively enter the features when the tool is opened within ArcMap.

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//002w0000004t000000

However, when I open the tool in ArcMap, the Clip Features input does not allow the interactive drawing of features, but rather only allows the user to select an already existing feature class. ArcGIS help and the above link indicate that under this setup the FeatureSet should by default point to an in_memory location where the user can temporarily draw features.

Here's how I have it set up in ModelBuilder:

Here's the resulting tool interface:

What am I missing here? Any ideas why this is not working as expected?


texts: Changing editing capabilities on MapInfo *.tab file using QGIS?

texts:
When I load a .Tab file into QGIS I am unable to make changes to it.
The information on the MetaData Tab reads the following:

"Editing Capabilities of this layer: Add Attributes, Fast Access to Features at ID"

How can I change the Editing Capabilities so that I can edit the features with this .TAB file.


texts: Clustering of Spatial Data in R or Python

texts:
I am new to GIS.

Problem Domain :-
I have a stoppage data of different vehicle i.e. the point where vehicle stopped.There are some specific vehicle which stops at some particular point. So, I am trying to find out such specific point(anomaly detection) in the midst of random stoppages.
I am assuming that Cluster analysis would let me know about the Percentage of vehicle at that stoppage point. If only a particular vehicle stops at that point , then obviously that's an anomaly point.

I have a large set of (latitude and longitude)spatial data. I am confused on which clustering method to adopt. I have came across two density based approach: DBSCAN and OPTICS. I am doubtful about the two approach since I don't have a particular minPts( 1 in my case). Obviously, I can't use K-means approach , K is unknown.

I have come across one more approach i.e X-means.


texts: Constructing boundary polygon around raster image when mask layer is ill-defined using GDAL with Python?

texts:
I have a raster image and I want to construct a new vector layer that has the boundary of the raster, using python. Below is a sample of the image.

The problem is that the nodata value is set to "None" and the white pixels that constitute empty space are coded as 255. So that means the mask is not going to mask out the whitespace as areas of nodata.

Just to confirm I checked rasterio I found the following:

import rasterio
src = rasterio.open('../image.tif', mode='r')
data_mask = src.read_masks(1)
data = src.read(1)

data looks like showing that the upper left region of the image is white and empty with the 255 value.

array([[255, 255, 255, ..., 254, 255, 255],
       [255, 255, 255, ..., 254, 255, 255],
       [255, 255, 255, ..., 254, 255, 255],
       ...,
       [255, 255, 255, ..., 248, 252, 248],
       [254, 254, 254, ..., 255, 235, 253],
       [255, 255, 255, ..., 255, 255, 245]], dtype=uint8)

And the data_mask looks like the white in the upper left of the image is valid. These values should read as 0 but instead they are 255--meaning they are valid.:

array([[255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255],
       ...,
       [255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255]], dtype=uint8)

So to extract the boundary polygon I tried to use the code from another SE post.

gdal_translate -b mask -of vrt -a_nodata 0 test.tif test.vrt
gdal_translate -b 1 -of vrt -a_nodata 0 test.vrt test2.vrt
gdal_polygonize.py -q  -8 test2.vrt -b 1 -f "ESRI Shapefile" testdata.shp

as well as the alternative code:

gdalwarp -dstnodata 0 -dstalpha -of GTiff foo1 foo2
gdal_polygonize.py foo2 -b 2 -f "ESRI Shapefile" foo3

But both of these approaches produce a boundary polygon that reflects the mask, which includes the whitespace.

So I need to know how to adjust this code to first mark the pixels that are 255,255,255 as nodata and then I should be able to run the code above right. The other issue is that the resulting polygon will probably have little holes in it corresponding to areas of white in the image.

Is there a good way to filter those out when creating the polygon?


texts: Can points not in a feature be roped into an existing feature in Trimble Pathfinder?

texts:
I am editing a Trimble Pathfinder .cor file to reduce standard deviation on collected points in features and meet government standards.  Some features have a large standard deviation between contributing points, but not many points collected as part of the feature.  However, there are numerous points right near the feature which are not part of any feature, but it's pretty obvious that they were collected in the same place as the feature (i.e. they're all <3m away from the feature).  The feature's points spread out across a 10m area and average to a point right near these unassigned points, but not many points were collected as part of the feature so if I delete the outliers to tighten the cluster I don't meet government standards for # of points collected.

I would like to tell Pathfinder to take the unassigned points and add them to the existing feature, thereby increasing the number of nearby points in it to average, and allowing me to delete the points furthest away to create a tighter cluster (and smaller standard deviation) on that feature.

Anyone know if this is possible?


texts:
How do you generate a species distribution probability map from GAM models? Or in other words, how do you set up a GAM equation in GIS to produce probability of species distribution?

It is simple to apply generalized linear models (GLM) or equation in GIS due to the independent variables have a liner relationship with dependent variable but for GAM, I am not able to do it because the parameters in GAM are non-linear predictors and also their estimates (β) are not be calculated.


texts:
I want to show different types of surveyed land parcels which are available in  shapefile format, in our country map.  That map will be published in a web site and it will be available for surveyors.

I am using QGIS version 2.18.9 and PostgreSQL database. Can I do this using QGIS or should I use QGIS server version?


texts: Getting centroid of multipart polygon using ArcPy

texts:
I am trying to get the center of a multipolygon, I have tried doing:
with arcpy.da.SearchCursor(fc, ['SHAPE@']) as cursor:
    for row in cursor:
        multipoly = row[0]
        centroid = multipoly.centroid

The feature class is using WGS84 coordinates.
This only seems to grab the centroid of only one of the parts. Ideally, I would like to get the center point of all the parts in the multipolygon.
Any suggestions for approaches I can pursue?


texts: Topographical Elevations for Central Africa and Cameroon

texts: I am using QGIS and I have been faced with the problem of topographical maps. From searching online I could only find data for the United States. Google maps offers elevation contours, however I am looking for a solution for Central Africa, Cameroon that is more specific. I am also looking for any other resources that can equally show the different geographic basins in the target area.


texts: Unable to use GPS coordinates for GEOTiff map

texts:
I have been using a map to calculate some distances at sea. The map had a bit low resolution so I got another one with better resolution (both GEOTiff), but with the amounts of distances I need to calculate the algorithm quickly grew to require a run time of 6 months to complete. I want to try to make cropped versions of the map with the point I am calculating distances for in the center, and a box with sides of 200km around it. I found some example code using rasterio to do this:

MIN_LAT = 57.75
MAX_LAT = 71.35
MIN_LON = 6
MAX_LON = 31.45
r = 200
c = (22.3542, 70.14145)

dlon, dlat = 150*0.015060, 100*0.008983

def check_val(val, latorlon):
    if latorlon == "lat":
        if val < MIN_LAT:
            return MIN_LAT
        elif val > MAX_LAT:
            return MAX_LAT
    elif latorlon == "lon":
        if val < MIN_LON:
            return MIN_LON
        elif val > MAX_LON:
            return MAX_LON

    return val

points = [(check_val(c[0] - dlon, "lon"), check_val(c[1] - dlat, "lat")),
           (check_val(c[0] + dlon, "lon"), check_val(c[1] - dlat, "lat")),
           (check_val(c[0] + dlon, "lon"), check_val(c[1] + dlat, "lat")),
           (check_val(c[0] - dlon, "lon"), check_val(c[1] + dlat, "lat"))]

geoms = [{'type': 'Polygon', 'coordinates': [points]}]

with rasterio.open('map_100x100_original.tif') as src:
    out_image, out_transform = mask(src, geoms, crop=True)
out_meta = src.meta.copy()

out_meta.update(
    {
        "driver": "GTiff",
        "height": out_image.shape[1],
        "width": out_image.shape[2],
        "transform": out_transform
    }
)

with rasterio.open("map_100x100_cropped.tif", "w", **out_meta) as dest:
    dest.write(out_image)

The code works as intended with the old map, however, with the new map I get an error message from rasterio saying: "WindowError: windows do not intersect" followed by "ValueError: Input shapes do not overlap raster".

Using some code I found here I got the following information on the two maps:

old map

None
WGS 84

new map

WGS 84 / UTM zone 33N
WGS 84

Using some other code I found here I got the corners of the two maps

old map

(-180.0, 90.00000000000001) (180.00000000007202, -90.000000000036)

new map

(1121948.79, 6426051.97) (-99551.21, 7962751.97)

Using QGIS I got some coordinates by hovering on the map which I could plug in to the code and it would work. For example, the GPS coordinates below

(59.48129, 5.89213)

Would be something like

(-60045, 6609466)

But I need to be able to use GPS coordinates, so is there a way I could fix the new map?


texts: st_startpoint(the_geom) returning empty values

texts:
When I upload a shapefile to PostgreSQL using the QGIS tool I can create the routing topology following this tutorial without any problem.

But when I use the PostGis Shapefile Import/Export Manager or the SQL generated by the shp2pgsql.exe I can't get any values for the x1,x2,y1 and y2 fields.

UPDATE edge_table SET x1 = st_x(st_startpoint(the_geom)),
  y1 = st_y(st_startpoint(the_geom)),
  x2 = st_x(st_endpoint(the_geom)),
  y2 = st_y(st_endpoint(the_geom)),
cost_len  = st_length_spheroid(the_geom, 'SPHEROID["WGS84",6378137,298.25728]'),
rcost_len = st_length_spheroid(the_geom, 'SPHEROID["WGS84",6378137,298.25728]');

Only cost_len and rcost_len have values after running that script.

And because of that I can't get any route at all. What am I missing?


texts: Correlation/relationship between map layers in R?

texts:
I have a point shapefile concerning asthma instances, a point shapefile concerning rainfall, and a polygon shapefile with population information i.e gender, ages etc.

I would like to research if there is a spatial relationship between the location of asthma and rainfall, and also if there is a spatial relationship between instances of asthma and  the population data. Ultimately I would to be able to say:

  There are/are not higher instances of asthma in areas that receive more rainfall and/are there are/are not higher instances of asthma in areas with more men/women.

I have done research on the areas of spatial correlation and spatial autocorrelation which I assumed would be related to this area. However, it seems these only correspond to data within one map layer i.e. you cannot technically look for spatial correlation between two map layers.

Is there any way in R to assess correlation/a relationship between two point layers, or a point layer and a polygon layer?


texts: Adding M values to a shapefile so the distance along the line is displayed as a derived attribute in QGIS

texts:
I'm new to QGIS which I'm using to create data for a road network. I have a road network setup as a polyline shapefile, but need to identify the chainage long the network sections to use in other datasets. I have loaded a PolylineM shapefile and in that the distance along the line is given when I use the identify features tool as the derived attribute Interpolated M. I have tried using the 'Set M Value' tool to create a new layer which stores the lengths as the M value, but can't get that to work so that it reports the distance along the line the same way the PolylineM shapefile does.
Does anyone know if this method would/should work, or how I can go about setting this up?


texts: Determine if point is within an irregular polygon using Python

texts:
I would like to create a mask of a grid where points are evaluated as either 1's or 0's depending on if they lay outside of an irregular (i.e. not rectangular) polygon.

I have a set of high-latitude coordinates in lat/lon that define the boundary of an area I would like to make a polygon.

From this polygon I would then like to evaluate another set of lat/lon points as to whether they lay within the polygon.

Is there a library best suited for creating such a polygon (given a list of lats/lons that define it)? shapely appears to be geared towards 2D/projected points.

What would be the best method of evaluating other points as to whether they lie within or out of this polygon?

Thanks!


texts: Flattening multiple fields using Pivot Table in ArcPy?

texts:
I'm trying to flatten a table, but I want to keep/relate certain fields tied to each other. For context, I ran a Generate Near Table to try to find the traffic volumes near intersections.

I am trying to keep the "Route_Name" field tied to the "Curr_Vol" field, while flattening it to the input field "IN_FID". (IN_FID is intersections, so we need to know what traffic volumes are near the intersection, and what route name is related to each traffic volume)

I've tried doing a pivot table, but I can't flatten the table to IN_FID completely while keeping the Route_Name and Curr_Vol fields.

Original data:

After trying to run Pivot Table:

What data needs to look like in the end:

The result I want is a table where each row will have; The unique IN_FID value, the Route_Name & Curr_Vol values from every row in the original table, arranged as name1, vol1, name2, vol2, and so on.

I'm wondering if an ArcPy script would be the best solution?

However, I don't know how to write that complex of code.


texts: Set proximity requirement for ESRI Collector observations

texts: I am trying to create a Collector map that allows field crews to update a "status" attribute for features in a feature class.  What I would like to do is limit the features that can be edited by the field crew to the features within 100 yards of the user's current device location based on GPS.  Is that something that is supported by the Collector app?


texts: How to find the mutation area and time according to Mann-Kendall mutation test. When using mutation test SCF trend, we can get the SCF trend of the region according to the trend chart, but how to get the mutation time of a specific area?


texts: Remove black borders when converting GeoTiff to MBTile with gdal_translate

texts:
I am using gdal_translate to convert a fairly large geoTiff to MBTiles.

As the geoTiff is reprojected, it is not square anymore in Google Mercator and the parts that are not covered appear black when I would like them to be transparent.

I think it could have something to do with the -a_nodata parameter and I have tried with "-a_nodata 0", but it did not change things.

I have tried with both PNG8 and PNG. It did not change things either.

Here is what my command looks like:

gdal_translate -a_nodata 0 -of MBTILES -co TILE_FORMAT=PNG8 France.tif france.mbtiles

Any idea what I am doing wrong?


texts: Combining two search cursors (one for table and other for feature class)?

texts:
I have a problem from combining two search cursors (one for a table and other for a feature class), there are missed records in the first table, that have not been considered by the first search cursor, I have tried many trials for error handling like join both tables and using try and except syntax,  but all of that have not worked, also after reading answers on previous questions here, I have tried (cursor.reset() function), but it gives no end for the script.
 here you are the part of my code:

cursor3= arcpy.da.SearchCursor(table1, [field1, field2])
for row in cursor3:
    try:
        cursor4= arcpy.da.SearchCursor(FeatureClass1, ['OBJECTID'])
        for row1 in cursor4:
            try:
                if row1[0] == row[0] or row1[0] == row[1]:
                    FeatureClass2= arcpy.env.workspace+"\\FeatureClass2"
                    where = "\"%s\" = %i OR \"%s\" = %i" % ('OBJECTID', row[0],'OBJECTID',row[1])
                    arcpy.Select_analysis(FeatureClass1, FeatureClass2, where)

                    # Local variables:
                    FeatureClass2 = "FeatureClass2"
                    FeatureClass2_Intersect = "FeatureClass2_Intersect"
                    # Process: Intersect
                    arcpy.Intersect_analysis(FeatureClass2, FeatureClass2_Intersect, "ALL", "", "LINE")

                    arcpy.AddField_management(FeatureClass2_Intersect,'Value',"DOUBLE")
                    cursor2= arcpy.da.UpdateCursor(FeatureClass2_Intersect, 'Value')
                    for row2 in cursor2:
                        try:
                            if row2[0] ==6:
                                row2= cursor2.next()
                                if row2[0] ==7:
                                    row2[1]=float(42)
                                    cursor2.updateRow(row2)# give a value to value field
                        except:
                            del cursor2

            except:
                row1= cursor4.next()

    except:
        row= cursor3.next()
        del cursor3


texts: Importing CSV (via VRT) to PostGIS is slow with -skipfailures?

texts:
I am importing about 60 million records of csv data (about 25gb) to Postgresql server using ogr2ogr

After googling million times, I found out that I can make VRT file and uploading my csv file using ogr2ogr.
Here's the code I used to importing data and vrt file

<OGRVRTDataSource>
    <OGRVRTLayer name="combine">
        <SrcDataSource>D:\combine.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <LayerSRS>WGS84</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="longitude" y="latitude"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

ogr2ogr -overwrite -progress -gt 999999999 -skipfailures -f "PostgreSQL" PG:"host=000.000.000.000 port=0000 dbname=myDB user=me password=youtellme" -lco GEOMETRY_NAME=shape -t_srs EPSG:3857 D:\input.vrt -nln output

I thought I could adjust -gt option to speed up the importing time, yet it doesn't help much.

It has been 3 days since I started to importing data to Postgresql server and only 20 million data are uploaded so far. (means that It will take almost a week to upload)

Seems like my ogr2ogr code (or VRT) can only import 100 - 200 records per second depending on the server status (see the screen capture below)


texts: QGIS - Import CSV column names as numbers

texts:
I want to import a CSV table (no geometry) to QGIS. The column names which are numbers (e.g. 2004) are not getting imported.

This is my test CSV file:

Name;Eigenschaft;2004;test
derName;1;21;21
sdfsfaf;11;11;15

This is how it looks in QGIS:

The third column is not called 2004 like I want but field_3

How to solve this?


texts: can the geometry element be accessed within CSVLayer in ArcGIS for JS 4.x

texts: I was wondering if there is a graphic or geometry element that can be access within CSVLayer in arcgis JS 4.7, I have been unable to extract the geometry which I want to save as a JSON or array that I can add to my featurelayer using applyEdits()


texts: How to add read-only SQL Textbox calculated from a SQL statement in Qgis Form field

texts:
I have two tables which have a relation one to many, because the second table stores the editing history of the first table which has geometry, and they are connected through an ID.

I want to create a read-only Sql Textbox in the first table Qgis forms which value is calculated through a sql statement.

For example if i click on a geometry the field search on the second table for the rows that has the current ID and display the desire data from a certain column.

Like in this screenshot:

Where at the end of the form you can see the displayed records from the other table through a sql statement.


texts:
In my database table there are four fields(lat,lon,sensid,value).i am displaying a vector layer(point feature) by retrieving the lat, lon values from database.
Now i want to display the value, whenever i click on the feature based on id.
So how to get the feature id(database sensid) on click.I am getting lat lon values on click, but i need to get sensid.

my code is

var  str_array =["8541652,3397019,TATM1", "8786758,3536158,TATM2"]
 str_array = jsstr.split(';');
         var a,b,c;
    for (i = 0; i < str_array.length; i++) {
   var text = str_array[i] ;
    var  str_coord = text.split(',');
     a = parseFloat(str_coord[0]);
     b = parseFloat(str_coord[1]);
  c = parseFloat(str_coord[2]);
var feature=new OpenLayers.Feature.Vector(new OpenLayers.Geometry.Point(a,b));

           vectorLayer.addFeatures(feature);
     mapPanel.map.addLayer(vectorLayer);
var selectCtrl = new OpenLayers.Control.SelectFeature(vectorLayer);

    function createPopup(feature) {

        popup1 = new GeoExt.Popup({
            title: 'My Popup',
            location: feature,
            width:200,
            html: '
'+feature.geometry.x +','+feature.geometry.y+'
',
            maximizable: true,
            collapsible: true,
            anchorPosition: "auto"
        });
       // unselect feature when the popup
       // is closed
        popup1.on({
            close: function() {
                if(OpenLayers.Util.indexOf(vectorLayer.selectedFeatures,
                                           this.feature) > -1) {
                    selectCtrl.unselect(this.feature);
                }
            }
        });
        popup1.show();
    }

   //create popup on "featureselected"

    mapPanel.map.addControl(selectCtrl);
   selectCtrl.activate();
}
vectorLayer.events.on({
        featureselected: function(e) {

            createPopup(e.feature);

        }
    });


texts: How to insert non-spatial features into a Polygon shapefile - FME 2015

texts:
Here is the issue:

I have the original data (a shapefile with spatial a non-spatial features) which I need to transform: (highlighted those non-spatail features)

I need to make some modifications on the data (and for that I'm using FME Desktop 2015.1).
Once performed, I need to write again to a shapefile, trying to keep also those non-spatial features... but NO luck with this structure...

Here we can see those highlighted with a feature type: Null. So, once I try to write into a Shapefile, it is just written those spatial ones. But, what I need is also those non-spatial.

------- Edited to add log: -----------------

08 08:21:36|  37.6|  0.0|WARN  |Polygon feature must have at least 4 coordinates...rejecting
08 08:21:36|  37.6|  0.0|WARN  |REJECTING BELOW FEATURE:
08 08:21:36|  37.6|  0.0|INFORM|+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
08 08:21:36|  37.6|  0.0|INFORM|Feature Type: `GCCSA_intermediate'
08 08:21:36|  37.6|  0.0|INFORM|Attribute(64 bit real)          : `AREA' has value `0'
08 08:21:36|  37.6|  0.0|INFORM|Attribute(64 bit real)          : `Dwellings' has value `42'
08 08:21:36|  37.6|  0.0|INFORM|Attribute(encoded: windows-1250): `GCC_CODE11' has value `79799'
08 08:21:36|  37.6|  0.0|INFORM|Attribute(encoded: windows-1250): `GCC_NAME11' has value `Migratory - Offshore - Shipping (NT)'
08 08:21:36|  37.6|  0.0|INFORM|Attribute(64 bit real)          : `Per_usu_re' has value `111'
08 08:21:36|  37.6|  0.0|INFORM|Attribute(encoded: windows-1250): `STE_CODE11' has value `7'
08 08:21:36|  37.6|  0.0|INFORM|Attribute(encoded: windows-1250): `STE_NAME11' has value `Northern Territory'
08 08:21:36|  37.6|  0.0|INFORM|Attribute(string)               : `fme_basename' has value `MB_2011_added_INFO'
08 08:21:36|  37.6|  0.0|INFORM|Attribute(string)               : `fme_dataset' has value `J:\Analysis\Census\Australia\3_Int\01_MB_with_metadata\MB_2011_added_INFO.shp'
08 08:21:36|  37.6|  0.0|INFORM|Attribute(string)               : `fme_feature_type' has value `MB_2011_added_INFO'
08 08:21:36|  37.6|  0.0|INFORM|Attribute(string)               : `fme_geometry' has value `fme_aggregate'
08 08:21:36|  37.6|  0.0|INFORM|Attribute(string)               : `fme_rejection_code' has value `AREA_REQUIRED'
08 08:21:36|  37.6|  0.0|INFORM|Attribute(entangled: string)    : `fme_type' has value `fme_area'
08 08:21:36|  37.6|  0.0|INFORM|          entangled to [SHAPE_GEOMETRY, SHAPE_GEOMETRY]
08 08:21:36|  37.6|  0.0|INFORM|Attribute(string)               : `multi_writer_id' has value `0'
08 08:21:36|  37.6|  0.0|INFORM|Coordinate System: `_LL-GDA94_0'
08 08:21:36|  37.6|  0.0|INFORM|Geometry Type: IFMEMultiArea
08 08:21:36|  37.6|  0.0|INFORM|Number of Areas: 0

Any suggestion?


texts: Labelling problem for Shapefile with Indic Fonts

texts:
On my computer (Ubuntu 14.10, 64 bit), I am using QGIS 2.6.

I have a shapefile, where the names of the features are written in Indic Script (using Google Transliterate). When I wish to label them, the labels appear as ??????? even if I use a Unicode font like (Lucida Sans Unicode).

I created the shapefile as a ISO-8859-1 type.


texts: how to get all features from a layer in openlayers?

texts:
I deal with geoserver and openlayers,
i use WMS layer for displaying due to loading time is shorter than vector layer.

My problem i couldn't get all features of a layer without event,
For sure i tray "getfeatureinfo" but it's return selected point not all
And i try to use vector layer for getting the features, But it not work unless i add it to the map ,
which will take long long time in loading .!

could any one help?


texts: I would like to use a hybrid metric to draw isochrones in sDNA. what are the possible ways, and are there any difference in the results?


texts: How to add legend to layer downloaded from arcgis.com?

texts:
I have downloaded a layer from arcgis.com that has a legend but when i view the layer in ArcMap, the legend is not there. How do i add the legend to the layer in ArcMap?
The map i have downloaded is: https://maps.owrb.ok.gov/arcgis/rest/services/Layers/Geology/MapServer
The legend is: https://maps.owrb.ok.gov/arcgis/rest/services/Layers/Geology/MapServer/legend


texts:
I've just observed something quite odd.

I have created a table with data for routing, and this table contains data for a graph over six different hour conditions. So when performing shortest_path, I select only the segments that corresponds to the "hour" column that I currently want.

I have created an index on the hour column but not on the source/target columns, and afterwards performed an analyze.

But when querying for shortest_path, I was NOT getting the shortest path, but a path where the route was lead through edges with extremely high weights (1,000,000 (created to tell not to go this way) where other weights are in the range 1-100).

After a lot of debugging (and blaming both my own graphs, pgrouting and everything else), I tried to perform a "Vacuum Analyze" on the table, and suddenly the results are correct.

So it seems like PostgreSQL had some problems that lead to pgRouting determining the wrong shortest path.

Now I ask you, if you have seen something like this before, and how to make sure this does not happen? The creation task of the graph table is automatized and I can see that analyze was actually performed after creation of the table.


texts: Display "values" or "description" of Value Map in Attribute Table

texts:
I am using a geopackage with a number of tables, some of which I have added in a Value Map to the fields.  The value map is to ensure consistent entry of codes to the data to comply with a specified schema by the client.

When editing via the Attribute Table, I would like to be able to flip between displaying either the "Value" or the "Description".

The reason is that if I have a few edits the drop down list of "Descriptions" is handy.  However if I want to use the update bar at the top of the table to do a bulk edit (eg Update All), I cant recall the a corresponding "values" and it create errors if I get it incorrect.

Is there a way to flip the display of the attribute table so that it alternates between the two options (values or description).

My usage is the creation of work as executed data.  I convert the surveyors informaiton into GIS and then add attributes to it.  Initially I would do an "bulk" update all bar (using the "values") to get the information say 90% correct, and then use the "Descriptions" to fine tune the remaining features.


texts: Storing scale & offset values from gdal_translate

texts: When using the -rescale option of gdal_translate, the output file does not have the scale or offset values stored in its metadata and therefore cannot be automatically unscaled when read. Is this intended behaviour or a bug?


texts:
I am trying to perform road network analysis with sDNA. I have successfully generated a shapefile with all centrality measures I need using the following command in the python shell

C:\Python27\python.exe "c:\Program Files\sDNA\bin\sdnaintegral.py" --im
"net=C:\...\inputfile.shp" --om "net=C:\...\outputfile.shp" start_gs=;
end_gs=;metric=ANGULAR ;radii=n;arcxytol=; arcztol=;custommetric=;disable=;;
origweight=;destweight=;weight_type=Link;cont;nojunctions;nohull

But what I really want is to use the python interface explained in the runcalculation.py file (which comes within the installation folder). In it, a function called runcalculation is defined. It takes six arguments: env, name, config_string, input_map, output_map and dll. I guess the following would work for the last four items:

config_string = 'start_gs=;end_gs=;metric=ANGULAR; radii=n; arcxytol=;
arcztol=; custommetric=;disable=;;origweight=;destweight=;weight_type=Link;
cont; nojunctions;nohull’
input_map = {‘net’:’ C:\...\inputfile.shp’}
output_map = {‘net’: C: \outputfile.shp}
dll=’’

but I can not infer the contents of env (environment object) and name.

Would you please provide me with an example to see how those two parameters should look like?
Will I be able to obtain all centrality measures at once as with sdnaintegral.py through this function runcalculation?


texts:
I work on a wood fire simulator project in university that uses QGIS.
I have to draw various symbols on the map.

As I'm not able to design it with qgis functions so I made it with in a stand alone python program using PYQT5 like that:

Can I use my PYQT drawing functions in QGIS?


texts: Adding multiple layers in GeoServer WMS URL

texts:
I have a URL that I am passing into an API call to bring back a static image of a map. This URL currently works:

http://computername:8080/geoserver/Packt/wms?service=WMS&version=1.1.0&request=GetMap&layers=Packt:countries&styles=&bbox=-181.800003051758,-90.8669281005859,181.800003051758,84.4676055908203&width=768&height=370&srs=EPSG:4326&CQL_FILTER=continent%20%3D%20%27North+America%27&format=application/openlayers#

In the book I am reading it says the layers parameter can 'insert a comma-separated list of layers' but when I try to add another layer like so:

http://computername:8080/geoserver/Packt/wms?service=WMS&version=1.1.0&request=GetMap&layers=Packt:countries,Packt:ne_10m_railroads&styles=&bbox=-181.800003051758,-90.8669281005859,181.800003051758,84.4676055908203&width=768&height=370&srs=EPSG:4326&CQL_FILTER=continent%20%3D%20%27North+America%27&format=application/openlayers#

the browser does not accept it. My end goal is to call two layers and filter one of the layers with the CQL syntax.

I am running GeoServer 2.11.0 on Apache Tomcat 8.5 on Windows Server 2012.


texts: Merge large number of polygons using QGIS

texts:
I have a series of polygons that are the output of another GIS program - I need to merge the polygons into a single polygon.  I have tried a couple of different methods with no luck.  I have tried to use the "Dissolve" tool, using the Dissolve all - however this led to an essentially empty output (most probably because there isn't a common attribute between all the polygons).  I have tried to use the calculator in the attribute table (there are around 350 000 polygons) to add a dummy column but it doesn't show in the Dissolve tool.

I am finding the merge tool does the job - to a point, but ultimately is extremely time consuming.  Is there a way to merge all polygons in a layer without too many steps?  Ultimately I just need the geographic area, the values/names etc are entirely irrelevant for what I'm currently doing.  I have tried to select all within the attribute table by using the "invert selection" option (a recommendation on another site) but this doesn't appear to do anything.  I have also read about a SelectPlus plugin, but no luck finding this one either.

Any suggestions greatly appreciated.

EDIT: Sorry it's a single file that I'm wanting to merge the polygons in (want the outline of all of them together, there are no gaps).

EDIT2: and while it's very slow, I've now found how to select all the polygons at once (using the radius selection option), which at least seems to be working if I do in stages.


texts: How to populate QTableWidget via plugin?

texts:
I am attempting to populate a QTableWidget by adding rows, this is done by the user clicking on the add button:

I have the following code to try and insert rows but nothing is happening:

def addButton_clicked():
    # objectName for add button = addButton
    # objectName for QTableWidget = tableWidget
    model = dockwidget.tableWidget
    rowCount = self.model.rowCount()
    self.model.insertRows(rowCount, 1)

add_button = self.dockwidget.addButton
add_button.clicked.connect(addButton_clicked)

The Group and Order columns were created from Qt Designer. I am trying to create rows in which the user can select a group from the ToC (in the Group column) and assign an integer (in the Order column). Is something like this possible?

I am open to other ideas.

The following code from @WKT creates a table with columns and x number of rows. I can populate the table with the names of the groups in the ToC. However, I only want to populate the first column. The code I run at the moment populates all columns. Here is the code:

qTable = self.dockwidget.tableWidget
data = []
group1 = root.findGroup('Group1')
group2 = root.findGroup('Group2')
for child in group1.children():
    data.append(child.name())
for child in group2.children():
    data.append(child.name())
nb_row = len(data)
nb_col = 2
qTable.setRowCount(nb_row)
qTable.setColumnCount(nb_col)

for row in range(nb_row):
    for col in range(nb_col):
        item = QTableWidgetItem(str(data[row]))
        qTable.setItem(row,col,item)

qTable.setHorizontalHeaderLabels([u'Column1',u'Column2'])
qTable.resizeColumnsToContents()


texts: How do you create a map that has choropleths that will change with a timeslider?

texts: I have data for different countries over time that I would like to display in a choropleth map with a time slider that plays automatically on click.  What is the best way to do this.


texts: Output polygon data from postGIS and read it on Google Maps API

texts:
So, I have a group of polygons stored in PostGIS. I'm trying to build a web application with Google Base maps (using Google Maps API). Now, How do i output the polygons in the right format so that Google Maps API can read it?

I know about ST_Askml which outputs the data as KML, but what is the Google Maps function to read a kml on the fly? Or is there an alternative to outputting it as kml? Can i pull out the co-ordinates of each point on a polygon and load it onto the API on the fly?


texts: Variations in the number of Sub-Watersheds delineated by various tools

texts:
I have demarcated sub-watersheds using different tools like ArcHydro, ArcSWAT, and Hydrology tool (ArcGIS). The same input data have been used in all tools. There are 35 watersheds derived using ArcSWAT, 44 from ArcHydro and 50+ from the hydrology tool (ArcGIS).

Why are there are variations in the number of watersheds?

Which one is correct?

The area of the watershed is 1960 square Km.

How many sub-watershed can be there in this size of the watershed?

Is there any way to define the required number of sub-watersheds with the specific area range?

Because the 30+ sub-watershed for this watershed with 1960 square KM. seems too many.


texts: Clipping a road layer by using region attributes

texts:
I've been trying to clip a road layer into several layers by using region polygons. So in detail: I've got a layer with all roads of Germany and a second layer with all states and regions of Germany. And now I'm trying to divide this road layer by German states and regions and save them as shapefiles named after that region. Meaning: Loading this new shapefile only contains the roads of one region. It only works for one feature at a time, selecting the one region and going into the geoprocessing tools and clip the roads layer.

Is there any easy one-step-way to do it for all features? I could not find anything similar in the forum, as I don't know where and what to search for exactly.

If you need a more specific explanation, I'll do my best.


texts: How do I convert Geo-Coordinates to Pixel-Coordinates on a Map?

texts:
I have an image file that is a Robinson Projection Map of the world. I have a list of geo-coordinates. I would like to find the pixel-coordinates on the map that correspond to my the geo-coordinates I have. What is the best way of doing this?

Note that I'm developing a tool for this purpose --- I would like to be able to find the corresponding pixel-coordinate for any geo-coordinate. So I'm trying to create a function to project geo-coordinates to pixel-coordinates. I'm wondering if there's a standard way of doing this, or if I have to adjust for particularities in map distortion, etc.

Can anyone help?


texts: geojson data not rendering with openlayers in geodjango

texts:
map.html javascript code
var map = new ol.Map({
        layers: [
          new ol.layer.Tile({
            source: new ol.source.OSM()
          })
        ],
        target: 'map',
        view: new ol.View({
          center: [0, 0],
          zoom: 2
        })
      });

       var vectorSourceJsonp  = new ol.source.Vector({
        format: new ol.format.GeoJSON(),
        loader: function(extent, resolution, projection) {

            $.ajax({
                        <!-- url: "{% url 'geojson_provinces' %}", -->
                        url:'/geojson_provinces/',
                        dataType: 'json',
                    });
                },
                strategy: ol.loadingstrategy.createTile(new ol.tilegrid.XYZ({
                    maxZoom: 19
                })),
                projection: 'EPSG:4326'
            });
             // Executed when data is loaded by the $.ajax method.
            var loadFeatures = function(response) {
                vectorSourceJsonp.addFeatures(vectorSourceJsonp.readFeatures(JSON.parse(response)));
            };
            // Vector layer
            var vectorLayerJsonp = new ol.layer.Vector({
                source: vectorSourceJsonp,
                projection : 'EPSG:4326',
                style: new ol.style.Style({
                    stroke: new ol.style.Stroke({
                        color: 'green',
                        width: 2
                    })
                })
            });
            var mapJsonp = new ol.Map({
                target: 'mapJsonp',
                renderer: 'canvas',
                layers: [osmLayer, vectorLayerJsonp],
                view: new ol.View({
                    center: [-75.923853, 45.428736],
                    maxZoom: 19,
                    zoom: 11
                })
            });

views.py
def geojson_provinces(request):
    conn = psycopg2.connect(dbname="geodjango",host='localhost',user='postgres', password='postgres', port=5433)
    dict_cur = conn.cursor(cursor_factory=psycopg2.extras.DictCursor)
    cur=conn.cursor()
    res = dict_cur.execute("SELECT ST_AsGeoJson(geom) AS json FROM pft_pft LIMIT 10;")

    points=dict_cur.fetchall()
    print (points)
    json_res  = []
    for row in points:
        json_obj  = row['json']
        print (json_obj)
        json_res.append(json_obj)
        print(json_res)
    #json = simplejson.dumps(data)

    dict_cur.close()
return render_to_response('map.html',
                          {'json_res': json.dumps(json_res)},
                         context_instance=RequestContext(request))

urls.py
url(r'^geojson/$', geojson_provinces, name='geojson_provinces'),

only map is getting displayed . When I use the debugger in chrome this is what i get

and finally i get Uncaught TypeError: ol.tilegrid.XYZ is not a function
please help! I even tried changing the url to "{% url 'geojson_provinces' %}", but not working


texts: How to select features that equals to 159 different values

texts:
I am new to QGIS. My question is, e.g, I have a .shp that has 10,000 features: .

how to select features where the column of id that equals to 159 different values, e.g., selecting
id=63703,63706,63708,63713,63733,63783,63705,63702,63745,65713,66733,64792....?


texts: Generating Image from two vector files in sequence using Python

texts:
I have two folders in my system containing multiple vector files in it, i have named them as A1,A2, A3........., B1,B2,B3..... and so on. I want to map the A1 file of first folder with the A1 file of second folder but both files should have different symbology and then to export them as a image in full extent.
I have tried writing a Python code for that but it is giving the last file only as output.
In first folder it has file which contain overlay vector file in a grid and in second folder it has file which contains within vector file of the grid
My code is this:
> import os
import processing
>
> for withinfiles in os.listdir (directory):
>     if withinfiles.endswith('.shp'):
>         inputwithinfile = directory + "/" + withinfiles
>         #splitting the filename for naming folders
>         inputwithinfilename, file_extention = os.path.splitext(withinfiles)

    for overlayfiles in os.listdir (directory):
        if overlayfiles.endswith('.shp'):
            inputoverlayfile = directory + "/" + overlayfiles
            #splitting the filename for naming folders
            inputoverlayfilename, file_extention = os.path.splitext(overlayfiles)

            if inputwithinfilename == inputoverlayfilename:
                withinlayer = QgsVectorLayer(inputwithinfile, "Within layer", "ogr")
                QgsProject.instance().addMapLayer(withinlayer)
                withinlayer.loadNamedStyle(qml file directory)

                overlaylayer = QgsVectorLayer(inputoverlayfile, "overlay layer", "ogr")
                QgsProject.instance().addMapLayer(overlaylayer)
                overlaylayer.loadNamedStyle(qml file directory)

                vlayer = iface.activeLayer()
                settings = QgsMapSettings()
                settings.setLayers([withinlayer , overlaylayer])
                settings.setBackgroundColor(QColor(255, 255, 255))
                settings.setOutputSize(QSize(800, 600))
                settings.setExtent(vlayer.extent())

                image_location = directory  + "/" + inputwithinfilename + '.png'
                render = QgsMapRendererParallelJob(settings)

                def finished():
                    img = render.renderedImage()
                    # save the image; e.g. img.save("/Users/myuser/render.png","png")
                    img.save(image_location, "png")

                render.finished.connect(finished)

                render.start()
            else:
                continue


texts:
How do you automatically split all loop links in a table with line objects?

A loop being one where first and last coordinate is the same.

Should be split at 50%.


texts: Using GeoWebCache as offline layer server?

texts:
I'm trying to use a Geoserver external WMS service in a offline environment.

First I put my own geoserver online and configure a WMS store pointig to the desired source.  Next I publish the desired layer and configure Create a Cached Layer for this Layer. Next I go to GWC and seed all tiles from a bounding box from zoom level 0 to 18. I got near 10G of images.

Now I just turn off my router and go to:

geoserver/gwc/demo/myws:mylayer?gridSet=EPSG:4326&format=image/png8

but got this error:

Failed to execute request http://ORIGINAL_SOURCE_WMS/?REQUEST=GetCapabilities&VERSION=1.3.0&SERVICE=WMS`
Caused by: org.apache.commons.httpclient.ConnectTimeoutException: The host did not accept the connection within timeout of 30000 ms

As I can see the Geoserver still depending the connection.

EDIT
I found this post where Ian said:

  If you create your tiles using GeoWebCache GWC (and other tile caches
  I expect) and run the preseeding operation to create all the tiles on
  disk you can then copy those tiles to a remote machine and use the
  OpenLayers XYZLayer to access the tiles.

This can help me too  since I've done all steps he describe but I'm using Cesium instead OpenLayers. Some tip about this?

Added the Cesium tag.


texts: Print functions for data projection and rename file

texts:
I have been struggling with functions in this script as I would like to see the projection in the file and then rename. The problem is that it causes and definitions are not done properly

Here is the long code:

# -- coding: utf-8 --

# Exercise 1 - rename and redefine

# ===================================== Imports ======================================= #

# import all necessary libraries here
import arcpy
import os
import sys
import glob
import shutil
from os import path
from arcpy.sa import *

# ===================================== Functions ======================================= #
# define your function(s) here

# finding the british grid projection
def british_grid_projection(in_shp):
  # input data is in NAD 1983 UTM Zone 11N coordinate system
  input_features = r"C:\Users\Leon\prog_gis_coursework\part_a_gis_functions\test_data\London_Borough_Excluding_MHW.shp"
  # output data
  output_feature_class = r"C:\Users\Leon\prog_gis_coursework\part_a_gis_functions\test_dataLondon_Borough_Excluding_MHW.shp"
  # create a spatial reference object for the output coordinate system
  out_coordinate_system = arcpy.SpatialReference('')
  # run the tool
  arcpy.Project_management(input_features, output_feature_class, out_coordinate_system)
  print ( "brit_grid_projection_processing")

def rename_files (brit_grid):
    # retrieve file names from the list
  file_list = os.listdir(r"C:\Users\Leon\prog_gis_coursework\part_a_gis_functions\test_data")
  print (file_list)

    # renaming the file

def main(duplicate):
    # make a duplicate of an existing file
  if path.exists("London_Borough_Excluding_MHW.shp"):
# get the path to the file in the current directory
    src = path.realpath("London_Borough_Excluding_MHW.shp");

# rename the original file
    os.rename("London_Borough_Excluding_MHW.shp","ldnboroughs_britgrid.shp")

# ====================================== Main script ====================================== #
# the execution starts here.
print ("Script starting")

# set up environment
arcpy.env.overwriteOutput = True
arcpy.env.workspace = r"C:\Users\Leon\prog_gis_coursework\part_a_gis_functions\test_data"

# using the try tool

try:
    # set local variables
    in_dataset = "London_Borough_Excluding_MHW.shp"

    # get the coordinate system from the shapefile
    dsc = arcpy.Describe("London_Borough_Excluding_MHW.shp")
    coord_sys = dsc.spatialReference

    # run the tool
    arcpy.Project_management("London_Borough_Excluding_MHW.shp", "ldn_britgrid")

    # print messages when the tool runs successfully
    print(arcpy.GetMessages(0))

except arcpy.ExecuteError:
    print(arcpy.GetMessages(2))

except Exception as ex:
    print(ex.args[0])

# Calling all 3 functions here

british_grid_projection('in_shp')
rename_files('brit_grid')
main('duplicate')

print( "Script finished!")

Here is the error codes

Script starting
Failed to execute. Parameters are not valid.
ERROR 000735: Output Coordinate System: Value is required
Failed to execute (Project).

Traceback (most recent call last):
  File "C:\Users\Leons\Documents\part_a_gis_functions\folder for export\part a\Question 1 - new draft.py", line 84, in <module>
    british_grid_projection('in_shp')
  File "C:\Users\Leons\Documents\PGIS\part_a_gis_functions\folder for export\part a\Question 1 - new draft.py", line 29, in british_grid_projection
    out_coordinate_system = arcpy.SpatialReference('NAD 1983 StatePlane Nevada West FIPS 2703 Feet')
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\arcobjects\mixins.py", line 949, in __init__
    self._arc_object.createFromFile(item)
RuntimeError: ERROR 999999: Error executing function.


texts: Which coordinate system is used in "gstat" package in R

texts: I am confused about which coordinate system is used in gstat package in R-language and also in geoR package in R-language.


texts: Creating perpendicular line transects in PostGIS

texts:
I have a river network (Line with arrow as shown in image below) and I need to create perpendicular lines to the river at every 15 meters as you walk along it (see image below).
Is there a function in PostGIS capable of doing this?


texts:
I'm using version 4.0 of the ArcGIS API for my maps.
I have this problem with maps whenever the page has to be scrolled down to view the map.
My map is somewhere in the middle of the page so the user has to scroll down to view it. Then if a certain point was clicked on the map or if a certain area in dragged to zoom in, the pointers the map shows as selection are not the exact points the mouse is pointing at. it points to a location above the mouse pointer.

Here is a sample code.

    <!DOCTYPE html>

  <meta charset="utf-8">
  <meta name="viewport" content="initial-scale=1,maximum-scale=1,user-scalable=no">
  <title>Add the Compass widget to a basic 2D map - 4.0</title>

  <link rel="stylesheet" href="https://js.arcgis.com/4.0/esri/css/main.css">

If you run this code you'll see if you scroll down and drag an area (by holding shift and dragging the mouse) it will drag an area above your selection.


texts: Attribute names being truncated during extraction

texts:
I have 30-40 raster layers and each raster has approximately 40 columns. The attribute/column names are very long, and the beginning of each column name is the same. During the extraction process, these column names are truncated, losing a lot of information regarding the variable/column name.

I am, for example, using Extract multi values to points to extract from .tiff to a shapefile point feature class.

For example:

vpr_Avg_1970to2010_m01Jan_800m is truncated to vpr_Avg_19
vpr_Avg_1970to2010_m02Feb_800m is truncated to vpr_Avg_20

.... vpr_Avg_21

Is there any way to keep the entire column name?

I'm working with ArcMap 10.1


texts: Which coordinate system/ projection would you use to measure lines globally distributed?

texts: Which geographic coordinate system /projection would you recommend to measure the length of hundreds of lines which are globally distributed (southern and northern hemisphere)? The lines have a length up to 200km.


texts:
I am trying to add layers from the fields of a reference layer using select by location and a cursor.  Also, want to prompt the user with a message box asking if they want to add the selected by location.  I am able to print string value from the reference layer when the user responds with "Yes" with the message box, but get an error when passing the string value as a path using "arcpy.mapping.addLayer".

def onMouseDownMap(self, x, y, button, shift):
    mxd = arcpy.mapping.MapDocument("CURRENT")
    df = arcpy.mapping.ListDataFrames(mxd,"*")[0]
    #layer path
    lyr1 = arcpy.mapping.Layer("C:\ArcGis\Aerial_Imagery.lyr")
    #Add Layer
    arcpy.mapping.AddLayer(df, lyr1, "TOP")

    point = arcpy.Point(x, y)
    ptGeometry = arcpy.PointGeometry(point)
    #Select layer by points/mouse down on map
    arcpy.SelectLayerByLocation_management("Aerial_Imagery.lyr", "INTERSECT", ptGeometry)
    #cursor to display selected layers on field name "Layer_Path"
    cursor = arcpy.SearchCursor("Aerial_Imagery.lyr")
    field = "Layer_Path"
    for row in cursor:
        selpath = (row.getValue(field))
        message = "Would you like to add this Imagery:" + selpath
        Addbutton = pythonaddins.MessageBox(message, "Imagery Finder", 4)
        if Addbutton == "Yes":
            addlyr = arcpy.mapping.Layer(selpath)
            arcpy.mapping.AddLayer(df, addlyr, "TOP")#this breaks

            print selpath #this works

        else:
            pythonaddins.MessageBox("this is the answer to NO", "no answer", 0)

   #Remove Layer
    for lyr in arcpy.mapping.ListLayers(mxd, "Aerial_Imagery.lyr"):
        #print "Remove Layer", lyr
        arcpy.mapping.RemoveLayer(df, lyr)

    #Refresh view and TOC
    arcpy.RefreshActiveView()
    arcpy.RefreshTOC()


texts: Retaining UNIQUE attribute fields after dissolving a shapefile

texts: I tried the solution to "Dissolving shapefiles but retain attribute fields" but my original shapefile (3000 rows in attribute table) has many duplicate values in the field that I want to retain so that the merge overruns the field width. Is there a way to merge only unique field values in the merged attribute field?


texts: Returning state of point from postgis / mapnik pgsql db

texts:
I realize that I can confirm whether or not a point lies within a state like this:

            SELECT ST_Contains(poly,pt) As testr
            FROM
            (
                    SELECT ST_GeomFromText('POINT(-11817883.86 4424094.21)',900913) As pt,
                    (SELECT way FROM planet_osm_polygon WHERE admin_level='4' and name='New Mexico') As poly
            ) As foo;

However, I'm looking for a more direct way to find what state a point lies in.

Any suggestions?

Taking @Paul's advice, I constructed the following query that worked surprisingly fast:

SELECT name
FROM planet_osm_polygon As poly
WHERE admin_level='4' AND ST_Contains(poly.way,ST_GeomFromText('POINT(-11820276.11 4441642.75)',900913));


texts: How to convert IPUMS-NHGIS Shapefiles to GeoJSON files using R, and then display using Leaflet?

texts:
The NHGIS website has a fantastic wealth of US boundary shapefiles. These files can be filtered and downloaded here.

My problem is that after bringing these shapefile(s) into R and then converting them into GeoJSON files, the GeoJSON files do not display when I add them to a Leaflet object/map.

I have gone out and found other shapefile data of the same census year, used the process mentioned above, and then successfully displayed the GeoJSON file on a Leaflet map.

So I am confused as to why the NHGIS shapefiles are not converting correctly, while other shapefiles from other sources are. Both sources have all the necessary files for shapefiles (.DBF, .PRJ, .SBN, .SBX, etc...).

My code (that does not display on a Leaflet) is listed below, along with an image of a successful map with counties in 1790.

A helpful shapefile or GeoJSON converter can be found here; however, after using this GUI converter I still can't get the NHGIS shapefiles to convert correctly.

I am using a shapefile from NHGIS Data Finder and filtering to 1790 and requesting the STATE shapefile.

library(tidyverse)
library(magrittr)
library(rgdal)
library(spdplyr)
library(geojsonio)
library(rmapshaper)

# Import Shapefile into R.
us = rgdal::readOGR(dsn="nhgis0004_shape/nhgis0004_shape/nhgis0004_shapefile_tl2000_us_county_1790",
                layer="US_state_1790", stringsAsFactors = FALSE)
# Filter to get only up to 1790    us %<>%
dplyr::mutate(START_YEAR = lubridate::year(us[["START_DATE"]]),
             END_YEAR = lubridate::year(us[["END_DATE"]]))
us %<>%
   dplyr::filter(END_YEAR > 1790 & START_YEAR <= 1790)

# Convert SP Data Frame to GeoJSON.
us_json <- geojson_json(us)

# Simplify the geometry information of GeoJSON.
us_json_simplified <- ms_simplify(us_json)

# Keep only the polygons inside the bbox (boundary box).
us_json_clipped <- ms_clip(us_json_simplified, bbox = c(-170, 15, -55, 72))

# Save it to a local file system.
geojson_write(us_json_clipped, file = "~/Downloads/us.geojson")

leaflet::leaflet() %>%
   leaflet::addTiles() %>%
   leaflet::setView(lng = -98.583, lat = 39.833, zoom = 3) %>%
   leaflet::addGeoJSON(us_json_clipped)


texts: geojson very slow to load in openLayers3

texts:
I'm new to openLayers. I have parcel data stored in a geojson file, and it's pretty huge - there are about 170,000 records. I would load the data as a WMS, but I want the features to highlight in a different color on mouseover. I'm following this example:

http://openlayers.org/en/v3.2.1/examples/vector-layer.html?q=geojson

I'm assuming I can't apply this functionality with a WMS layer? But the performance is just too slow with the geojson. I've also tried setting maxResolution and minResolution for the vector layer but it doesn't seem to make a difference - it also doesn't seem to work no matter what I set the values to (but I could be doing something wrong - I don't yet know how the resolutions correspond to zoom levels). Does anyone have any suggestions on how I can make this work?


texts: Iterating rasters from two different folders using ArcPy?

texts:
I have created a simple script like that, which is working fine

import arcpy
from arcpy.sa import *
arcpy.CheckOutExtension("spatial")
arcpy.env.workspace = r'C:/2017/wgs84/project'
direction = arcpy.ListRasters("*WIND*","TIF")
print direction
speed = arcpy.ListRasters("*SPEED*","TIF")
print speed
for k in range(len(direction)):
  outCon = Con(Raster(direction[k]) == 1, Raster(speed[k]), 0)
  outCon.save('speed0' + str(k) + '.TIF')
  print ('done')
print ('done final')

I have different wind direction files(rasters) in different folders and speed files in another one. This code works for what I need, only if I have pasted the speed files and the wind files in the same folder.
What if I wanted to create a list from one folder and a list from another so as not having to copy paste every time?

I tried to define also a different workspace for direction and different for speed but apparently its not working(maybe the last env.workspace is "hold" in the memory?)

arcpy.env.workspace = r'C:/speed'
speed = arcpy.ListRasters("*SPEED*","TIF")
print speed
arcpy.env.workspace = r'C:/2017/wgs84/project'
direction = arcpy.ListRasters("*WIND*","TIF")


texts:
I found a file that contains about 85,000 messages (3 minute live AISHub feed), but I can't make heads or tails of it.

http://www.aishub.net/nmea-sample.html

!AIVDM,1,1,,A,13aEOK?P00PD2wVMdLDRhgvL289?,0*26
!AIVDM,1,1,,B,16S`2cPP00a3UF6EKT@2:?vOr0S2,0*00
!AIVDM,2,1,9,B,53nFBv01SJ<thHp6220H4heHTf2222222222221?50:454o<`9QSlUDp,0*09
!AIVDM,2,2,9,B,888888888888880,2*2E

Is there a schema somewhere?


texts: Trimble pathfinder data dictionary default value exclusion

texts:
Is it possible to have a numeric attribute in trimble pathfinder without having to use a default value? I am trying to make a numeric attribute compulsory to fill in BUT because the data dictionary editor forces the use of a default value such as zero, the software automatically treats the attribute as being logged instead of alerting the field capturer if he or she has missed a field if they accidentally do not answer the question in the field.

This means that setting a numeric attribute to "required" doesn't mean anything in the field because if the value is not filled in it just accepts the default as a field entry anyway.


texts: Splitting Line.shp at Point.shp with Python/ogr

texts:
I'm trying to modify some ESRI Shapefiles using Python with gdal/ogr module.
What I want to do is quite simple, but I didn't find any solution by now.

I have two shapefiles: (1) Lines.shp and (2) Points.shp
The features of the Points.shp mark intersections of the lines in the Lines.shp

Is there any way to split the line-features at the given points?

edit: It should work like the "Split Line at Point"-tool from ArcMap
(Data Manangement Tools --> Features --> Split Line at Point)

edit2: I already knew about the other topic, but I don't want to use shapely. Also I don't want to save my Points into a list. I just wanted to know, if tehre is any function/method like ESRIs "Split Line at Point"


texts: Creating a buffer around a raster in Google Earth Engine

texts:
I want to create a 10km buffer around mangrove patches from the mangrove forest raster. I've tried a number of approaches but can't get anything to work. I tried converting the raster to vector but at the resolution (30m) I need it is too big and the computation times out. I also tried to the focal.max tool but this seems to be dependent on the zoom level I'm using and therefore capture only certain patches

//import Mangroves of the World
var MangroveImage = ee.ImageCollection('LANDSAT/MANGROVE_FORESTS')

//select the image from the mangrove dataset
var mangrove = ee.Image(MangroveImage.first());

//buffer around mangroves
var mangrove_buffer = mangrove.focal_max({
radius: 10000,
kernelType: 'circle',
units:'meters',
iterations: 1
})


texts: Self-hosted WeoGeo-style geospatial library software

texts:
I've been hearing a lot about WeoGeo Library recently, but (a) it costs and (b) it is hosted 'in the cloud', which is not acceptable for the data my organisation is using because of data usage regulations.

Is there anything similar which I can run on my own server within my organisation? The sort of features that I'm looking for are basically managing a library of geospatial data (in various formats - geodatabase, shapefiles, images) and allowing them to be searched, queried and visualised.

Any ideas?


texts: Importing large number of shapefiles into feature class using ArcPy with column mapping?

texts:
I'd like to write an ArcPy script to import a large number of shapefiles into an EXISTING feature class. I also need to be able to perform column mappings as there may be different column names. I need to write a script so I can automate this process.

I'm not looking for a complete script but rather advice as to where I should be looking to perform this task. I'm fairly new to ArcGIS.


texts: Is it possible to somehow load feature layer from REST (like: ArcGIS Example Service) to ArcMap? No need to edit that layer, we just need to select data from it and then create a local FeatureClass.


texts: Where can I get aerial maps for offline use?

texts:
I'm looking for somewhere I can get aerial map files for exclusive offline use. Essentially I need aerial files (JP2K, GeoTIFF, etc.) that can be purchased or acquired and hosted on a private offline server.

I have used NAIP aerial files in a past and they have worked great but they don't have coverage in the area I am interested in. I need the files to have at least the equivalent resolution as NAIP. It also can't be a subscription service it needs to be a one time payment to acquire the files that I need.

Preferred resolution is equivalent to NAIP so 1-meter ground sample distance (GSD), I would like the files to be <5 years old but that can be a little flexible, and the area of interest is Aberdeen Proving Grounds outside of Baltimore.


texts:
I am very new to R, so please bear with my explanation. I have a raster list made up of 671 objects that I am hoping to clip to a shapefile. I don't know if this can be done, so any suggestions would be really appreciated.

At this link I explain how I uploaded the original .bin files and got help making a sub-list of just the March and April dates (I'm working with sea ice concentration rasters from NSIDC): Create a new list from a raster list in R

Now, I am hoping to crop that new sub-list. Here is what I have done (where "list_mar_apr" is the sub-list and "boundary" is the shapefile):

# Import and check shapefile
boundary <- readOGR("file_folder/shapefile.shp")
plot(boundary)

# Crop subset to the shapefile
list_mar_apr_crop <- crop(list_mar_apr, boundary)

I keep getting an error message that reads: "unable to find an inherited method for function 'crop' for signature 'list'.


texts: How do you read GeoTiff-file in SAGA GIS Program?

texts:
I'm very new with all GIS-stuff and I got advice that SAGA GIS would be a good program to start with. I would like to read a GeoTiff-file in SAGA, but I can't find an "import file"-function x)

Any advice?


texts: Convert bbox in meters (stereographic projection)

texts: I need to convert the bbox H-SAF area [25-75°N lat, 25°W-45°E long] in meters in the stereographic projection (EPSG:3413 WGS 84 / NSIDC Polar Stereographic North). How can I solve this?


texts: Can the OpenLayers QGIS plugin be modified to have better map loading performance? I use only Google satellite/maps.


texts:
I am starting with maps but im a beginner. The answer is probably easy, but i could not find it anywhere.
I cannot figure out why to use maxExtent. If I delete it from my application, its behavior doesn't change or I do not notice. So why should I use it? I understand purpose of restrictedExtent but not of maxExtent. Can you please somebody please explain it to me, thank you.

var extent = new OpenLayers.Bounds(-2003750.34, -2003750.34, 2003750.34, 2003750.34);

map = new OpenLayers.Map("map", {
        controls: [
                        new OpenLayers.Control.Navigation(),
                        new OpenLayers.Control.PanZoomBar(),
                        new OpenLayers.Control.ScaleLine(),
                        new OpenLayers.Control.Attribution(),
                        new OpenLayers.Control.LayerSwitcher()],
        maxExtent: extent,
        //restrictedExtent: extent,
        maxResolution: 156543.0399,
        numZoomLevels: 19,
        units: 'm',
        projection: new OpenLayers.Projection("EPSG:900913"),
        displayProjection: new OpenLayers.Projection("EPSG:4326")
    });

var osm = new OpenLayers.Layer.OSM.Mapnik("Mapnik", { 'attribution': 'Provided by blop' });
var gmap = new OpenLayers.Layer.Google("Google Streets", { visibility: false });

map.addLayers([osm, gmap]);


texts: MODIS-Sphere to WGS84 coordinate system in ArcGIS

texts:
I want to overlay two global raster datasets on land use which use different coordinate systems (MODIS-Sphere & WGS84) using ArcGIS Desktop. There seems to be no direct transformation option from MODIS-Sphere to WGS84. It seems to me that I would have to use a custom transformation. However, for that I would need some sort of transformations dialogue/script which I don't have. Or am I missing something?!

Does anyone have any idea on which steps to take?

The dataset I want to transform is http://urban.yale.edu/sites/default/files/data/setoguneralphutyra_uf.zip

Format: ESRI GRID
Projection: Goode's Homolosine Equal-Area Projection
Geographic Coordinate System: MODIS-Sphere

I want to transform this dataset so that it matches the other one (which uses GCS_WGS_1984 as spatial reference and D_WGS_1984 as datum)


texts: How to create/add features without geometry?

texts:
I have a simple SpatiaLite geodb loaded in QGIS with one table. There is a POINT geometry field and some alphanumeric fields.

In some cases I can record data about a specific site, but the coordinates are not available at this time and will be retrieved and recorded later.

How do I add features to this layer without adding their geometries?


texts: QGIS atlas: make only current polygon visible

texts:
I have a layer with lots of polygons and want an atlas that goes through each polygon. However, the atlas page should only show the current polygon and not its neighbors. In the following image only the rectangle with the red lines in it should be visible and the others should be invisible.

Any ideas how to do that?


texts: Change GeoServer admin password via REST

texts:
I have an architecture where I deploy GeoServer instances automatically and I would like to reset the admin password to something random for each instance.  I am running GeoServer 2.16.2 on Debian 10.3. The REST documentation suggests that I can change passwords through a push to the /rest/security/self/password path. When I do this, the GeoServer log seems to indicate that the password was changed, yet I can still log in with the default password but can't with the new password.

For example, first I change the password:

>>> curl -u admin:geoserver -X PUT http://localhost:8080/geoserver/rest/security/self/password -H  "accept: application/json" -H  "content-type: application/json" -d "{  \"newPassword\": \"test\"}"

The log indicates the change was successful:

29 Apr 20:44:48 INFO [geoserver.security] - Start reloading user/groups for service named default
29 Apr 20:44:48 INFO [geoserver.security] - Reloading user/groups successful for service named default
29 Apr 20:44:48 INFO [geoserver.security] - Start reloading user/groups for service named default
29 Apr 20:44:48 INFO [geoserver.security] - Reloading user/groups successful for service named default
29 Apr 20:44:48 INFO [security.xml] - Successful lock: security/usergroup/default/users.xml.lock
29 Apr 20:44:48 INFO [geoserver.security] - Start storing user/groups for service named default
29 Apr 20:44:48 INFO [geoserver.security] - Storing user/groups successful for service named default
29 Apr 20:44:48 INFO [geoserver.security] - Start reloading user/groups for service named default
29 Apr 20:44:48 INFO [geoserver.security] - Reloading user/groups successful for service named default
29 Apr 20:44:48 INFO [geoserver.security] - Adjusted last modified for file: security/usergroup/default/users.xml
29 Apr 20:44:48 INFO [geoserver.security] - Start reloading user/groups for service named default
29 Apr 20:44:48 INFO [geoserver.security] - Reloading user/groups successful for service named default
29 Apr 20:44:48 INFO [geoserver.security] - Adjusted last modified for file: security/usergroup/default/users.xml
29 Apr 20:44:48 INFO [geoserver.rest] - Changed password for user admin

But when I do something that requires authentication with the new password, like list the layers:

>>> curl -u admin:test -X GET http://localhost:8080/geoserver/rest/layers -H  "accept: application/json"

I get an error on the server:

29 Apr 20:46:27 WARN [geoserver.security] - Failed login, user admin from 172.17.0.1
29 Apr 20:46:27 INFO [geoserver.security] - Brute force attack prevention, delaying login for 1385ms

But the original password works fine:

>>> curl -u admin:geoserver -X GET http://localhost:8080/geoserver/rest/layers -H  "accept: application/json"

{"layers":{"layer":[{"name":"tiger:giant_polygon","href":"http:\/\/localhost:8080\/geoserver\/rest\/layers\/tiger%3Agiant_polygon.json"},{"name":"tiger:poi","href":"http:\/\/localhost:8080\/geoserver\/rest\/layers\/tiger%3Apoi.json"},{"name":"tiger:poly_landmarks","href":"http:\/\/localhost:8080\/geoserver\/rest\/layers\/tiger%3Apoly_landmarks.json"},{"name":"tiger:tiger_roads","href":"http:\/\/localhost:8080\/geoserver\/res...

What can I do to change my GeoServer admin password programatically?


texts: Multiple Ring Buffer Tool Stopped Working

texts: The Multiple Ring Buffer tool in ArcGIS Desktop 10.6.1 stopped working.  It worked fine, now won't work at all.  It doesn't matter where I send the output, I get the "Error 000210: Cannot create output".  Again, it doesn't matter where I send the output, what I call it or anything.  Even the defaults don't work.


texts:
I am using the PostGIS Manager Load data to PostGIS tool on a Windows install and am getting the following error:

Something's wrong:
["'shp2pgsql' is not recognized as an internal or external command,\n", 'operable program     or batch file.\n']

Does anyone have any ideas as to what is wrong?  I am trying to append to an existing table.

Thanks...


texts: Loading locations faster in ArcGIS Network Analyst?

texts:
I am currently using ArcGis Network Analyst to calculate an OD cost matrix.

I was wondering whether there is a faster way to load my locations (origins, destinations and barriers) because my current rate is around 10 locations loaded per second (which is quite slow for I have about 7.000 points to load).

Moreover, is there a way to save these locations in order to use them later?

I think of origins/destinations but also of the barriers. I didn't find anything about this in the ArcGIS documentation.


texts:
I'm trying to determine a way to accurately produce groundwater contours in ArcGIS. The industry standard is using Surfer to construct the contours from monitoring well data. I would imagine ArcGIS Spatial Analyst or maybe Geostatistical analyst could perform this function as well.

Has anyone found a method to construct groundwater contours using ArcGIS? How accurate is the process?*


texts: Problem with the attribute table values, they are replaced with ************

texts: I suddenly have a problem with my attribute table. When I change a value the whole line is messed up because values shift and the value I add is replaced by ***********. How can I avoid this from happening?


texts: Esri ArcGIS JavaScript - printTask.execute issues with error

texts:
I'm using the ArcGIS JavaScript API 3.18.

The print service is configured to be async, with a timeout of 10 minutes.

The code to call the job is as follows:

  printTask.execute(printParameters,
      function (result) {
          $log.debug('Esri portion of print task complete...');
          var exportedMapUrl = result.url;
          onComplete({ url: exportedMapUrl });
      },
      function (errorResult) {
          $log.info(errorResult);
          onError({ message: errorResult, details: [] });
      });
  });

If I intentionally cause an error, i.e. by disabling the print service, the error callback is triggered and all is good, however if I execute a complex print task, which takes >10 minutes (that's another story) and the task times out, nothing is triggered. i.e. if I look in the network tab of Chrome dev tools, I can see that the arcgis api is polling every second or so, and eventually the status
jobStatus:"esriJobFailed"
is returned, however this seems to be eaten by the arcgis api, and nothing is raised in my code - how can I trap a timeout? This seems to be a pretty fundamental "error" that just disappears and it's difficult for me to surface this to the user.

Edit

This is just using the standard print service, with slightly modified templates. The process succeeds when A4, A3 & A1 are selected, and also when a less complex A0 map is printed (i.e. less layers).