texts:
I am creating pie charts of what fish have eaten at several stations, and unfortunately the colours of the different slices are not the same across the different stations, i.e. they get assigned a new colour scheme each time despite the variables being the same at each station.

Is there a way of standardising so that all pie charts have the same color for the same variable, it seems like a very obvious requirement in many cases, as they all are being filtered from the same original dataset?

the only workaround i have found so far is to format one of the stations the way i want it, and creating copies of this and specifying a different data range afterwards, but this seems very backwards.

I am not the best at the coding language to explain data structure, but:

three categories,

Mackerel
Herring,
Blue whiting

Each of these has several years, and each year has several stations

Each station is one row, and contains the data about what one species has eaten at the station.

It ends up looking like this, with space indicating a new column

species year stationID prey1 prey2 prey3 preyetc.
mackerell 2010 52 0.2mg 5.0mg 25.5mg etc..
macekrell 2010 53 etc etc etc
mackerell 2010 54 etc etc etc..
herring 2010 52 etc etc etc


texts: QGIS 3: QgsAuxiliaryLayer with virtual layers

texts:
I'm using QGIS 3.10.1. I've created virtual layers from an ODBC data source and using the (relatively) new feature "label editing". I was able to show/hide/move labels of a layer. But just after saving and refreshing/reloading the layer all 'auxiliary_storage_' labeling attribute values are gone. The same happens if I restart the QGIS project. Saving the project as qgz (compressed, including QGD) doesn't provide any remedy.

The workaround would be to add a custom Python action for/at the virtual layer scope, pointing to an expression like this (it will temporarily re-populate all stored feature attributes from the joined, in-memory auxiliary layer to the target, non-editable vector layer), credits go to @Detlev:

# -*- coding: utf-8 -*-

from qgis.core import *
from qgis.gui import *
from qgis.utils import *

@qgsfunction(args='auto', group='My')
def update_aux(feature, parent, context):
    """
    Updates the auxiliary layer's attributes (if any) of the selected layer

Example usage:

update_aux()

    """
    layer = QgsProject.instance().mapLayer(context.variable("layer_id"))
    auxiliaryLayer = layer.auxiliaryLayer()
    if auxiliaryLayer is not None:
        joinField = auxiliaryLayer.joinInfo().joinFieldName()
        joinKey = auxiliaryLayer.joinInfo().targetFieldName()
        QgsMessageLog.logMessage("Updating auxiliary layer({}) values, join:  ({}->{})".format(auxiliaryLayer.id(), joinField, joinKey), "My", level=Qgis.Info)
        fnm = auxiliaryLayer.dataProvider().fieldNameMap()
        fim = {index: name for name, index in fnm.items()}
        for auxiliaryFeature in auxiliaryLayer.getFeatures():
            filterExpression = "{} = {}".format(joinKey, auxiliaryFeature[joinField])
            featureRequest = QgsFeatureRequest().setFilterExpression(filterExpression)
            layerFeature = QgsFeature()
            layer.getFeatures(featureRequest).nextFeature(layerFeature)
            if layerFeature is not None:
                for i in fim:
                    propertyDefinition = auxiliaryLayer.propertyDefinitionFromIndex(i)
                    auxiliaryFieldName = QgsAuxiliaryLayer.nameFromProperty(propertyDefinition, True)
                    auxiliaryFieldIndex = layer.fields().indexFromName(auxiliaryFieldName)
                    if auxiliaryFieldIndex >= 0:
                        QgsMessageLog.logMessage("changeAttributeValue({}, {}, {})".format(layerFeature.id(), auxiliaryFieldIndex, auxiliaryFeature[i]), "My", level=Qgis.Info)
                        layer.changeAttributeValue(layerFeature.id(), auxiliaryFieldIndex, auxiliaryFeature[i])
            else:
                QgsMessageLog.logMessage("Couldn't find a feature with {} in the current layer!".format(featureRequest.filterExpression()), "My", level=Qgis.Info)
        if iface.mapCanvas().isCachingEnabled():
            layer.triggerRepaint()
        else:
            iface.mapCanvas().refresh()
    else:
        QgsMessageLog.logMessage("There is no auxiliary layer!", "My", level=Qgis.Warning)

Does somebody have an idea, what's going wrong here? I've spent countless hours formatting the labels using the rule-based engine for printing, it doesn't seem to be a proper solution for map production.


texts:
My goal is to create summary information to eventually present rain gauge totals attributes such as 1 day, 7 day, 10 day, 24 day and YTD totals.

I was wondering how to go about this in ArcMap or ArcGIS Pro so that I have 5 attributes with site_id (see picture) and receive time for that day having the maximum rain gauge accumulation which would equal 11 in sensor_class?

I assume I would have to query out the table for the given time periods and use the max statistic (receive_time) for that day to get the rain accumulation.

I am not quite sure how to query that either! I also would need a ModelBuilder for this so it can be easily replicated again or can be used for future readings.


texts: Since ArcMap Projects data "on the fly" why aren't these data sets aligning?

texts: I have a polygon feature class in NAD_1983_UTM_Zone_13N (Green in map).  I projected it to NAD_1927_StatePlane_New_Mexico_East_FIPS_3001 (Pink in map).  Why are these not aligning?  Does this mean the original definition was wrong?  Since ArcMap projects different CRSs "on the fly", I can't figure out why they aren't perfectly aligned.


texts: QGIS - SAGA Majority Filter creating huge output files

texts:
I am trying to perform a workflow that I have made in QGIS which operates on very large rasters (38244x62124).

Leading up to this, I have been able to run various GDAL tools if I specify that the output image is to be compressed via LZW, and if I don't, then the image balloons from about 200MB to 20GB.

I am hitting the same problem when running the "Majority Filter" tool from SAGA. The image is exploding up to the order of GB, which quickly fills up the drive I have to work on and halts the entire workflow.

I have limited space on the machine that I am using. Is there a way to tell SAGA to make a compressed image as an output, similar to how you can compress the image using GDAL?


texts: Model Builder - Selections Add Up at the End

texts:
I am creating a model so that I can help colleagues visualize the workflow I used on a problem.

My task was find all manholes in a specific area that connect to pipes >= 36". Simple enough, I did this outside of model builder and was able to generate the answer. When I put it in Model Builder and run it, the final solution selects all manholes and all pipes (collectively) in that area. That includes all pipes, regardless of size (1 inch, 6 inch, 12 inch, 18 inch, etc.). In the tools I used the derived data from the steps and used pipes that intersect manholes on a new selection. Not sure what to do!

Gravity Mains --> Select by Attributes * where Diamater >= 36

Pilot Boundaries --> Select by Attributes * where Name = NN05

Manholes --> Select by Location * Where Centroid is in NN05

The last selection is Pipes >= 36" that intersect Manholes that have their centroid in NN05


texts: Query with dimension parameter as a metadata in Mapserver's client

texts:
I've read the documentation to use mapserver as a client but I can't see how to query a layer which get a dimension parameter using a metadata value.
The only way I found is to add the dimension parameter/value in the connection string this way :

CONNECTION      "http://host?DATE=2014-08-20&"

Can't it be done with metadata? Is it possible to add custom metadata that will be included in the url or mapserver only understand predefined ones?


texts: Clip OR Overlay in ArcGIS output feature empty

texts: I use the Clip and Overlay to deal with download data. But they do not work and result in failure. Maybe the data is too large for ArcGIS Desktop to calculate. I do not know where I can check whether the data is OK for ArcGIS Desktop.


texts: What is the default method in raster package R resample?

texts:
What is the default method in the resample command from the raster package in R?
I see that you can chose some methods like bilinear or nearest neighbor, but what is the default? This is important for me so I can stick to one method also when I resample with another software. Is the default method 'bilinear' because it is used in the example in the package documentation?

https://www.rdocumentation.org/packages/raster/versions/2.6-7/topics/resample


texts: How to add all layers to OpenLayers from a single GeoServer WMS?

texts:
I want to add all layers to OpenLayers from a single GeoServer WMS. I tried this but it didn't work for me:
<script>
  var layers = [
    new ol.layer.Tile({
      source: new ol.source.OSM()
    }),

    new ol.layer.Image({

      source: new ol.source.ImageWMS({
        url: 'http://localhost/geoserver/wmsExamples/wms?',
        params: {'LAYERS': '*'},
        ratio: 1,
        serverType: 'geoserver'
      })
    })

  ];

</script>


texts: Number a subset of entries based on current sorting

texts:
Is there a way to number a subset of entries in the attribute table based on the current sorting in Qgis?
e.g. select all entries with attribute X, sort them by abttribute Y, create a new column and fill the column (only for the selected and sorted values) whith 1 to n; not rownum or id.


texts: Using arcpy.da.editor on a transactional version raising errors

texts:
I've run up against a wall trying to get an edit session to work on a version. These are the scenarios Ive gone through in the code and their error mesages:

Setting up an edit session not as a context manager: 'workspace already in transaction mode'

target = r'D:\GIS-DATA\bay-arcgis-pro_ParcelEditing.sde\sde.SDE.Staging\sde.SDE.SPRCodeEnforcementCasesStaging'
try:
    edit =  arcpy.da.Editor(r'D:\GIS-DATA\bay-arcgis-pro_ParcelEditing.sde')
    edit.startEditing(False,False)
    edit.startOperation()
    with arcpy.da.UpdateCursor(target,'*') as cursor:
        for row in cursor:
            row[3] = 'test'
            cursor.updateRow(row)
finally:
    edit.stopOperation()
    edit.stopEditing(True)

Traceback (most recent call last):
File "D:\GIS-DATA\editortest.py", line 9, in <module>
for row in cursor:
RuntimeError: workspace already in transaction mode

Setting up an edit session as a context manager: 'The requested operation is invalid on a closed state'

target = r'D:\GIS-DATA\bay-arcgis-pro_ParcelEditing.sde\sde.SDE.Staging\sde.SDE.SPRCodeEnforcementCasesStaging'
with arcpy.da.Editor(r'D:\GIS-DATA\bay-arcgis-pro_ParcelEditing.sde'):
    with arcpy.da.UpdateCursor(target,'*') as cursor:
        for row in cursor:
            row[3] = 'test'
            cursor.updateRow(row)

Traceback (most recent call last):
File "D:\GIS-DATA\editortest.py", line 8, in <module>
cursor.updateRow(row)
RuntimeError: The requested operation is invalid on a closed state

The data is versioned on a SQLServer geodatabase.


texts: Selecting features based on similar name

texts:
I am building a script to integrate into a model that iterates through data. I've been able to get two submodels to iterate through the original data and create large numbers of features in file geodatabase. One submodel creates point features the other creates polygons. I would like to select features (one point and one polygon) based on similar names from the file geodatabase. The point features are civic addresses that fall within a catchment area of each facility. The polygon is the 25 min drive time (service area) from each facility. The point ID of the facility has been integrated into the names of the points and the polygons.
I was hoping to select the features from the geodatabase based on the actual name of the feature class which contains the point ID as seen in the screen shot below. Ex: SA207900045_01 polygon the point ID is 207900045. The corresponding point file is nscafPL207900045_01. The point and polygon features do not contain the point ID information. I was hoping something could be done with 'list features' in python then 'select where' the names are the same but I'm having no luck finding this function. I want to be able to loop this selection. So I would select the point and polygon that have similar names then feed them into a tool in model builder then loop back and find another pair of point and polygon features and feed them into the tool.
Is there anyway to do this in python or model builder? I'm new to both. Below is an image of some sample data.
I'm using ArcGIS for desktop 10.2.2


texts: Using cursor to delete rows with  values in table using ArcPy?

texts:
I am trying to delete any record in my table that has a null value in a certain field (the Grid_Code field) and the code I have does not return any errors and seems to run fine, but doesn't actually delete the records with null values.

Here is my code:

with arcpy.da.UpdateCursor(fc, 'Grid_Code') as cursor:
    for row in cursor:
         if (row[0] == '<Null>'):
            cursor.DeleteRow(row)
         #I also tried adding this and got the same result:
         else:
             pass

Any ideas about what I am doing wrong?


texts: postgis ST_Intersect on 2 valid geometries throws TopologyException

texts:
Using Postgis 2.1 on PostgreSQL 9.3.4, I'm getting an error when I try to intersect two valid geometries:

TopologyException: found non-noded intersection...

These geometries are created on the fly and are just 2 of a much larger table so while the usual tricks (a tiny buffer around either geoms) fixes it in this case, they are not practical in production with large geometry sets.

How can I avoid this in a consistent, performant manner?

See here for code snippet with geometries.


texts: How to query some properties of a list of typenames in a GeoServer WFS

texts:
I want to query selected propertynames from a list of Typenames in a GeoServer service.

I could find an example from OGC WFS 1.1.0 Implementation Specification ( PDF document page 80 Example 10) which is not working for me. Here is the example:

Example 10

Query some properties of a list of feature types. In this case, the attributes WKB_GEOM and TILE_ID are fetched for the INWATERA_1M feature type and all the attributes of feature type BUILTUPA_1M are fetched.

http://www.someserver.com/wfs.cgi?
SERVICE=WFS&
VERSION=1.0.0&
REQUEST=GetFeature&
PROPERTY=(INWATERA_1M/WKB_GEOM,INWATERA_1M/TILE_ID)(BUILTUPA_1M/*)&
TYPENAME=INWATERA_1M,BUILTUPA_1M&

When I try to implement this I got the following error:

"Requested property: PropertyName1 is not available"   The possible propertyName values are:[PropertyName1, , , , , ]"

It would be  great if someone give me a working example for GeoServer.


texts: QGIS export DEM to SAGA grid - no coordinate system

texts:
I have created a DEM using survey points from a GNSS system. Works like a charm so far. I can extract elevation lines etc. However, when I export it as a SAGA grid it looses the coordinate system. I am able to open it in SAGA and run every module and terrain analyzing tool, but the grid is missing the coordinate system.

What am I doing wrong? When exporting to grid I choose the coordinate system that I want it to be in.

QGIS 3.8.0 Zanzibar
SAGA GIS 7.4.0


texts:
I created a Mosaic Dataset (Raster Mosaic) using ArcCatalog 10.3.1 and then added a number of existing Raster Mosaics to that Mosaic Dataset. I then edited the ZOrder of the constituent Raster Mosaics.

When I view my Mosaic Dataset in ArcMap the ZOrder does not appear to be having the desired effect.  My smaller Raster Mosaics are visible when I zoom right in (large scale) but zooming out a bit (smaller scale) causes the smaller Raster Mosaics to disappear as they are covered by my large 'basemap' Mosaic.  It's my understanding that this shouldn't happen as the smaller Raster Mosaics have a smaller ZOrder than the large basemap Raster Mosaic.

How do I get my smaller constituent Raster Mosaics to display at all scales without being covered up?

tl;dr after creating a Mosaic Dataset, adding rasters, and then editing the ZOrder, what is the next step to ensuring the raster display order actually follows the ZOrder?


texts: How to supress 'Start With' Menu in ArcPad 10.2?

texts:
I am fairly new to ArcPad in general and I am attempting to suppress the start with menu as seen in the attached photo. I have a script that is called in the event: onstartup, that loads the .APM I need. I find this menu unnecessary and have been looking for a solution via the APX in ArcPad Studio or via an object/keyword for that particular menu. I have searched geoNet, Stackexchange, and the ArcPad Dev Help but have come up short.

Has anyone had experience in suppressing this menu? or could point me in a more appropriate direction?


texts: Issues with raster data on postgresql / error:Attempting to get pixel value with out of range raster coordinates

texts:
So here is the issue, i am trying to extract values from a big raster divided in 200 or so tiles and applying them to a selection. The points are startpoints, endpoints and midpoints from vector items, lines created using flow tools on SAGA. Everytime i try to run the query it returns me the following error:

  NOTICE:  Attempting to get pixel value with out of range raster
  coordinates:  CONTEXT:  fonction PL/pgsql st_value(raster,integer,geometry,boolean)

Here is my function so far, tried to use a "where" query in order to force him to read only the raster data where it intersects with my lines but that didn't fix a thing and the query failed.

SELECT

    tronc.id AS id,
    tronc.geom AS geom,
    ST_value(raster.rast, ST_startpoint(tronc.geom), TRUE) AS alt_startpoint,
    ST_value(raster.rast, ST_endpoint(tronc.geom), TRUE) AS alt_endpoint,
    ST_value(raster.rast, ST_line_interpolate_point(tronc.geom, 0.5), TRUE) AS mid_point,
    (ST_value(raster.rast, ST_startpoint(tronc.geom), TRUE) + ST_value(raster.rast, ST_endpoint(tronc.geom), TRUE) + ST_value(raster.rast, ST_line_interpolate_point(tronc.geom, 0.5), TRUE))/3 AS mean

FROM

    sinuosity.channels as tronc,
    rasters.tr_rge_alti_5m_2016 AS raster

WHERE

    ST_intesects(tronc.geom, raster.rast, 1)


texts:
I have successfully been using the trace upstream in ArcGIS with 10,000+ points as flags to find what polygons are upstream (using a select of the polygons than intersect the upstream trace).

Now I need to know how far upstream a polygon is. That is, I know which polygons are upstream but I do not know the distance from the flag along the network to the polygon. If I cannot get distances even sequence would do (this is the first poly upstream, this is the second and so on) although distance upstream is preferred.

ArcGIS Geometric solution.
I have close to 100,0000 flags (points) in total.
Trace upstream works great but now I need to add the distance upstream to each feature.


texts: Exporting ArcGIS Metadata as PDF using ArcPy?

texts:
I've been asked to automate the exporting of metadata from ESRI feature classes into .pdf files. We don't want the .xml files, rather a .pdf version of what you see in the Description panel in Catalog. Manually we can achieve this by using the print button in the description panel tools and then using the cutepdf print drive. But I don't see a print option for metadata in ArcPy, all I see are options for exporting to .xml files.

In ArcObjects the IGxView object which looks like it could do what I need.

Is there a way to do what I want in ArcPy or should I use ArcObjects and its IGxView instead?


texts:
I produce printed maps out from OSM-Data. While I work with the data, I find some missing places, names or something to fix and I would like to change/add it to OSM directly. But for several reasons I have to do my work with a local OSM-Copy (osm2pgsql).

My work flow:
I download the OSM-Area of interest and import it into a local PostGIS database. Then I do some SQL queries to fit my needs. I style and edit with QGIS. If I have to correct something, I do it on my local DB. At last the result will be printed...

I'm looking for a new work flow, so that my changes and adds are uploaded to OSM, so that anyone else can use it!

I could do my changes directly to OSM and wait at least one week to download the updated data. For some large-scale projects I do that already, but my small scale works have to be just in time - I often can't wait for the next OSM extract.


texts: How do you export a .cou file without opening the Attribute Table?

texts:
We sometimes receive .cou files which are essentially just .csv files that contain address data, UPRNs etc. for processing.
As I have access to QGIS in my organization I was wondering how to produce these files myself but I couldn't find anything when looking up .cou files.
If possible I'm looking for a way to export this file without having to open up the Attribute Table as it contains around 70,000 rows of addresses


texts: Choosing number of points in Create Fishnet using ModelBuilder?

texts:
Is there a way of getting the analysis "create fishnet" to generate a grid of x number of points in ModelBuilder?

For example, I wish to generate a grid consisting of for example 80 points within a polygon (see picture).


texts: Downloading SoilGrids data globally using R

texts:
I'm trying to download the ISRIC SoilGrids data via R (Version 3.6.2 on macOS Catalina Version 10.15.5) as described in the tutorial here:
https://www.isric.org/explore/soilgrids/soilgrids-access
Following this guide I can download the data for the small example bounding box. However, I'm interested in the global data so my modified R script looks like this:
library(rgdal)
library(gdalUtils)

bb=c(-20037500,-6729000,20037500,8600750) #bounding box parameters (in Homolosine) for whole globe, as suggested on the website
igh='+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs' #proj string for Homolosine projection

gdal_translate(
  '/vsicurl?max_retry=3&retry_delay=1&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/ocs/ocs_0-30cm_mean.vrt',
  "./crop_roi_igh_r.vrt",
  of="VRT",
  tr=c(25000,25000),
  projwin=bb,
  projwin_srs =igh,
  verbose=TRUE
)

gdalwarp(
  "./crop_roi_igh_r.vrt",
  "./crop_roi_ll_r.vrt",
  s_src=igh,
  t_srs="EPSG:4326",
  of="VRT"
)

gdal_translate(
  "./crop_roi_ll_r.vrt",
  "./crop_roi_ll_r.tif",
  co=c("TILED=YES","COMPRESS=DEFLATE","PREDICTOR=2","BIGTIFF=YES")
)

Running
gdal_translate(
  '/vsicurl?max_retry=3&retry_delay=1&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/ocs/ocs_0-30cm_mean.vrt',
  "./crop_roi_igh_r.vrt",
  of="VRT",
  tr=c(25000,25000),
  projwin=bb,
  projwin_srs =igh,
  verbose=TRUE
)

aborts with the following error message:

Checking gdal_installation...
GDAL version 2.4.2
GDAL command being used: "/usr/local/Cellar/gdal/2.4.2_3/bin/gdal_translate" -tr 2500 2500 -projwin -20037500 -6729000 20037500 8600750 -of "GTiff" -projwin_srs "+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs" "/vsicurl?max_retry=3&retry_delay=1&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/ocs/ocs_0-30cm_mean.vrt" "./crop_roi_igh_r.tif"
ERROR 1: -projwin_srs ignored since coordinate transformation failed.
Input file size is 159246, 58034

What I tried so far:

incrementally increasing the bounding box from bb=c(-337500.000,1242500.000,152500.000,527500.000)
up to bb=c(-19861874.000, 1242500.000, 19861874.000, -1500000.000)

This spans a belt from -179.xx to 179.xx Longitude but I can't increase Latitude parameters even to 20°N/S before encountering the error.

Going for a different parameter (mean clay content, which I'm actually more interested in) via

gdal_translate(
  '/vsicurl?max_retry=3&retry_delay=1&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/clay/clay_0-5cm_mean.vrt',
  "./crop_clay_igh_r.vrt",
  of="VRT",
  tr=c(25000,25000),
  projwin=bb,
  projwin_srs =igh,
  verbose=TRUE
)

This gives the error

ERROR 1: Error: Computed -srcwin 71784 36135 16030 -6132 has negative width and/or height.
Input file size is 159246, 58034


texts: QGIS shows incorrect $area in label, but correct in field calculator

texts:
I tried showing $area of polygons in QGIS, and it always showed planimetric square degrees instead of metric units. No matter what I set in project general settings, the value remained the same.

QGIS version is 3.4.3, layer CRS is 4326, and I tried different CRS in the project, which made no difference.

I calculated a new field with the formula: $area / 10000 (to get hectares, which is more readable). And looking it up in Wikipedia, the city area is close to that and is rather correct.

So, you can get $area but have to use field calculator. Is there a workaround to make $<function> work correctly in labels?

Following suggestions by @MAP: the derived tab in object properties is always correct...

...But the reason I'm asking is that the function description promises to work this way too:

...while simple area(geometry) is always planimetric.


texts: Need information about the Envisat satellite

texts: I am new to remote sensing and I want to know about sensor design of Envisat satellite. Can anybody help me or preferably direct me to some good references?


texts: QGIS - Labels on print layout move or don't appear on exported PDF

texts: QGIS 3.16.1 on Windows 10 Pro.  I often manage to get a print layout with label placement just the way I want it, but when I export to PDF, some labels either don't appear or appear in a less desirable location.  Is there some way to control or minimise this effect?


texts: Features disappear from Leaflet after zooming in

texts:
I have a Mapbox studio template which already looks the way I want it to, but when I put it in a Leaflet app, after I zoom in past level 17, a lot of the features disappear.  Zooming in inside Mapbox studio, everything looks fine. It's only in the Mapbox app. Another peculiar thing is that generating the static image also makes the features disappear.

All the layers have visibility from zoom 0-22, so I don't know what the problem is.

Any thoughts?


texts: Accuracy of shape area measurements when GCS shapefile is source of layer in PCS Data Frame?

texts:
I would need to deliver a result shapefile in WGS 84 (GCS: Geographic –lat/lon- coordinate system), the shapefile is originally in WGS84; but in the process of creating this result shapefile I need to measure areas (the ArcGIS Calculate geometry tool doesn’t allow calculating areas with the Data frame coordinate reference system set to a GCS).

The solution I find in ArcGIS is: set the Data Frame coordinate system to a proper projected coordinate system for that zone, work with the shape layer, in WGS 84 (GCS) whereas the Data frame is in a Projected coordinate system, and once finished working and obtained the result shape layer, Export the shapefile by doing Data > Export Data and , for the option “Use the same coordinate system as:”, choose “this layer’s source data”.

Is it possible to make an accurate area measurement in a layer in meters, hectares, acres, etc, when the layer’s coordinate reference system is a GCS (though the Data frame has a projected coordinate reference system defined)?


texts: How to fix The procedure entry point squlites3_open_v2 could not be located in the dynamic link library squlite3.dll?

texts:
Possible Duplicate:
  How to fix QGIS error “Entry Point could not be located”?

I downloaded QGIS to my laptop (Window 7 64b) yesterday but when I click on the desktop logo I get the following warning (!) massage: " The procedure entry point squlites3_open_v2 could not be located in the  dynamic link library squlite3.dll"

I already uninstalled the program once and reinstalled it but the problem persists


texts: Pairing more than two lists using Google Earth Engine

texts:
I would like to know how to pair more than two lists.
I know this is easy to do using .zip function, but how do I pair a third or more lists?
var list1 = ee.List(['20190101T160509_20190101T160509_T17RNJ',
                     '20190104T160649_20190104T161509_T17RLH',
                     '20190104T160649_20190104T161509_T17RMJ',
                     '20190104T160649_20190104T161509_T17RMK']);
var list2 = ee.List(['17RNJ','17RLH','17RMJ','17RMK']);
var list3 = ee.List([1,2,3,4]);

var paired = list1.zip(list2);
print(paired);


texts: Cesium terrain tiles too high by roughly 30m

texts: I'm using the NED 1/3 arc second data set to generate tiles to zoom level 15 for Cape Canaveral in Florida, and finding that the terrain looks accurate in appearance, but the elevation is too high by roughly 30m. The Cesium AGI STK World Terrain set uses NED data as well, and I can see the exact same artifacts and characteristics in the terrain I'm generating as that in the AGI data, except for the elevation difference, so it looks like we're working from the same data.  Is there a parameter or step I'm missing when generating terrain tiles?


texts: Add a vector GeoJSON layer in OpenLayers 5.1.3

texts:
I am trying to add a simple capitals GeoJSON vector layer to my map for basics learning, without managing to get it to work. I spent hours of roaming the documentation API. As follow is my script:

new ol.layer.Vector({
        source: new ol.source.Vector({
            url:'res/world_countries.geojson',
            format: new ol.format.GeoJSON({
                featureProjection: 'EPSG:3857'
            }),
            }),
        })


texts: R coordinates on line nearest to point - multiple points

texts:
In R, I have a collection of points (~15,000) within a polygon.

For every point, I want to calculate the location on the polygon boundary nearest to the point's location. I'm not bothered about the point's distance from the boundary (i.e. gDistance etc.).

I have a method that works (below). I used nearestPointOnLine from maptools but this only works for 1 point. I wrapped it up in apply but with my spdf being very large, this is taking a while;

require(maptools)
require(sp)
require(rgeos)

# create a dummy polygon: coordinates, promote to Polygons and SpatialPolygons
coordsPol = cbind(c(1,2,3,4,4),c(3,2,2,1,3))
verts <- rbind(coordsPol, coordsPol[1, ])
pol <- SpatialPolygons(list(Polygons(list(Polygon(verts)), ID="a")))

# create spatialpointsdataframe (they all fall inside the polygon)
p = data.frame(x=c(1.5,2.5,3,3.5,3.8),y=c(2.75,2.25,2.55,2.9,1.6),ID=c("a","b","c","d","e"))
coordinates(p) <- ~x+y

# Extract the point on the polygon boundary that is nearest to each spdf point.
# The polygon needs to be a matrix (ggplot::fortify) and maptools::nearestPointOnLine
# only works 1 point at a time so requires apply

v2.x <- apply(p@coords,1,function(x) nearestPointOnLine(as.matrix(fortify(pol)[,c(1,2)]),x))[1,]
v2.y <- apply(p@coords,1,function(x) nearestPointOnLine(as.matrix(fortify(pol)[,c(1,2)]),x))[2,]

# The output from nearestpointOnLine can be converted to a SpatialPoints object
v2m <- matrix(cbind(v2.x,v2.y),ncol=2)
v2p <- SpatialPoints(v2m)

Are there any functions (i thought rgdal would have one) that take a polygon and many points and give resulting nearest location points on a boundary (instead of apply and nearestPointOnLine)?


texts:
I like to summarize of the ID-Field with the ModelBulider because in the end I have to build a tool. The problem I have is that the field “GRIDCODE” should get in the new table comma separate. So that it is possible to see the combinations, for example the ID-Field 244 has the GRIDCODE 7, 10, 10. The Field area should get added.

Background: The Field ID-Field is from a fishnet lying over a shape with a classification “GRIDCODE” and got intersect. Now I need just a table with the ID form the fishnetcell and the combinations of the classification-code and the area.
 Does anyone have an idea how to do this with the Modelbuilder?
Thanks a lot.

Here is the last part of the Model so far.

The final table should look like this; in the field area should be the sum area from the combination.


texts: How to merge attribute table properties from one layer to another in QGIS?

texts:
I am using QGIS 2.2.3 (Lyon). I have a vector layer with a linked attribute table. In order to do multiple tests on the layer, i created a copy in a new layer. The attribute table on this new layer is empty. I want this new layer to have the same layout/properties than the original.

Is there a way to merge the layout from one layer to another ?


texts: Replace multiple fields in a column, not characters but words, in R

texts:
I have a dataframe with some thousands of entries each, I've different province names. I need to substitute these names (which are repeated) for others.

Quick example, I have this column (I don't know how to show it like a column here):

**Provincia**   Álava Álava Albacete Albacete Albacete ... Zaragoza Zaragoza

And I need to change for:

**Provincia**
ALAVA_1
ALAVA_1
ALBACETE_1
ALBACETE_1
ALBACETE_1
...
ZGZ_1
ZGZ_1

I've tried the following:

old_names <- unique (df$Provincias)

  Álava Albacete ... Zaragoza

new_names <- c ("ALAVA_1", "ALBACETE_1",... "ZGZ_1")

df$Provincias <- gsub (old_names, new_names, df$Provincias)

  ALAVA_1 ALAVA_1 Albacete Albacete Albacete ... Zaragoza

  1: In gsub (new_names, old_names, df$Provincias) :
    argument 'pattern' has length > 1 and only the first element will be used

  2: In gsub (new_names, old_names, df$Provincias) :
    argument 'replacement' has length > 1 and only the first element will be used

So it only changes the first term (Álava for ALAVA_1) all the times, but it doesn't change the rest of Provincias

Of course I know that I can do this one by one, for the 52 provinces, but then I have more dataframes to do the same... I'm sure there's a better way.


texts: Get specific interaction from the map in OpenLayers 5

texts:
I am developing a map with OpenLayers v5.2.0, and I am having problems for the interaction change.

Currently, I have a function that adds all the controls to the map (selection, draw, etc). Which when it changes to a specific interaction eliminates the other interactions of the map.

The main problem is that I need several maps in the project, but each map has specific interactions. So, I need to separate the assignment of interactions separately, so that I can generate a map with the selection and drawing of circles (for example).

How can I distinguish the interactions of the map when doing map.getInteractions ()? To know which interaction to deactivate and which does not.

`addControlButtons : (map, vectorSource) ->
draw = undefined
snap = undefined
select = undefined
modify = undefined
window.app = {}
app = window.app

#SELECT
app.SelectControl = (button) ->
  changeSelectInteraction = ->
    addSelectInteraction()

  button.addEventListener 'click', changeSelectInteraction, false
  button.addEventListener 'touchstart', changeSelectInteraction, false

  element = OpenLayers.createElementDiv('select-control', button)
  ol.control.Control.call this, element: element
  return

ol.inherits app.SelectControl, ol.control.Control
selectControl = new app.SelectControl(OpenLayers.createControlButton('S', 'Selección'))
selectControl.setMap(map)

#DRAWPOINT
app.DrawPointControl = (button) ->
  changeDrawPointInteraction = ->
    addDrawInteraction("Point")

  button.addEventListener 'click', changeDrawPointInteraction, false
  button.addEventListener 'touchstart', changeDrawPointInteraction, false

  element = OpenLayers.createElementDiv('draw-point-control', button)
  ol.control.Control.call this, element: element

ol.inherits app.DrawPointControl, ol.control.Control
drawPointControl = new app.DrawPointControl(OpenLayers.createControlButton('P', 'Dibujar Punto'))
drawPointControl.setMap(map)

#DRAWPOLYGON
app.DrawPolygonControl = (button) ->
  changeDrawPolygonInteraction = ->
    addDrawInteraction("Polygon")

  button.addEventListener 'click', changeDrawPolygonInteraction, false
  button.addEventListener 'touchstart', changeDrawPolygonInteraction, false

  element = OpenLayers.createElementDiv('draw-polygon-control', button)

  ol.control.Control.call this, element: element

ol.inherits app.DrawPolygonControl, ol.control.Control
drawPolygonControl = new app.DrawPolygonControl(OpenLayers.createControlButton('D', 'Dibujar Polígono'))
drawPolygonControl.setMap(map)

addDrawInteraction = (drawType) ->
  map.removeInteraction select
  removeDrawInteractions()
  modify = new ol.interaction.Modify (source: vectorSource)
  map.addInteraction modify

  draw = new ol.interaction.Draw(
    source: vectorSource
    type: drawType)
  map.addInteraction draw

  snap = new ol.interaction.Snap (source: vectorSource)
  map.addInteraction snap

  draw.on 'drawstart', (e) ->
    vectorSource.clear()

removeDrawInteractions = ->
  map.removeInteraction draw
  map.removeInteraction snap
  map.removeInteraction modify
  vectorSource.clear()

addSelectInteraction = ->
  removeDrawInteractions()
  select = new ol.interaction.Select
  map.addInteraction select`

This is my function.


texts: How can I omit points when filtering OSM data with Osmosis

texts:
I am trying to get data for rivers/waterways around the world. I have downloaded the planet file from OSM (as well as a North America file from Geofabrik), and am trying to use Osmosis to filter out only the ways tagged as rivers, streams, etc. The command I am using is the following:
osmosis --read-pbf planet-latest.osm.pbf --log-progress --tf accept-ways waterway=river,stream,canal natural=water --used-node --write-pbf planet-latest-waterways.osm.pbf
When I run this on the North America file, it completes and results in a 2GB file. However, when loading the file in QGIS, the points take quite a while to process and if I choose not to display them, QGIS runs for a long time before displaying any lines (presumably silently iterating through and rejecting all the points). Many of the points are irrelevant (e.g. bus stops). Is there a way to not include the points using Osmosis, so that QGIS jumps straight to rendering the lines and polygons?
If I am going around this in the wrong way, please let me know! I am happy to learn new tools, I just am not sure which ones are the best for this type of job.


texts: Visualizing street shape data for London?

texts:
I am trying to visualise some street-level data I got, but since I'm very new to GIS I don't know the correct terminology to Google. My data looks like this:

Street, Incident Count
Congo Road, 12
Roydene Rd, 21
....

Using a combination of geocoding and Leaflet.js, I managed to get something like this:

However, with geocoding I only get the geometric center of each street, so the visualisation isn't very informative. I'd like to get an image, where each street is colour-coded according to how many incidents occurred in it, e.g (excuse my poor Paint skills)

I was thinking I can use OSM to get the coordinates of each node that a street consists of, so I'd be able to represent a street as a collection of points, which would make the visualisation a bit better.

Any other suggestions?


texts: PYQGIS tool to get Raster XYZ data from inside Polygon?

texts: I am working on automating a process to overlay a polygon on a raster layer and get the data inside polygon for calculations. Is there a method to extract the XYZ data of the raster that is inside the polygon?


texts:
In ogr2ogr, I can successfully read *.DBF files and write them to PostgreSQL DB:

SET DBCONNECT=host=bla port=blue dbname=blub user=me password=easy active_schema=now
SET PGCLIENTENCODING=UTF8
SET SOURCE=\\fileserver\test\
ogr2ogr -f "PostgreSQL" PG:"%DBCONNECT%" %SOURCE%DBF\countries.DBF  --config SHAPE_ENCODING "CP850" -overwrite
ogr2ogr -f "PostgreSQL" PG:"%DBCONNECT%" %SOURCE%DBF\taxes.DBF  --config SHAPE_ENCODING "CP850" -overwrite
SET PGCLIENTENCODING=LATIN1

However, I'm trying to find a simple solution to create (and update) an additional table in the same schema that logs the following values:

file found or not
Data directory
username of import
Date of import
modification date of the specific file

As I am pretty new to ogr2ogr.


texts:
From MS Access, I'd like to create a link to the attribute table of an SDE feature class in a spatial database.  For example, I have the following connection string which works for an OLE DB connection:

ODBC;DRIVER=SQL Server;SERVER=WCARCSDE00;Trusted_Connection=Yes;DATABASE=MyDb

I want the equivalent connection string for a Spatial Database Connection.  Something like the following (which does not actually do what I want because the SERVICE parameter is ignored):

ODBC;DRIVER=SQL Server;SERVER=WCARCSDE00;SERVICE=sde:sqlserver:wcarcsde00;Trusted_Connection=Yes;DATABASE=MyDb

Is this possible?  Do I have some fundamental misunderstanding about how this all works?


texts: Cannot Open VRT of Sentinel2 image created by AutoCreateWarpedVRT

texts:
I am trying to create a warped VirtualRaster of a Sentinel2 image. Here the Python code I tried:

from osgeo import gdal, osr

srs = osr.SpatialReference()
srs.SetFromUserInput('EPSG:32632')

gdal.UseExceptions()
d = gdal.Open('SENTINEL2_L1C:S2A_OPER_MTD_SAFL1C_PDMC_20160506T205836_R122_V20160506T100527_20160506T100527.xml:10m:EPSG_32632')
vrt = gdal.AutoCreateWarpedVRT(d, None, srs.ExportToWkt())
gdal.GetDriverByName('VRT').CreateCopy('test.vrt', vrt)

The AutoCreateWarpedVRT command actually succeedes, i.e. the returned dataset can be effectively used to read data and metadata.

The problem arises with the last command, which raises the following exception:

RuntimeError: No such file or directory

Please notice that:

The test.vrt file is actually created and filled, but gdal.Open raises the very same error as above (I suspect that the error raised by CreateCopy is actually raised by an internal gdal.Open)
If I inject some error in the XML, gdal.Openraises a different exception, which is actually related to the injected error.
gdalinfo, gdal_translate and QGIS can successfully read test.vrt.
The procedure does work as expected if, instead of a SENTINEL2 image, I use LANDSAT8 (i.e. plain GeoTIFF) data.

I am using Ubuntu 14.04 with GDAL-2.1.0 compiled by myself with the flags --with-dods-root=/usr --with-ogdi=no (of course, gdalinfo, gdal_translate and QGIS are also compiled and ran with the same build).


texts: Select by Location in ArcGIS with numerous polygons?

texts:
Please help me with this question about "Select by Location". I'm currently using ArcGis 9.3.1

I have multiple "Selector" polygons. Ideally, I would like to select a points shapefile with "Select by Location" using the selector polygon. Then I create a column in the attribute table and put a value of "1" in only the rows that are selected.

So the process is farily straightforward, except that I have about 30 selector polygons, so it's hard to repeat the process of select by location, create column, tag the rows that are selected, over and over.

I've tried to think of a few ways of doing this, but they don't seem to accomplish exactly what I might need. I could try a model in Modelbuilder, but my polygons are constantly changing, and I'm not sure of a way for Modelbuilder to grab all the latest polygons. Other way I thought was merging all the selector polygons into one shapefile, then maybe doing a spatial join? However I accomplish this, I really need separate columns for each selector polygon in the table, so the points file can show any overlapping polygons.


texts: Should Legend entries be singular or plural?

texts:
In a map legend should the entries be singular or plural?  What is correct?  Road or Roads?
A Google image search for "map legend" lead me to conclude there is no consistent way of having a legend entry be singular or plural. The more authoritative maps I've seen have legend entries in singular. (USGS, National Geographic, etc)
So my question is, in what cases is a plural legend entry correct?
This guide says singular:
http://gisgeography.com/map-elements-how-to-guide-map-making/
entry 32: Features in Legend are Present and Singular
Example of a legend with singular, and plural entries that I made.


texts: Uploading an asset into asset folder Earth Engine Python?

texts:
Would like to know if there is a way to upload a TIF as an Earth Engine asset, but have it ingested and stored in a specific folder in my Earth Engine assets via Python. I have a bunch of images to upload and would like to batch upload them into a specific folder in Earth Engine assets rather than upload all images to the root directory.

I was originally using the earthengine command line code:

earthengine upload image --asset_id=users/username/asset_id gs://bucket/image.tif

To copy images from Cloud Storage to Earth Engine and it works fine. But again, I would like to upload to a specific folder in my Earth Engine assets.


texts: I can see the source code of the core QGIS plugins from the repository in github. But is there any repository of source codes for QGIS toolbars also?


texts: How to link up my GeoServer layers with my standalone geowebcache instance

texts: How do I get my GeoServer layers show up on the demo page of gwc. Do I have to edit the config XML manually for each layer?


texts: I need to calculate shortest distance from an inputted point to the nearest lines. Nearest 4 lines

texts:
I have a huge dataset of points and roads. The point file is of Residential houses or Comercial shops/ markets. And the Line network is the roads around them. Ideally, I want to take 1 point (input) and find out its distance from nearest 4 lines (roads). The roads shapefile also includes some attributes (including name and Width). The final output should include Distance(meters) from inputted point to the nearest 4 roads and the attributes of these 4 roads.
roads = gpd.read_file (r'fRoads.shp')
long = float(input("Enter longitude E value ))
lat = float(input("Enter latitude N Value ))
radius = float(input("Enter Radius value in km: "))
p1 = Point((long, lat))
df = pd.DataFrame({'a':[lat,long]})
po = gpd.GeoDataFrame(geometry = [p1], crs = roads.crs)
bufo = createbuf(po,radius)
cliped = createClip(bufo,roads)
fig, ax= plt.subplots(figsize = (15,7))
roads.plot(ax=ax, color = 'black')
bufo.plot(ax=ax, color = 'None', edgecolor = 'red')
po.plot(ax=ax, color = 'blue')

fig, ax= plt.subplots(figsize = (15,7))
bufo.plot (ax = ax, color ='None', edgecolor = 'Black')
cliped. plot (ax = ax, color = 'red')
po.plot (ax = ax, color = 'Black')
print (cliped.length)
dist  = po.distance(Point (roads))
print (dist)

I am using the Distance Function from shapely. but it gives the distance between 2 points not between a point and line.

I need something equivalent to the following code but without Arcpy and that should be Open Source.
    import arcpy
    from openpyxl import load_workbook

    arcpy.env.workspace = r'C:\Users\MJ\Desktop\populationFunc\roads func\froads2'
    arcpy.env.overwriteOutput = True

    long = float(input("Enter longitude N value e.g: 77.xxxxx: "))
    lat = float(input("Enter latitude E Value e.g: 28.xxxxx:  "))
    radius = float(input("Enter Radius value in km: "))
    spatial_reference = arcpy.SpatialReference(32642)
    p1 = arcpy.Point(long, lat)
    pnt_geometry = arcpy.PointGeometry(p1, spatial_reference)
    roads = 'fRoads.shp'
    point = r'C:\Users\MJ\Desktop\populationFunc\roads func\propData\propPoints.shp'
    search_radius = '1500 Meters'
    location = 'NO_LOCATION'
    angle = 'NO_ANGLE'
    out_table = "out_table.dbf"
    closest = 'ALL'
    closest_count = 4
    arcpy.GenerateNearTable_analysis(pnt_geometry, roads, out_table, search_radius, location, angle, closest, closest_count)
    fp = arcpy.SearchCursor(out_table)
    for f in fp:
        print f.getValue("NEAR_FID")
        print f.getValue("NEAR_DIST")


texts: Intersecting raster layer with grid to find average within each cell in QGIS

texts: I have a raster of population density and would like to intersect it with a polygon layer and have each polygon carry the average for the underlying raster cells. How would I go about this? I feel like they need to be intersected but what else?


texts: Merge Two Polygons Using Editor merge in arcGis

texts:
i want to create a ArcEngine  function using .net that allow to do the same thing like (see picture):

select 2 polygons
click merge

Is there somebody who can help please!!
Sincerely


texts: Expression does not evaluate to column or table name?

texts:
This is my code:

Include "MAPBASIC.DEF"
Include   "MENU.DEF"

Dim s_filename As String
s_filename = FileOpenDlg("","","TAB"," Open Table Interactive")
Browse * From s_filename

I want to open a file from my computer.


texts: Selecting points within polygon AND sum up point layer's attribute

texts:
i have two layers:
1. a point layer. It contains a value.
2. a polygon layer.

I want to sum up the values of the points within each polygon.

I tried to do it with "Count points in polygon(weighted)" - but it's extreml slow. Not even 1 % finished after 10 minutes (and it's only ~ 50000 points!). Is that correct?!


texts: Get attribute table from WMS GeoServer in OpenLayers

texts: I have WMS so I can easily display it in the map, but I want to get the attribute table of the layer from WMS. All solutions I found are about getFeatureInfo which is giving properties when clicking on a place on the map. Is there a way to get the attribute table of a layer using WMS?


texts: Select Layer By Location (Management) gives different result in ArcPy and ArcMap?

texts:
I have two line Shape files. When I use the Select layer by location tool to select two intersecting lines, it works fine and selects only lines which are actually intersecting.

But when I use this in Python Script by below code, it select extra lines which are not actually intersecting the Source layer line.

arcpy.SelectLayerByLocation_management("Split_Line","INTERSECT","Line")

Why this tool works different in Python script?

I am actually trying to make an offset of a given Shape line. For this I created a buffer in memory, converted it into lines, split at vertices and delete the intersecting Feature to original line.

 #specifying buffer distance
    buffer_feaure = "in_memory" + "\\" + "buffer" # temporary feature
    line_lyr = "in_memory" + "\\" "line" # temporary feauture
    #Make buffer along line
    arcpy.Buffer_analysis(InputLine, buffer_feaure, "20 Meters", "LEFT", "FLAT", "ALL")
    #Convert buffer to polyline
    arcpy.PolygonToLine_management(buffer_feaure, line_lyr)
    #Line Layer Split at vertices
    Splitline_lyr = "in_memory" + "\\" "Splitline" # temporary feauture
    arcpy.SplitLine_management(line_lyr,Splitline_lyr )
    # Selecting layer to intersecting to the line inpput
    arcpy.SelectLayerByLocation_management(Splitline_lyr,"INTERSECT",InputLine)
    # Deleting Selected Features
    if int(arcpy.GetCount_management(Splitline_lyr).getOutput(0)) > 0:
        arcpy.DeleteFeatures_management(Splitline_lyr)

It runs until it split lines at vertices and give below output.


texts: QGIS isn't bringing all the features from an WFS URL

texts:
I have a geoserver installed and some layers published in it. There is a specific layer that has 70 features. The problem comes out when I try to load this layer in QGIS 2.14, only 19 features are shown.

If I download the shapefile of this layer, there are all the 70 features and even when I try to load the shapefile in QGIS, all the 70 features are shown.

The problem only occurs when I try to load the WFS URL (instead of 70 features, only 19 are shown).

If I try to load the WMS URL, everything is shown correctly.

May someone give a clue?

I'm using geoserver 2.4.4 and QGIS 2.14.

I've already tried using WFS URL with maxFeatures setting it to 100, and no success.


texts: How to make ArcGIS Pro use a cloned Python project environment when it keeps resetting to default?

texts:
Just getting started using ArcGIS Pro 2.1.0, am following instructions on this page:
https://solutions.arcgis.com/water/help/cctv-manager/get-started/install-python-package.htm
However after I clone the python environment I can see the clone in the "manage environments" window, and I can click the radio button to activate it, I see the "installed packages" window switch to the cloned version, but as soon as I click the Close button to close the Manage Environments window, the Python Package Manager window resets to the default python environment.  If I quit and relaunch ArcGIS Pro like it says to, I still see the default Python environment being used.
Any ideas how to make the clone stay the active environment so I can install packages into it?


texts: Enabling a grid overlay in QGIS Georeferencer

texts: Is there any way on enabling on QGIS Georeferencer a temporary grid overlay? I am having a very hard time georeferencing several map images that only have coordinate crosshairs at the edges of the image, and a grid will greatly help in finding the latitude and longitude intersection.


texts: Storing features in NoSQL database like Cassandra?

texts:
I want to make a new contrib development based on NoSQL databases like Cassandra.

Does anyone know of a University or Research Group, or anyone else, who works in "Geo-Cassandra" currently?


texts:
The docuentation leaves me in the dark about this. We used to call static mapfiles directy by:

 http://www.domain.com/cgi-bin/mapserv?map=/var/www/mapfiles/mf.map?[parameters]

But now I want to call a mapscript instead of the mapfile directly.

How do I do this?

EDIT: i'm using PHP mapscript 4.10, on mapserver 5.03. PHP is not configured to run as CGI.


texts: ogr.driver.Open() not working on second execution

texts:
I have the following code.

                fn = 'C:/temp/file.shp'
                layerfields = QgsFields()
                layerfields.append(QgsField('ID', QVariant.Int))
                layerfields.append(QgsField('ELEV', QVariant.Double))

                writer = QgsVectorFileWriter(fn, 'UTF-8', layerfields, QgsWkbTypes.LineString,
                                         QgsCoordinateReferenceSystem('EPSG:4326'), 'ESRI Shapefile')

                del(writer)

#Open Layer in OGR

                driver = ogr.GetDriverByName('ESRI Shapefile')
                ds = driver.Open(fn, 1)
                srs = osr.SpatialReference()
                srs.ImportFromEPSG(4326)
                Dst_Layer = ds.GetLayer()

                layer_iface = iface.addVectorLayer(fn, dst_layer_name, 'ogr')

The code works perfectly fine when it is executed for the first time. However, if I execute it again, it produces an error:

AttributeError: 'NoneType' object has no attribute 'GetLayer'

The problem seems to be in driver opening the file the second time, when it is already opened. The error goes away when I delete the opened layer (due to layer_iface expression). However, deleting the same layer through the python script prior to the above code block does not fix the error.

Furthermore, I tried to execute this block twice using a for loop. Lets say first execution is execution 1 that has this block loop twice A and B. Then second execution 2 that had this block loop twice A and B.
I expected the code to run perfectly on run 1A but encounter an error on run 1B, 2A and 2B. However, the runs 1A and 1B run smoothly, and the error occurs on 2A, 2B.

I'm not sure what exactly is the problem, but I do think it's got something to do with closing the driver.Open(fn, 1) line.

Could anybody help me out here? TIA


texts: Is it possible to estimate vegetation ( crops) height using free of charge satellites, like sentinel or Landsat, with QGIS?


texts:
I'd like to display the same vector data in different ways as different layers. While it's very handy to be able to copy and paste the layer style in 1.8+, it would be nice to be able to select a layer and right-click to copy or duplicate it directly.

Any easy way of doing this, like a Python command?

Update: @dakcarto's suggestion works perfectly. Since I needed to add many duplicate layers, I adapted it to loop as many times as needed:

for n in range(1,5):
    iface = qgis.utils.iface; vl = iface.activeLayer(); iface.addVectorLayer(vl.source(), vl.name() + "_clone" + str(n), vl.providerType())

The only problem is that the new layer becomes the activeLayer, so it adds "_clone" and the number to the previous name, ending up with example_clone1, example_clone1_clone2, up to the unwieldy example_clone1_clone2_clone3_clone4_clone5. Need to figure out how to duplicate from the originally selected layer, now.


texts: How to shift vertices of shapes in shapefile using FOSS python library?

texts: I need to shift each vertice of shape in a shapefile (of any type:point, line, polygon, multiline etc.) with Python. I need to loop in all shapes, get the vertices, get it's coordinates and for each X,Y update with new values. Which Free and Open Source Software (FOSS) library should I use (ogr, fiona shapely etc.)? Any example?


texts:
CRS is wgs 84 epsg:4326

Still the map and the data on screen is showing tilted.


texts: Finding which shape contains coordinates (over function) using R?

texts:
I'm working in R and I am new to maps and GIS and I have no idea how it works. But I am trying to return which part of the city specific coordinates are located.
I got the shapefile from here:

  https://opendata.arcgis.com/datasets/ac4e22c0d8444b8b81a0ee9af3273f40_2.zip

I read the file and transform it:

library(rgdal)
library(mapview)
nhood <- readOGR(dsn = "hverfahlutar.shp")
nhood2 = spTransform(nhood, CRS("+proj=longlat +datum=WGS84"))

the dataframe I am working with is this:

schools <- structure(list(school = skolar <- c("Vesturbæjarskóli", "Landakotsskóli",
                      "Barnaskóli Hjallastefnunnar í Reykjavík ", "Tjarnarskóli", "Austurbæjarskóli",
                      "Langholtsskóli", "Vogaskóli", "Suðurhlíðarskóli", "Brúarskóli",
                      "Hlíðaskóli"), long = c(64.1497778, 64.1478234, 64.1287022, 64.1406399,
                                              64.1417735, 64.1380849, 64.131696, 64.1216781, 64.1256253, 64.130929
                      ), lat = c(-21.959836, -21.9517121, -21.9311303, -21.9365311,
                                 -21.9237807, -21.8623288, -21.8606524, -21.9040609, -21.9130514,
                                 -21.9122735)), .Names = c("school", "long", "lat"), row.names = c(NA,
                                                                                -10L), class = c("tbl_df", "tbl", "data.frame"))

I can plot it with points so the data seems to be okay:

plot(nhood2)
points(schools$long ~ schools$lat, col = "red", cex = 1)

What to do next I do not know. I think I need to use the over function but it has not been working for me.

What I need to do is, for each coordinates in each row of school. I need to return the corresponding value in nhood2$HEITI (the name of the neighbourhood.)

I think I need to prepare the data in the schools dataframe for this, but I don't understand how... or why.


texts: Calculate daily precipitation data from hourly GSMaP over several years

texts:
I am still new to the Google Earth Engine and try to calculate the daily precipitation based on the hourly GSMaP data over a longer period of time. I know, that there are already a couple of questions about that topic, so I put an eye on this and this. However, as I need to calculate the daily precipitation over several years, I use a list of Dates. It seems to work in general, but the first image of the image collection is always empty and it seems like, the daily sum of the precipitation data for the other images is not calculated properly. When I compare my results with the GSMaP webmap, my image seems not to cover all raining events of that day. Has somebody an idea why the first image has no band and the rest might not be calculated correctly?

//Add Dates
var iniDate = ee.Date('2019-05-01')
var endDate = ee.Date('2020-05-19')
//Add ImageCollection
var filterGSMAP = GSMAP.filterDate(iniDate, endDate).select(['hourlyPrecipRateGC']);
//print(filterGSMAP)

var difdate = endDate.difference(iniDate, 'day') //calculate steps/days in between

var createList = ee.List.sequence(0, difdate) //make list with the needed number of entries

var listdates = createList.map(function(day){ //change list into list of dates
  return iniDate.advance(day, 'day')
})
print(listdates)

//Summarize hourly precipitation data to daily precipitation data
var gsmapImageCollection  = ee.ImageCollection.fromImages(listdates.map(function(summarize_day){
    var year = ee.Date(summarize_day).get('year'); //get the year of the list entry
    var doy = ee.Date(summarize_day).getRelative('day', 'year'); //get the "day of year" of list entry
    var filteredYear = filterGSMAP.filter(ee.Filter.calendarRange({ //filter the Image collection for the year
    start: year,
    field: 'year'
    }));
    var filterDay = filteredYear.filter(ee.Filter.calendarRange({ //filter the Image Collection for DOY
    start: doy,
    field: 'day_of_year'
    }));
    var currentDate = ee.Date(summarize_day)
  return filterDay.sum().copyProperties(filterDay).setMulti({
    Date: currentDate
  });
}))
print(gsmapImageCollection)

Map.addLayer(filterGSMAP.first(), GSMAPvis, 'GSMap');

Here is the Code in GEE


texts:
I recently downloaded QGIS on a Mac OSX 10.7.4 (lion) and purchased a master map shape file, of a particular area.

I am trying to divide it into polygons following features on the master map vector file (.gml 2.shp) and it just won't snap. I have been to snapping and digitising options box, and ticked and changed tolerances far and wide but it just will not snap.

My map layer and topographic options are ticked. I've tried the general options tolerances and still it does not want to snap. I have tried changing vertex to segment and both to no avail, as well as pixels to map.

If I am working at 1:2500 map scale in metres what would a typical tolerance be?  (tried 1- 50)

I think there must be a simple solution to enable snapping or have i down loaded file in the wrong format? In MapInfo one had to use keyboard strokes 'S' for snap and 'T' for trace, but this is not mentioned and you can be sure I have tried!


texts:
My DEM elevation values varies from -255 to 81.

What does the negative values imply?

I have downloaded the dem from BHUVAN.

About 90 % of values fall in negative category.

How can I convert these values to positive?


texts:
I am working with Geoserver for some weeks now and were able to deploy some services already. For my thematic map (e.g. temperature) I created a custom style displaying all tiles in specific colors (blue, yellow, red). Therefore I want to create a legend, but I really dont know how to do this. This legend should show 3 polygons in 3 different colors and some text right of it. Something like "[red polygon] > 30".

custom style sld:

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
    xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
    xmlns="http://www.opengis.net/sld"
    xmlns:ogc="http://www.opengis.net/ogc"
    xmlns:xlink="http://www.w3.org/1999/xlink"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <NamedLayer>
    <Name>Temperature</Name>
    <UserStyle>
      <Title>My Temperature Style</Title>
      <FeatureTypeStyle>
        <Rule>
            <PolygonSymbolizer>
               <Fill>
                 <CssParameter name="fill">
                   <ogc:Function name="Interpolate">
                     <ogc:PropertyName>dTA</ogc:PropertyName>

                     <ogc:Literal>0</ogc:Literal>
                     <ogc:Literal>#0000ff</ogc:Literal>

                     <ogc:Literal>15</ogc:Literal>
                     <ogc:Literal>#ffff00</ogc:Literal>

                     <ogc:Literal>30</ogc:Literal>
                     <ogc:Literal>#ff0000</ogc:Literal>

                     <ogc:Literal>color</ogc:Literal>
                   </ogc:Function>
                 </CssParameter>
               </Fill>
            </PolygonSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
   </NamedLayer>
</StyledLayerDescriptor>

any idea?


texts: Satellite Images for the Mediterranean Sea Oil Spill

texts: I am looking to obtain free satellite images to detect oil spills in the Mediterranean sea. What sensors would be most appropriate for this task and where can I download this type of data?


texts: I would like to add "http://www.arcgis.com/" layer to my maps but I'm unable to connect.  After adding the URL into connection details it gives me an error. "Failed to retrieve service cap":  Not sure what I'm doing wrong.


texts: Seeking cloud-based alternatives to Google Fusion Tables?

texts:
I'm a beginner in web-mapping. I was playing around with Fusion Tables, but I realized that I had too much data to display.

What are some alternatives to Fusion tables?

I have 9 layers of data (one layer is parcel data so there's over 500000 records) and I have to be able to query the data (i.e. how many parcels are within 150m of a school etc.)


texts: ArcGIS JavaScript API: map loads after pan

texts:
I'm using this.mapa.map.centerAndZoom(center, zoom); with 2 different zooms, 4 and 5.

    loadModules([
            'esri/geometry/Point',
        ]).then(([Point]) => {
            switch (evt.target.value) {
                case 'country_cl':
                    center = new Point([-70.71582015815409, -33.47133434146025]);
                    zoom = 4;
                    break;
                case 'country_co':
                    center = new Point([-73.63857136146157, 4.0920688808058]);
                    zoom = 5;
                    break;
                case 'country_mx':
                    center = new Point([-99.94867549215655, 20.55088183550196]);
                    zoom = 5;
                    break;
                }
            this.mapa.map.centerAndZoom(center, zoom);
        });

Chile has zoom 4, Colombia and Mexico has 5.

When I pan from Chile to Colombia or Mexico, the map gets blank, and it get loaded after the pan:

How can I load the map, with the new zoom (4 when comming from 5, and viceversa) before the centerAndZoom?


texts: Trying to utilize the contour extraction QGIS tool but having trouble creating smoothed contours. The smooth, generalize, and dissolve tools have not worked for me. Also the plugin for contours appears to no longer be working?


texts:
I have been learning as much web mapping as possible lately. I really enjoy it, as it seems what you need to know to keep up with the geospatial field. I do have a question that is probably easy for most people to answer.

What is the best way to query/filter/display/symbolize data based on attributes stored on a spatial database and displayed through something like OpenLayers. I know SQL is great for querying data. But how does this fit into the aspect of a web page?? Is this where PHP comes into play? Or Ajax? I am having a lot of trouble desciphering between where the serverside/clientside work ends and begins.

Open layers is great and easy to make maps. But I want to do more with these maps and having trouble what exactly I should learn next in order to actually do something with my data. Can I do all my querying within javascript if something is served through WMS? Can I do basic queries to my GeoJSON files? Or is that where filtering comes in? I saw with CartoDB, they return a url that has the SQL in it. Is this an httpGetRequest? If I have data in MapBox, how can I filter through it and display data based on a attribute, pick from a drop-down box?

Sorry if this is a lot of question. Haha. I am just trying to grasp the concept of how to actually manipulate the data I have displayed on OpenLayers. I understand the differences between the database/map server/APIs and their function. But tying them together in a way which can change how/what/where/when my data displays is whats getting to me. For example, if I want data to be displayed based on an an attribute I type into a box. Where is this work done? Or if I want to find the closest bar to a point I place on the map, what is it on the serverside or clientside that goes through my data and finds this?


texts: How to speed up slow performance on postgresql Postgis?

texts:
I'm using postgresql 9.6 as my database.
It takes too much time than the query itself.
The query take about 160ms which is not as good as the wanted performance.
My query is like :

SELECT DISTINCT code from table where ST_intersects((ST_Transform(ST_PointFromText('POINT(lon lat)', 4326),3857)),geom)

I used this sql query to create spatial indexes

CREATE INDEX table_geom_gist ON table USING GIST (geom);

but it still consume the same time :

VACUUM ANALYZE table

and also tried to update stats :

SELECT UPDATE_GEOMETRY_STATS(table, geom);

but it fires an error that says :

the column table doesn't exist !

Any advice on how to speed up my request and solve this issue ?


texts: How to select polygons in polygon and export to CSV?

texts:
I have two polygon layrers, each layer a table in PostgreSQL:
One layer with polygons areas in a city.
Another layer with buildings (polygons) that are within areas.

How can I to do to know what buildings are in each areas?
Can I have a list or csv with that information?

Something like: layer1(area) -> layer1, 3, 4 (buildings)

Thank you


texts: Iterate Feature Selection in my sub model is not recognizing data in my main model

texts:
I have an iterate feature selection in a sub model; however, when I go to add the input for the iterator, the data that I want to add to it from the main model is not showing as a choice. I have tried closing the windows multiple times, refreshing the models, creating the sub model and dropping it in the main model, etc.

Can anyone help me with this?

So with the sub model, when I go to add the Airfields_Border_Clip_Layer file for the Input Features, I cannot find it in the list. I hope this makes sense.


texts:
I have the following code snippet I run in QGIS 3.4 which goes through a list of shapefiles (my_wild) and runs a polygon to lines process. It performs the process and adds them to the map as expected, but how do I do work on them from there?

I want to perform some attribute work and do further processing on the outputs. The print function returns:

output_f8c71c3d_8781_4057_8295_3dd7099f214d

but then the next poly2line_result.startEditing() function returns:

AttributeError: 'str' object has no attribute 'startEditing'

for f in my_wild:
    #get perimeter line for the study area
    poly2line_result = processing.runAndLoadResults(
        'qgis:polygonstolines',
        {
            'INPUT':f,
            'OUTPUT':'memory:'
        }) ['OUTPUT']
    print(poly2line_result)
    #add attribute field and commit changes
    #poly2line_result.getFeatures()
    poly2line_result.startEditing()
    if poly2line_result.dataProvider().fieldNameIndex("Segments") == -1:
            poly2line_result.dataProvider().addAttributes([QgsField("Segments",QVariant.Double)])
    poly2line_result.commitChanges()


texts:
I am running a docker-compose (postgis+pgadmin) inside an EC2. Got everything up and running. But I got a problem using ogr2ogr.
So, when I have uploaded files (i.e., geojson) LOCALLY in my computer (NOT EC2), I had simply
ogr2ogr -f "PostgreSQL" PG:"host=localhost user=user dbname=db password=password" filename.geojson

However, when using an EC2, putting 'localhost' as the host does not seem to work and I get this error:

ERROR 1: PQconnectdb failed. could not connect to server: Network is
unreachable   Is the server running on host "********" and accepting
TCP/IP connections on port 5432?
ERROR 1: PostgreSQL driver doesn't currently support database
creation. Please create database with the `createdb' command.

Where can I find my host IP address?


texts: ArcGIS Layer Visibility Shortcuts or Tools

texts: I am wondering if there are ways in which we can enable the layer visibility when holding a key and turn the layer off when releasing the key (in ArcMap). Something like in some online games you can check stats by holding TAB.


texts: Creating Vector Map Tiles in ArcGIS Pro with attribution

texts:
I am looking at the ESRI flavour of vector map tiles as a possible replacement for traditional ESRI cache created on ArcGIS Server.

I have created two sample tile packages to test the workflow.  The source data used for one of these was in a postgis database.  The source data used for the other was in an esri fgdb.

Both vector tile packagaes were created in ArcGIS Pro 1.4.  The tiling Scheme used was "VTTS_27700_British_National_Grid" and the Tiling Format was set to "Indexed".

I am able to view both .vtpk files in ArcGIS Pro (v1.4) and in ArcGIS Online.

I have so far been unable to get the attribution to come through when I create the tile packages using the Create Vector Tile Package tool.

According to this ESRI help page I should be able to bring through the attribution, however, I have to specify the fields that I want to bring through.

I have followed the steps suggested in the article and ticked the "Highlight" box for all the fields I want to carry over to the Vector Map Tiles.

Once the tiles have been created, if I click on a feature I get no attribution & no pop-up.

Has anybody managed to do this?


texts: Keep attribute data while converting XML file to shapefile in QGIS?

texts:
I want a map of particular area in shapefile format from Google Earth. I created a map using polygon method and save it as KML file format, but when I convert it from ".xml" to "shapefile" format using QGIS Desktop software I am not getting attribute data in the converted shapefile.

What should I do to keep my attribute data in shapefile format?

I am attaching image which showing no attribute data in QGIS:


texts: Get shortest distance between two geometries in ArcPy

texts:
I'm using ArcPy and using a cursor to step through one polygon feature class. For each of these 'source' polygons, I am creating a SearchCursor to find all polygons in a second feature class that share a common attribute. For each of these matching polygons, I want to measure the distance between it and the search polygon. I also plan on measuring other geometry differences (to determine if it has just been shifted, rotated, scaled, etc).

The purpose of this work is to build up a change history between polygon datasets representing features at different times.

I've seen the ReturnDistance method mentioned in ArcObjects, but can't figure out how to run this using ArcPy. Is it possible, and if so, can anyone help with some rough sample code? If not, I'll probably resort to measuring the change in centroid position between the polygons.

Thanks, Steve


texts:
I have uploaded some data to arcgis online to create a web map, even converted to WGS84 from BNG before zipping, I am using QGIS to create my maps, I do have a desktop licence for ArcGIS but trying to use QGIS to not rely on ArcGIS desktop. You can see after loading to arcgis Online the features have shifted to the left and south.

See pictures below.


texts: Managing Conflicts in Scripted ArcGIS Reconcile of Versions?

texts:
We have some integration between our Assets system and our ArcGIS 10.3.1 system.  The ArcGIS SDE is versioned.  The assets system uses its own version (lets call it 'Assets'), and there is another version, 'QA', that all the editors use (actually, they use their own child versions of 'QA', but lets keep it simple for now).

The assets system only integrates non-spatial attributes with the GIS features, not the geometries of features.  In general, users are not permitted to change non-spatial attributes in the GIS, but only in the assets system.  (This is a simplification, but will do for the purposes of this question.)

However, they can (and must) change spatial data (geometries) in the GIS.

Every night, an automated Python script reconciles and posts both of these versions to the DEFAULT version.  In fact, the script runs twice every night so that the changes from each version get reconciled back to the other version, via DEFAULT.

The problem here is that if spatial data (geometry) is changed in the GIS and a non-spatial attribute is changed in the assets system on the same day, the reconcile/post will overwrite the changes in one of these (which ever gets reconciled last 'wins' and the changes in the other are lost).

The Python arcpy script uses the following line to run the reconcile/post:

arcpy.ReconcileVersions_management(workspace, "ALL_VERSIONS", "sde.DEFAULT", versionList, "LOCK_ACQUIRED", "NO_ABORT", "BY_OBJECT", "FAVOR_EDIT_VERSION", "POST", "KEEP_VERSION", reconcileLog)

(where versionList is a list of 'QA' and 'assets' versions).

I could change "NO_ABORT" to "ABORT_CONFLICTS", but that would just result in no reconcile/post occurring.

Is there any way to automatically log/notify of conflicts during this process?  If set to "NO_ABORT", conflicts don't appear to be logged at all!

Even better, is there some way to automatically manage conflicts (eg, rule-based of some kind)?

Better still, when a conflict is detected between versions, is there some way to automatically retain the spatial data of the feature from one version and the non-spatial data of that feature from the other version?


texts:
Possible Duplicate:
  How to install Pyton libs for XYtools to read Excel files

I was using XY Tools plug-in for QGIS that converts Excel files into attribute table and shape files etc, but once I select the file for conversion this appears:

There is a problem with loading a python module which is needed to read/write Excel files. Please see documentation/help how to install python xlw and xlrd libraries.

What should I do?


texts: How do you calculate latitude and longitude of a point in a static map?

texts:
Background: I have generated a static map (map.png) of a particular city using mapnik and osm data of the highest zoom level. I am able to zoom in and zoom out successfully using my Java program. It is an offline desktop application.

Question: I want to represent a point at a particular latitude and longitude on this map. How do I do it?

I went through the following references:
http://wiki.openstreetmap.org/wiki/Zoom_levels
http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Resolution_and_Scale
How to plot a point on a static Google map [png]?

I'm still not able to understand how to plot it on a static map by calculating a specific latitude and longitude at a particular zoom level with the image having a specific pixel width and height.


texts: Counting and rasterizing polygon overlaps in QGIS?

texts:
I have a large number of polygons that are overlapping. I'd like to create a heat map and plan on using raster calculation to do this. My plan: to give each polygon a raster value of 1, and adding them up. But, when I convert to raster, it will only the entire shape so that all overlapping polygons have the raster value of 1.

Can I make it convert the individual polygons? I'm doing this in QGIS but come from an ArcGIS background, so either explanation may give me enough to get this done.


texts: How do I fix an error starting QGIS 1.8.0 on Windows XP?

texts:
I note that the installation of qgis completed successfully. However, the log procedure returns this message:

  "...... Estrazione: wvhpgn.gsb Estrazione: wyhpgn.gsb Cartella di
  destinazione: C:\Programmi\Quantum GIS Lisboa Creato il programma di
  disinstallazione: C:\Programmi\Quantum GIS Lisboa\Uninstall-QGIS.exe
  Crea cartella: C:\Documents and Settings\All Users\Menu
  Avvio\Programmi\Quantum GIS Lisboa
  Impossibile accedere al file. Il file Š utilizzato da un altro processo. Running postinstall grass64.bat... Running postinstall
  msvcrt.bat... ......"

  System parameters: Quantum GIS Lisboa 1.8.0  Microsoft Windows
  Professional, Versione 2002, Service Pack 3

The message:

  Impossibile accedere al file. Il file Š utilizzato da un altro
  processo

translates as

  Unable to access the file. The file used by another process

I then get the following error:

  "Impossibile trovare il punto di ingresso locale_charset della
  procedura nella libreria di collegamento dinamico iconv.dll"

which translates as:

  Can not find the entry point of locale_charset procedure in the
  dynamic link library iconv.dll


texts: Strange banding in output raster after weighted overlay

texts:
I am trying to conduct a weighted overlay. I have two rasters, one is a continuous slope raster that I reclassified to have 6 classes (8-bit signed integer), and one is an elevation raster that I also reclassified to 6 classes (also 8-bit signed integer).

My inputs look like this:

My output looks like this:

Any idea what could possibly be going wrong?

I've never had an output like this before.


texts: Export to raster tool changes the output values

texts:
I am working with interpolated temperature data. The image's maximum and minimum values are 11.9 and -1.7, respectively.

I exported this image to raster using export to raster, but the max and min values of the output raster image are different from that of the interpolated temperature image.

I'm wondering why values of both images (temperature and raster) are not the same. They were supposed to be. Is there any way to make the values of the raster image be the same as the interpolated image?


texts: Getting daily-air temperature of several months of past years for Australia?

texts:
I really need to obtain the daily-air temperature for several months of 2011. In addition, the data I'm looking for must have the maximum spatial resolution of 1Km*1Km. Since, as I know, there isn't any reliable method to estimate air temp. from satellite imagery, I have to get the data from weather stations available on the internet. However, I haven't found any sites to help me what I need.

Does anyone know how can I get the data for Australia?


texts: How to add hypothetical values to a raster?

texts: I have to create Non-Regulatory Products for a few FEMA regions. Right now I am working on the 'Increased Flooding Scenarios' where you add a hypothetical value to the existing Base Flood (1', 2', and 3', specifically). I already created the Base Flood (or Coastal Depth Grid), now I am stuck on how to add a hypothetical value to it (1 foot) and have a 'plus 1' raster extending to the right elevation. I've tried various methods where the elevation increases by a foot, but the raster does not extends to the right elevation (Con, Plus, Raster Calculator). Attached is a  flowchart in FEMA's guidelines on how to create the flooding scenarios.


texts:
I'm using gdal.ReprojectImage() to reproject mod021km bands to geotif data. it takes one of the modis multiband subdatasets (e.g. EV_500_Aggr1km_RefSB with 5 bands) and saves it as a tif file with one band as a result (i did test it and get result). as far as i need to reproject every band, i tried to save each band of subdataset ('tmp_band.hdf') at first, read it by gdal.Open() and feed it to the gdal.ReprojectImage but i got these errors:

  RuntimeError: The transformation is already "north up" or a transformation between pixel/line and georeferenced coordinates cannot be computed for tmp_band.hdf. There is no affine transformation and no GCPs. Specify transformation option SRC_METHOD=NO_GEOTRANSFORM to bypass this check.

the code i used to save bands and feed them to gdal.ReprojectImage is:

def _reprojectOne(self, l, quiet=False):
    """Reproject a single subset of MODIS product
    l = complete name of input dataset
    """
    l_src_ds = gdal.Open(l)
    number_of_bands = l_src_ds.RasterCount
    meta = l_src_ds.GetMetadata()
    band = l_src_ds.GetRasterBand(1)
    bands = l_src_ds.ReadAsArray()
    if '_FillValue' in list(meta.keys()):
        fill_value = meta['_FillValue']
    elif band.GetNoDataValue():
        fill_value = band.GetNoDataValue()
    else:
        fill_value = None
    datatype = band.DataType
    try:
        l_name = l.split(':')[-1]
        out_name = "{pref}_{lay}".format(pref=self.output_pref,
                                             lay=l_name)
    except:
        out_name = "{pref}".format(pref=self.output_pref)
    if self.vrt:
        out_name = "{pref}".format(pref=self.output_pref)

    for i in range(number_of_bands):
        tmp_ds = l_src_ds.GetDriver().Create('tmp_band.hdf',
                                    l_src_ds.RasterXSize,
                                    l_src_ds.RasterYSize,
                                    1, datatype)
        m = tmp_ds.GetRasterBand(i+1)
        m.WriteArray(bands[i,:,:])
        tmp_ds.SetGeoTransform(l_src_ds.GetGeoTransform())
        tmp_ds.SetProjection(l_src_ds.GetProjection())

        m.FlushCache()
        tmp_ds = None

        try:
            tmp_src_ds = gdal.Open('tmp_band.hdf')
            dst_ds = self.driver.Create(out_name+'_B{b}.tif'.format(b=i+1),
                                        self.dst_xsize,
                                        self.dst_ysize,
                                        1, datatype)
        except:
            raise Exception('Not possible to create dataset %s' % out_name)
        dst_ds.SetProjection(self.dst_wkt)
        dst_ds.SetGeoTransform(self.dst_gt)
        if fill_value:
            dst_ds.GetRasterBand(1).SetNoDataValue(float(fill_value))
            dst_ds.GetRasterBand(1).Fill(float(fill_value))
        cbk = self._progressCallback
        # value for last parameter of above self._progressCallback
        cbk_user_data = None
        try:
            gdal.ReprojectImage(tmp_src_ds, dst_ds, tmp_src_ds.GetProjection(),
                                self.dst_wkt, self.resampling, 0,
                                self.error_threshold, cbk, cbk_user_data)
            if not quiet:
                print("Layer {name} reprojected".format(name=l))
        except:
            raise Exception('Not possible to reproject dataset '
                            '{name}'.format(name=l))
        dst_ds.SetMetadata(meta)
        dst_ds = None
        l_src_ds = None

where is the problem? how can i feed the specific band to the gdal.ReprojectImage?
should i save more information in 'tmp_band'(information from subdataset metadata)?

Note: the original code is a script named convertmodis_gdal.py that is used by modis_convert.py


texts: What are the encodedpath and simplepath properties of a GeoJSON geometry?

texts:
I have GeoJSON FeatureCollection with an encodedpath and simplepath property instead of a coordinate array.

For example:

    "geometry": {
        "type": "LineString",
        "encodedpath": "sdecF`noxMq@hBKbAH\\r@l@ZFfAKf@b@",
        "simplepath": "sdecFbnoxM}@hDF\\~DnA"
    }

When viewed in QGIS, I can view attributes but no geometry.

I don't see this in the spec, but I am guessing it is a newer draft?

What are these properties?


texts: Executing a standalone QGIS python script on Windows 10

texts:
In Windows 10, I currently have two working scripts (.BAT and .PY) which are used to generate PDFs from QGIS based on the solution here.

My goal is to be able to execute the python script without having to call the batch commands each time, but am having trouble getting this to work.

The createQGISmap.bat file contains:

REM Change OSGeo4W_ROOT to point to your install of QGIS.

SET OSGEO4W_ROOT=C:\Program Files\QGIS 2.18
SET QGISNAME=qgis
SET QGIS=%OSGEO4W_ROOT%\apps\%QGISNAME%
set QGIS_PREFIX_PATH=%QGIS%

CALL "%OSGEO4W_ROOT%\bin\o4w_env.bat"

REM Python Setup
set PATH=%OSGEO4W_ROOT%\bin;%QGIS%\bin;%PATH%
SET PYTHONHOME=%OSGEO4W_ROOT%\apps\Python27
set PYTHONPATH=%QGIS%\python;%PYTHONPATH%

ECHO OSGeo path is: %OSGEO4W_ROOT%
ECHO Getting QGIS libs from: %QGIS%
ECHO Python loaded from: %PYTHONHOME%

"C:\Program Files\QGIS 2.18\bin\python.exe" getScreenshot.py

And the beginning of the getScreenshot.py script:

from qgis.core import *
import sys
from qgis.gui import QgsMapCanvas, QgsLayerTreeMapCanvasBridge
from PyQt4.QtCore import QFileInfo, QObject
from PyQt4.QtXml import QDomDocument

gui_flag = True
app = QgsApplication(sys.argv, gui_flag)

# Make sure QGIS_PREFIX_PATH is set in your env if needed!
app.initQgis()

# project.qgs
project_path = 'C:/mydata/qgis/autoexport.qgs'

I tried manually specifying the environment variables in Windows (in C: - Properties - Advanced system settings - Environment variables...) as follows:

QGIS = C:\Program Files\QGIS 2.18\apps\qgis
QGIS_PREFIX_PATH = C:\Program Files\QGIS 2.18\apps\qgis
PATH = C:\Program Files\QGIS 2.18\bin
PATH = C:\Program Files\QGIS 2.18\apps\qgis\bin
PYTHONHOME = C:\Program Files\QGIS 2.18\apps\Python27
PYTHONPATH = C:\Program Files\QGIS 2.18\python

But when executing python getScreenshot.py from the OSGeo4W shell, I keep getting

from qgis.core import *
ImportError: No module named qgis.core

I tried following advice from How can I fix "ImportError: No module named qgis.core"? and here without luck. I should point out I have two installations of Python (the one bundled with QGIS and a Python 3 from Anaconda), so I am trying to run the script with the former installation.


texts: How to fix "import qgis.core ImportError: DLL load failed"

texts:
I am trying to develop a GIS application in QGIS using the Python language. Unfortunately, I encountered a problem from the outset is that of import and qgis.core qgis.gui. the Problem in the python console is as follows:

import qgis.core

Traceback (most recent call last):
  File "<pyshell#0>", line 1, in <module>
    import qgis.core
ImportError: DLL load failed: Le module spécifié est introuvable.

I am using:

OS : Windows 7 , 32 bits
QGIS : 1.8 ( Lisbona )
Python : python 2.7


texts: quite different results in spatialEco::tpi with circular and rectangular window

texts:
In R, I have an elevation raster. I want to calculate the topogographic position index using the spatialEco::tpi function. The CRS of the elevation raster is "+init=epsg:5677"; the CRS unit thus should be "meter". The resolution of the elevation raster is
> raster::res (elevation_raster)
[1] 20.9 30.9

This is the elevation raster:

spatialEco::tpi with win = "rectangle"  and scale = 5 should yield similar results as spatialEco::tpi with win = "circle" and scale = 150. At least, I thought so, because 5 * 30 = 150. However, I get completely different values and the plots also look quite different:
spatialEco::tpi (elevation_raster, win = "rectangle", scale = 5)

spatialEco::tpi (elevation_raster, win = "circle", scale = 150)

Notice that the two plots look quite different but also that there is a huge difference in the value ranges. Apparently, I am missing something here. However, I have no clue what. Perhaps someone can give a proper explanation?


texts: lidR: no slot of name "header" for this object of class "LAScatalog"

texts:
edit: This question does share an answer with an existing question, and the problem is the same: lastrees() does not work on LAScatalogs. However, I'm going to leave this question up, because I think the code and error message makes it easier to find.

I'm using lidR 2.0.3 (from the rocker/geospatial image). I would like to perform tree segmentation on an entire LAScatalog using the Li 2012 algorithm:

#### setup #####################################################################

# load packages
library(here)
library(lidR)

# load normalized point cloud catalog (just a pointer to the tiled point
# clouds)
normalized_als =
  readRDS(here::here('02-data',
                   'lidar',
                   'aerial',
                   'normalized_als_catalog.rds'))

#### configure catalog #########################################################
# set catalog output path
lidR::opt_output_files(normalized_als) =
  here::here('02-data',
             'lidar',
             'aerial',
             '02-segmented',
             '{ID}_segmented')

# don't use any points with Z below 0 when making the CHM
lidR::opt_filter(normalized_als) =
  paste0(lidR::opt_filter(normalized_als), ', drop_z_below 0')

# use n-2 cores to process
lidR::opt_cores(normalized_als) = as.integer(parallel::detectCores()-2L)

# use a 20m buffer
lidR::opt_chunk_buffer(normalized_als) = 20

# process by files - this may preserve the lasheader?
lidR::opt_chunk_size(normalized_als) = 0

# run the algorithm, writing the segmented point cloud to disc
segmented_als =
  lidR::lastrees(normalized_als,
               algorithm = li2012())

The last call to lastrees returns:

Error in algorithm(las):
    no slot of name "header" for this object of class "LAScatalog"

I've confirmed that the individual las files DO have valid headers, but it seems that lidR::li2012() is assuming that "las" is a LAS object, not a LAScatalog. The call fails at line 219 of lidR/R/algorithm-its.R:

 if (las@header@PHB$`Max Z` < hmin)

My question is: Is my code incorrect, or is this a bug in lidR::li2012()? I'm trying to follow the examples on the wiki, but I don't 100% understand how lidR maps LAS-based functions over LAScatalogs, and it's very possible that I am misusing the lastrees function. (Is there some catalog_apply() approach that I should be using?) The li2012() function seems to expect a LAS object specifically, not a LAScatalog (which doesn't have a comparable header slot). How can I pass an actual LAS object to li2012()? I would just iterate over the separate tiles, but would like to use the buffering functionality offered by LAScatalogs. Alternatively, was the data validation checking for li2012() written with a LAS object in mind, and causing it to break when given a LAScatalog or a LAScluster?

I would be happy to email a reproducible dataset, but I believe the issue should be reproducible with an arbitrary catalog: Either my code is written incorrectly, or li2012() can't currently be used on LAScatalogs.


texts: Intersect function is generating different outputs in ArcMap and ArcPy

texts:
I am currently running the following code in ArcPy:
in_dir_temp = "D:\\X\\"
x1= in_dir_temp + "x1.shp"
x2= in_dir_temp + "x2.shp"
y1= in_dir_temp + "y1.shp"
arcpy.Intersect_analysis([x1, x2], y1, "ALL", "", "INPUT")

Where x1 and x2 are polygons and so is y1.
When running this in ArcPy, it does not generate a shapefile in the output. It only generates .sbn .xml and .sbx files, despite the code outputting
<Result 'D:\\X\\y1.shp'>

However when I manually intersect the two polygons in ArcMap, it generates a polygon and shapefile.
I can't figure out why I am getting two different outputs?


texts: Generating properly-aligned shape from land survey description?

texts:
I'm looking through a variety of land surveys from construction permits in Florida, and I'm trying to generate a set of lat/lon points from them, in order to map them out. The surveys list a PLSS reference as a starting point, and give bearing and distance all around the shape.

From what I've read, I think I should be able to use the Vincenty formula to generate points based on that information. I've got PLSS section corner points from Geocommunicator to use as starting points, and I'm using the WGS-84 constants for the formula.

I've tried Geopy's implementation as well as a Python port of a JavaScript implementation, and they both suffer from the same problem: the overall shape looks right, but it's not aligned correctly. Here's what it looks like in JOSM, overlaying the resulting shape (including the path from the PLSS point to the survey's Point of Beginning) on aerial photography for reference.

Why is this happening? The shape looks right and the size fits, but it's mis-aligned. My first thought is that I had wrong PLSS data for my starting point, but the alignment problem isn't uniform. It's usually off by within 200-300 feet, and mostly east of where it should be, but sometimes a little west instead, and the north/south offset varies as well, but usually it's not as dramatic as the east/west misalignment.

Is it possible that I just need better PLSS points to start from? Or should I be using some other approach to calculate the points, besides the Vincenty formula? Do I need to reproject something? I'm still learning all the tools here, and I'm sure I'm just missing some step of the process.

Update:

In response to some requests for more information, the PLSS shapefile I have is in NAD83 coordinates (prj2epsg identifies it as EPSG:4269), which identifies the corners of each section within the county. I then used simple arithmetic to identify the midpoints where appropriate, as some of the surveys use those instead of the outer corners. I've tried reprojecting into EPSG:4326 using ogr2ogr, but it doesn't seem to change the coordinates at all.

Also, the survey data consists of bearings and distances, which I think is different from metes and bounds. For example, the permit identified by the above image is as follows:

Point of reference: South quarter corner of Section 11 of Township 24 South, Range 27 East off the Talahassee Meridian. According to my PLSS shapefile and midpoint math, that's 28.4056197, -81.58239195.

From there, here are the measurements:

N 89° 41' 52" E 237.44 ft
N 00° 00' 00" E  31.87 ft
Point of Beginning
N 74° 54' 56" W 160.57 ft
N 15° 05' 04" E  14.09 ft
N 74° 54' 56" W 239.99 ft
N 15° 05' 04" E  92.44 ft
S 74° 54' 56" E 213.61 ft
N 15° 05' 04" E  13.29 ft
S 74° 54' 56" E  36.43 ft
S 15° 05' 04" W  27.38 ft
S 74° 54' 56" E 150.53 ft
S 15° 05' 04" W  92.44 ft

I've also tried making sure I'm using survey feet, rather than standard feet, but the difference is minuscule at these distances.

Also, the image comes from Bing aerial photography, so I assume it's EPSG:3857. JOSM uses that as a fairly standard layer for tracing, so I also assume it's done whatever work is necessary to line it up with lat/lon points properly.