texts:
in ArcGIS 10.3 I have 50 simple models I would like to run one after another. Im confused on the submodel process, currently I run each clicking from arccatalog, 50 in all. not a big deal but Id like to just click once and let it fly.

Note: Each Model has 1 iteration process within the Model.

currently Ive only run each model at a time, Ive tried to reasearch how to string them together via various website & help, but the iteration it seems is the stumbling issue.


texts: Rotating point feature symbols by intersect polyline in ArcGIS Desktop?

texts:
I would like to rotate point features symbols by angle of polyline that they overlap. Similar to cartographic line symbol, that can rotate to follow line angle:

I´m using basic license of ArcGIS 10.1. I´m not able to use Points along one or more lines at a time, with rotation solution. It was made for ArcGIS 10.3. I think that method arcpy.AddGeometryAttributes_management can be used from ArcGIS 10.2.

I followed this approach:

Create buffer Find upper and lower intersection (buffer x line)
Calculate angle
Update point feature attribute with calculated angle
Rotate point by attribute with calculated angle

Solution could by something like this:


texts: Mapping raw GPS data to Openstreetmap's Node/Way/Relation

texts:
I have a bunch of raw GPS data points which constitute trips.
I want to add semantic information to those data points. For example I'd like to be able to say that a trip consisted of taking street A, B and C.

With Openstreetmap I can download nodes and ways for a bounding box, which are again described through a bunch of nodes (lat/lng) points.

I actually have all the information I need to map my points to those of Openstreetmap. But reality is messy and my GPS data is not perfect, has outliers and goes off the road.

What is a good approach to calculate which ways or even relations belong to my GPS data?

What is the research done in this area?

Is there already a library that does all the work for me?


texts: SVG created by QGIS Profile tool won't scale to include the vertical units!

texts: I exported an svg image from the QGIS Profile tool but when I go to use it in a Picture on a map layout it will not scale to show both the X and Y units.  It's like they are not considered as part of the svg image for zoom or stretch purposes.  Any thoughts?  Photos attached.


texts: Querying subset of large ArcGIS server dataset via REST

texts:
I'm currently set up as follows:

A large wetlands polygon dataset is stored in an arcgis.com account as a feature service (~395,000 polygons).

I'm using the Esri-leaflet javascript library in conjunction with leaflet to dynamically load polygons into my map by querying using the map bounds and limiting by zoom level (in order to avoid exceeding the feature limit). I set a minimum zoom level of 14 this purpose.

When I use a smaller dataset that is also stored on arcgis.com (~17,000 features) it works flawlessly but when I switch the URL to the larger dataset (~395,000) features it returns an empty dataset.

Here is a snippet of the working code with the "small" dataset test:

var wwi = L.esri.featureLayer("http://services.arcgis.com/F7DSX1DSNSiWmOqh/arcgis/rest/services/WAWFA_Small_test/FeatureServer/0",{
    precision: 4,
    style: function(feature){
        return {
            color: layerColorRamp(feature.properties,'wetlands','SP'),
            weight: 1.5,
            opacity: 1,
            fillOpacity: 0.7
        }
    },
    simplifyFactor: 0.3,
    minZoom: 14
}).addTo(map);

The large dataset uses this url: http://services.arcgis.com/F7DSX1DSNSiWmOqh/arcgis/rest/services/WAWFA_Large_test/FeatureServer/0

but should otherwise be identical. When I look at the network tab of Chrome's web tools, I see that the actual REST query looks like this:

http://services.arcgis.com/F7DSX1DSNSiWmOqh/arcgis/rest/services/WAWFA_Small_test/FeatureServer/0/query?returnGeometry=true&where=1%3D1&outSr=4326&outFields=*&inSr=4326&geometry=%7B%22xmin%22%3A-88.11035156249999%2C%22ymin%22%3A43.100982876188546%2C%22xmax%22%3A-88.08837890625%2C%22ymax%22%3A43.11702412135048%2C%22spatialReference%22%3A%7B%22wkid%22%3A4326%7D%7D&geometryType=esriGeometryEnvelope&spatialRel=esriSpatialRelIntersects&geometryPrecision=4&maxAllowableOffset=0.00003328443215608169&f=geojson

which produces a working geojson output but when I switch the prefix to http://services.arcgis.com/F7DSX1DSNSiWmOqh/arcgis/rest/services/WAWFA_Large_test/FeatureServer/0/    for the larger dataset I get an output like this:

{ "type": "FeatureCollection",
  "crs": { "type":"name",
           "properties": { "name": "EPSG:4326" }
         },
  "features":[] }

Any ideas on why this might be happening and whether there is a way to address this?


texts: Where can I find UK boundary data shapefiles where features have attributes pointing to their parent subdivision?

texts:
I know similar questions have been asked before but none have given me exactly what I need. I am trying to find boundary shapefile data for the UK which splits different subdivisons into different layers (the exact number of subdivisions doesn't matter), with each feature of a layer referencing it's parent in the above layer.

I would like the data to be 'clean' in that each layer only has appropriate regions, e.g. a county and city aren't on the same layer, so that all features of a layer are of a similar size and no features are surrounded entirely by another feature on the same layer (an enclave).

The closest I have got so far is at gadm.org, which provided shapefiles with perfect attributes, however Unitary Authorities were listed separately to their parent county on the same layer (e.g. Nottingham was an enclave within Nottinghamshire on the county layer, and thus in the next layer down Nottingham's parent layer was Nottingham, not Nottinghamshire). I could merge these and change the data for the parent, but this would have to be done manually which isn't preferable.

Also, this is not a duplicate of Seeking administrative boundaries for various countries?, I am asking for much more specific data.


texts: Create line/points at contacts between polygons of certain attribute?

texts:
I have a shapefile of polygons representing geologic map units. I want to create line/point features along the contact between two particular units (polygons identified by attribute code) so I can get those elevations off a DEM.

How can I create line features only along the contact between specific polygon attributes?


texts: reduceRegion not working in Google Earth Engine

texts:
I want to get the NDVI value for each point of my featureCollection vector. I had used reduceRegion for that in codeeditor javascript of GEE. I when i using the same function in python api of GEE, it is showing 'Unrecognized argument type to convert to a FeatureCollection'. May be some thing need to change in python api for the reducer function. Please help me in this regards.

import ee

import time
import datetime
import json
from ee.batch import Export

ee.Initialize()

startyear=2018
startmonth=01
startday=01

endyear=2018
endmonth=02
endday=01
def maskS2clouds(image):
    qa = image.select('QA60');
    cloudBitMask = ee.Number(2).pow(10).int();
    cirrusBitMask = ee.Number(2).pow(11).int();
    mask = qa.bitwiseAnd(cloudBitMask).eq(0) and(
             qa.bitwiseAnd(cirrusBitMask).eq(0));
    return image.updateMask(mask).divide(10000);

print("test")

boundary = ee.FeatureCollection('ft:1T-XrxMki2k49fhQdhMxesWKBoyyUJiKN70pErh9z')
#print(boundary)
gtpoints = ee.FeatureCollection('ft:1ZvFeFO5jHGo3VExlKNc-kdL-3hmJSeyfLzqLvHhy')
print (gtpoints)
image = (ee.ImageCollection('COPERNICUS/S2').
         filterBounds(boundary).
         filterDate(datetime.datetime(startyear,startmonth,startday),
                     datetime.datetime(endyear,endmonth,endday)).
         sort("CLOUD_COVERAGE_ASSESMENT").
         map(maskS2clouds).
         limit(20))

mosaic = image.mosaic();
    #median = mosaic.median();
palette = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
               '74A901', '66A000', '529400', '3E8601', '207401', '056201',
               '004C00', '023B01', '012E01', '011D01', '011301'];
clip = mosaic.clip(boundary)
#median = clip.median()
#print(clip)
ndvi = clip.normalizedDifference(['B8', 'B4'])
NDVI=ndvi.multiply(100).uint8()
#print(clip)
palette = ','.join(palette)
vis = {'min': 0, 'max': 300, 'bands':'nd', 'palette': palette}
toexport = NDVI.visualize(**vis).addBands(NDVI);
#print(toexport,'toexport')
points_ndvi = ndvi.reduceRegions({
  'collection':gtpoints,
  'reducer':ee.Reducer.first(),
  'scale':10
})
properties = ee.Feature(points_ndvi.first()).propertyNames()
new_names = properties.replace('first', 'ndvi')
points_ndvi = points_ndvi.select(properties, new_names)
print(points_ndvi,"test")


texts: Loading a KML Layer Underneath a Fusion Table Layer

texts:
I am working on a map that will incorporate both Fusion Tables and KML layers, parsed through geoXML3. As my Fusion Table Layer is a point file I would like it to always be on top of the KML's which will be polygons. By default when the page first loads the Fusion Table will be loaded, while the KML's are turned on by the user through some check boxes. It seems like because the Fusion table loads first it always gets stuck on the bottom. I tried to fix this by turning it off when parsing/ displaying the KML data and then turning it back on after the call to parse the data but this does not seem to work. Does anyone have an idea of how I could make it so the Fusion Table Layer is always on top?

Below is a sample of my code right now which doesn't work:

    function initialize() {
var myLatlng = new google.maps.LatLng(43.705199, -79.421707);
var myOptions = {
zoom: 12,
center: myLatlng,
mapTypeId: google.maps.MapTypeId.ROADMAP
};
map = new google.maps.Map(document.getElementById('map_canvas'), myOptions);

kml_FOS = new geoXML3.parser({map: map, zoom: false});
kml_Neighbourhood = new geoXML3.parser({map: map, zoom: false});
kml_Municipal = new geoXML3.parser({map: map, zoom: false});
kml_SSI_University_Aggregate_V01 = new geoXML3.parser({map: map, zoom: false, suppressInfoWindows: false});

roadFusionLayer2 = new google.maps.FusionTablesLayer(fusionTableID,{suppressInfoWindows: false});
roadFusionLayer2.setMap(map);
}

function boundaryToggle_FOS() {
if (boundaryIsLoaded_FOS == 0)
{
boundaryIsLoaded_FOS = 1;
    if (firstLoadFOS == 0)
        {
        //turn off fusion table layer
        roadFusionLayer2.setMap(null);
        //Parse the kml
        kml_FOS.parse('FOS_Boundary_v03.kml');
        //Turn Fusion Table ON
        roadFusionLayer2.setMap(map);
        firstLoadFOS = 1;
        }
    else
        {
        roadFusionLayer2.setMap(null);
        kml_FOS.showDocument();
        roadFusionLayer2.setMap(map);
        }
}
else
{
    kml_FOS.hideDocument();
    boundaryIsLoaded_FOS = 0;
}
}


texts: Reading/writing local datastore for QGIS?

texts:
What is a solution for reading/writing a local datastore for QGIS?

I've got data that's too large for a shapefile, I need lengthy column titles, but don't want to be logging into a remote server.
I've been giving Spatialite a good go but seem to bump up against limitations like not being able to delete columns plus it seems there is little development going on in this space.  I use the QSpatialite plugin.

I'm thinking I should be getting in db views, but don't see any workflows there for Spatialite.

Is there another option for a FOSS local datastore?


texts: How to restrict panning in QGIS to a given extent

texts:
I'm going to build a QGIS standalone application. According to the requirements full extent equals to the extent of a predefined rectangle (say bounding box of a parcel). The application should start in full extent view. When zoomed in the user is restricted to only pan inside the area given by full extent (using pan tool or sliding scrollbars), similar to the issue in this question.

Which is the appropriate QGIS API function or how can I mimic this behaviour using Python?


texts: Georeference image using center point and size?

texts:
I have looked through the forum extensively but found no answer. Apologies in advance if this has been addressed somewhere already...

I have >500 png image files (from a microscope) that I would like to display in their correct relative positions. Each image file has an associated xml file that contains:
1. the center XY position (in microns)
2. the XY size of the image (in microns)
3. The amount of rotation (if any, in degrees)

The images are in an arbitrary positive cartesian coordinate system.

I have looked at the Freehand Georeferencer and Georeference plugins, but they do not allow referencing based on the middle of the image plus image size (as far as I can tell).

Is it possible to use the information in the xml file to georeference the png files in QGIS?


texts:
I am using ArcMap.

I am trying to convert polygon to raster and specify cell size to 6, but created raster is empty. If I leave the default cell size (0.0015) everything works. Last week I was doing the same thing on the same data and everything was working.

What can be wrong?

Later I need to convert this raster (clipped with another one) to ASCII. I've tried to resample raster which I had created previously (cell size 16). Works good, but then have to change the projection and cell size is automatically changed. After projecting raster I can't resample it again to 6 (error with too small or too big cell size appears). Even, if I leave the raster with unchanged cell size I can't convert it to ASCII (the process stays at 0%).


texts: Seeking data of worldwide shore lines and bathymetry?

texts:
Is there a free database with all geolocations for shorelines?

Is there a free database for near-shore bathymetry data?


texts: Accessing Terrain Tiles using podpac - broadcast error

texts:
I am working through the examples from the podpac site. I have stated with the most basic example for using terrain tiles given here: https://github.com/creare-com/podpac-examples/blob/main/notebooks/5-datalib/terrtain-tiles.ipynb
from podpac.datalib.terraintiles import TerrainTiles
from podpac import Coordinates, clinspace

# create terrain tiles node
node = TerrainTiles(tile_format='geotiff', zoom=8)

# create coordinates to get tiles
c = Coordinates([clinspace(40, 43, 1000), clinspace(-76, -72, 1000)], dims=['lat', 'lon'])

# evaluate node
o = node.eval(c)

I get the following error: ValueError: operands could not be broadcast together with shapes (1000,1000) (0,0)
A similar error occurs for all coordinates/dimentions/zooms that I attempt.


texts: IO error on bulk MXD publishing to ArcGIS Server using Python

texts:
I have a folder of MXDs and I want to publish all of them in batch to a single ArcGIS Server. I've had success using Python on a single MXD, but when I add a loop into the code so that I can do multiple services at once, I run into errors that I'm unable to troubleshoot. I found the code below on ArcGIS Enterprise site and tried to modify it for batch processing.
import arcpy
import os

wrkspc = 'C:/Project/BULK/'
#mxd = arcpy.mapping.MapDocument('BULK_Test.mxd')

con = 'AGSconnectionfile.ags'

for (wrkspc, dirs, files) in os.walk(wrkspc):
    for fl in files:
        if fl.lower().endswith(".mxd"):

            mxd = arcpy.mapping.MapDocument(os.path.join(wrkspc, fl))
            # Provide other service details
            service = mxd
            sddraft = 'C:/Project/BULK/' + str(mxd) + '.sddraft'
            sd = 'C:/Project/BULK/' + str(mxd) + '.sd'
            summary = 'BULK_Test'
            tags = 'BULK_Test'

            # Create service definition draft
            arcpy.mapping.CreateMapSDDraft(mxd, sddraft, mxd, 'ARCGIS_SERVER', con, True, "Project", summary, tags)

            # Analyze the service definition draft
            analysis = arcpy.mapping.AnalyzeForSD(sddraft)

            # Print errors, warnings, and messages returned from the analysis
            print "The following information was returned during analysis of the MXD:"
            for key in ('messages', 'warnings', 'errors'):
                print '----' + key.upper() + '---'
                vars = analysis[key]
                for ((message, code), layerlist) in vars.iteritems():
                    print '    ', message, ' (CODE %i)' % code
                    print '       applies to:',
                    for layer in layerlist:
                        print layer.name,
                    print

            # Stage and upload the service if the sddraft analysis did not contain errors
            if analysis['errors'] == {}:
                # Execute StageService. This creates the service definition.
                arcpy.StageService_server(sddraft, sd)

                # Execute UploadServiceDefinition. This uploads the service definition and publishes the service.
                arcpy.UploadServiceDefinition_server(sd, con)
                print "Service successfully published"
            else:
                print "Service could not be published because errors were found during analysis."

            print arcpy.GetMessages()

I've tried as a standalone script and from within ArcMap 10.7.1. The error:
File "C:\Project\BULK\PublishAGSBulk.py", line 29, in <module>
IOError: Operation on file C:\Users\...\Temp\{67337BAB-1DDF-4A39-9FB4-B7C8EF8597B0}\<geoprocessing Map object object at 0x036534A0>.sddraft failed. Unknown error.

I am doing this with an active connection to the ArcGIS Server (10.8) in ArcCatalog.


texts: Rendering freehand drawn geometry from geojson on map with OpenLayers?

texts:
I want to render freehand drawn geometry from a geojson object on a map. I have the freehand drawing all set up and working. The geometry is LineString but I am adding the freehand: true property which works. Now   after getting the coordinates on 'drawend' event using

var shape = event.feature.getGeometry();
var coords = shape.getCoordinates();

I am saving the coordinates e.g [[-1.213123, 3.123123], [0.213123, 2.2312123]], so I can load it later.

To render it, I am doing the following

        var geojsonObject = {
            'type': 'FeatureCollection',
            'crs': {
                'type': 'name',
                'properties': {
                    'name': 'EPSG:3857'
                }
            },
            'features': [{
                'type': 'Feature',
                'geometry': {
                    'type': 'LineString',
                    'coordinates': [[-1.213123, 3.123123], [0.213123, 2.2312123]],
                    'freehand': true
                }
            }]
        };

        var vectorSource = new ol.source.Vector({
            features: (new ol.format.GeoJSON()).readFeatures(geojsonObject),
        });

to get my vector source, but it renders as an actual lineString instead of a drawn lineString. Any ideas?


texts: Closing QgsRasterLayer's file handler in Python script

texts:
I'm creating a QgsRasterLayer from temporary file, as result from some processing algorithm.
After that, I no longer can delete file, as expected. But, because It is only a temp file, I need delete It after all my processing is done.
I tried, without success:
   layer = QgsRasterLayer("myTempFile.tif")

#Do some work...

del layer

myTempFile.tif still can't be deleted, but only after closing QGIS.


texts: How to get National Geographic basemap in QGIS?

texts:
I'm proficient in ESRI's products and recently started using QGIS.

I've used the national geographic world map (NatGeo_World_Map) in ArcMap, but am struggling to replicate this in QGIS. The map is also available here.


texts: Randomly sample from geopandas DataFrame in Python

texts:
I am reading a shapefile as geopandas DataFrame and them using pandas subset method to select a region.

geodata = gpd.read_file(bayshp)

geodata.dtypes

geodata.head(10)

OBJECTID    FIPSSTCO    COUNTY  geometry
1   06001   Alameda (POLYGON ((6065941.393835935 2104148.464510527...
2   06013   Contra Costa    (POLYGON ((6143913.640835938 2209458.230510532...
3   06041   Marin   (POLYGON ((5879149.417835938 2203020.920510533...
4   06055   Napa    POLYGON ((6075700.362835937 2441916.530510533,...
5   06075   San Francisco   (POLYGON ((5990480.312835939 2123810.13351053,..

Subset the df:

# Subset based on the index
geosub = geodata.iloc[0:2]

I've got a function that accepts geopandas DataFrame and number of points to sample as arguments.

def sample_random_geo(df, n):

    # Randomly sample geolocation data from defined polygon
    points = np.random.sample(df, n)

    return points

However, the np.random.sample or for that matter any numpy random sampling doesn't support geopandas object type.

I am wondering if there is a way to randomly sample geocoordinates from the spatial region.


texts: Raster Calculation in Qgis using Python Script

texts:
I want to do raster calculation over raster layer for which i am using python script as:

from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
bohLayer = qgis.utils.iface.activeLayer()
entries = []
# Define band1
boh1 = QgsRasterCalculatorEntry()
boh1.ref = 'boh@1'
boh1.raster = bohLayer
boh1.bandNumber = 1
entries.append( boh1 )

# Process calculation with input extent and resolution
calc = QgsRasterCalculator( 'boh@1/10000',
                        'E:/data/abc.tif',
                        'GTiff',
                        bohLayer.extent(),
                        bohLayer.width(),
                        bohLayer.height(),
                        entries )

calc.processCalculation()

Sometimes this script run successfully and gives output but mostly it give error as

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

Is there any other way to do raster calculation over raster layer with python script.


texts: download ESRI World Imagery basemap for offline use?

texts:
So our search and rescue team would like to be able to view ESRI's World Imagery basemap through the free tool (ArcGIS Explorer).  Thing is, we will be offline when we are in the field.  So - is there a way to download our region of the ESRI World Imagery basemap, for use offline?

Note, I'm not a GIS professional, so am probably mixing some terminology.

I found the URL where the World Imagery basemap data lives, and the query form to get back a flat image (jpg I think).  But, we would really like to have it integrated as a basemap so we can use it seamlessly in ArcGIS Explorer.  (Is a basemap a different beast than just a raster layer?)

Is there any way to get this as a 'Data Appliance' from ESRI?  link

Our county staff uses the full paid ArcGIS tools so they might be able to get them for our SAR team's use - the issue is that we would want to view it offline in the free tool on our laptops in the field.

We really do want to use ESRI World Imagery, or an equivalent (it is 1 foot = 1 pixel in our area).

Thanks

UPDATE - 2013-10-01 - SOLVED:

gmapcatcher does the trick perfectly.  It doesn't read from ESRI World Imagery, but it does download tiles from Bing and others (0.3m in our area) for offline viewing.  It's not very feature-rich but that is a big part of the appeal!  Learning curve is just about zero.  There's probably a fairly easy way to turn those downloaded tiles into a correct directory structure for serving a tileset on MapServer, in order to use it with more feature-packed tools.  If anyone finds a quick way to do that tile-file-rearranging please post here, thanks!


texts:
I am trying to load raster in PostGIS from my Windows cmd using:

C:\Program Files\PostgreSQL\9.6\bin>raster2pgsql.exe -s 4326 -I -C -M "D:\9.Taji\data\1_Admin Boundary-20200211T030329Z-001\check\raster\road_raster.tif" -F -T 100x100 public.road  | psql -d postgres -U postgres -p 5432

But every time I run I get this error:

  BEGIN ERROR:  tablespace "100x100" does not exist ERROR:  current
  transaction is aborted, commands ignored until end of transaction
  block ERROR:  current transaction is aborted, commands ignored until
  end of transaction block ERROR:  current transaction is aborted,
  commands ignored until end of transaction block ERROR:  current
  transaction is aborted, commands ignored until end of transaction
  block ROLLBACK ERROR:  relation "public.road" does not exist

How can I fix this issue?


texts: Range of values in slope and aspect changes when I try to unproject my projected rasters

texts: I use the SRTM DEM and I have used the "Project Rater" tool in ArcGIS10.2 in order to project it from GCS_WGS_84 to a certain UTM zone in WGS84. With my projected elevation raster, I have calculated slope (min: 0, max:73 in degrees) and aspect (0-360). When I try to convert again the projected slope and aspect into GCS_WGS_84 (my analysis must be carried out in latitude, longitude) their values change, and I get for slope min:-1,24 max:74 (in degrees), and for aspect I get 1-460! Why does this happen?


texts:
I have a 3D line layer which I have made to a layout and tilted. I want to show how the view coresponds to XYZ and insert a XYZ axis.

Where can I find this in Pro?

Now it just shows overlapping lines. I'm using 2.2.


texts: Custom Styling OpenStreetMap without needing to host it yourself

texts:
I know it's possible to edit the way in which OpenStreetMap display, with CSS-like (like MapCSS) stylesheets.
Is possible to do this without having to host your own tile server?


texts: Point density from XY table in ArcGIS Pro

texts:
I have a large 1.1gb CSV table of points in Lat/Long WGS84 for which I'd like to generate a point density raster. I'm taking the following steps but the process is extremely slow. I welcome any advice on how to improve this or better understand where exactly my bottlenecks are.

In ModelBuilder I have the following steps

XY Table to Points - output to in_memory workspace;
Project - reproject in_memory output from LatLong WGS84 to WGS84 Web
Mercator (Auxiliary Sphere) - here I'm unable to store the output
in_memory (not sure why) so I write it to the gdb;
Point Density - calculate my point density raster from the output of
2.

I'm running on Windows Server 2019 with Intel Xeon Platinum 8175M CPU @2.5GHz (2 processors), 760GB RAM, 64-bit OS which I've spun up an AWS instance with 500GB hard disk.

Currently my ModelBuilder is on step 1 at 1% and it has been there for a while. My CPU usage is at 1% and Memory also at 1%

How can I speed up my process?


texts: Splitting line at point positions using QGIS?

texts:
I am a rooky with QGIS and i am stuck from a basic problems.

I have a polyline and i have some points (2) on the polyline.

I would like to know if there is a tool or script that i could use to split my original polyline into 3 polylines? (Starting point to point A, Point A to B and B to ending point of the polyline)?

I am looking for an automatic process rather splitting the  line manually.


texts: I am wondering if anyone has worked with Pictometry in GIS and their .pmi files.  I recently acquired imagery from them for my area and have discovered that I cannot add the pmi files to the webmap, I presume I would have to export them.  The problem is that there are literally thousands of images and the export would be unrealistic (I'm not very happy with their deliverables, wasn't involved in the acquisition process).  Has anybody worked much with their imagery and have any insights in how to best use it in GIS?


texts: Reclassify is not including one value in its range of one class?

texts:
Trying to reclassify data into a binary map of 0's for no data (in this case, no flooding) and 1's for any level of data (flooding). Range for 1's is about 0-85 feet. and 0's are NO DATA.

But, when I do the reclassification, the output is 3 classes - 0, 1, and 85. The 85's are 4 pixels with a value of 85. I don't understand why this is a separate class, when it is within the range for 1 (0.335534 - 85.22314).

This only happens sometimes, despite my using the same process again and again. I don't have any exclusions listed in the classify box.

Does anyone know why this is happening?


texts:
Using the leaflet minimap plugin to build a map with locator window.

The map starts zoomed at the state level, with minimap hidden on load via this:

var miniMap = new L.Control.MiniMap(layers, { minimized: true, zoomLevelFixed: 9 });

When a user selects a site from the drop-down list, the main map pans/zooms to that site using setView, and the mimimap is unhidden via this:

miniMap._setDisplay(miniMap.options.minimized = false);

However, I'm getting a conflict on initial load, where by the minimap will not pan/zoom to proper place when called like this (main map pans and zooms correctly):

map.setView([40.22087,-74.99483], 18);

Minimap will move to the proper spot, however when initially called in this manner:

map.setZoom(18);
map.setView([39.85019,-75.43627]);

However, calling it that way causes setzoom to break when toggling between points (if map has been panned/zoomed manually before switching points).

You can see the minimap behavior in this JSBin:

https://jsbin.com/lesokutona/edit?js,output

if you select "PH-07" on load. PH-67 uses the second method of setting location/zoom for comparison purposes. (to see the odd behavior when switching points, select PH-07 and zoom/pan, then switch to PH-67.

So there's a conflict here somewhere, but I can't seem to nail it down.

I'm thinking that changing the mimimap's display option is somehow interrupting its location change initially, but I can't seem to figure out how to prevent that (nor does that really explain why calling it via the second method works). I've tried setting it off in it's own function, but that doesn't seem to address the issue.

The second method's setzoom/setview conflict is fairly obvious, but I can't see how that would allow the minimap to plan correctly initially if the other one does not.

I thought maybe it was a bounds issue... but disabling that didn't affect it either.

Anyone have any ideas?


texts: How to build query expression to display text from one-to-many relationship in ArcMap?

texts:
I have two tables in an enterprise geodatabase stored in PostgreSQL and accessed via ArcSDE:

parcel and person
 1 or + parcel can have 1 or + person (relationship is pll_person)

I can display all  parcels with annotation "id"

What I want is to display in each parcel the Name of the person(s) in the related table - see picture please :

I want to display parcels proprietor name  as annotation please any suggestion ?


texts: "Extract by Mask" saving output as temporary file, instead of saving it in specific directory. How to remove this error?

texts:
I am using ArcGIS 10.0, while using ExtractByMask in arcpy output file is being stored as a temporary file, though the same file appears in the designated directory but it's status is temporary, which doesn't open up once we close the ArcMap.

Here is my script

import arcpy
from arcpy.sa import *

shape="G:/Path/To/Shapefile.shp"
raster="G:/Path/To/In_Raster.tif"
output="G:/Path/To/Out_Raster.tif"

out= ExtractByMask(raster, shape)
out.save(output)

How to save output in desired directory?


texts: Is resampling necessary while calculating Normalized Difference Built-up Index (NDBI) on Sentinel 2

texts:
The formula for calculating NDBI using Sentinel-2 bands are as follows:
Band 11 - Band 8 / Band 11 + Band 8

However, bands 11 and 8 have a spatial resolution of 20 metres and 10 metres respectively. Does this mean that I will have to resample the bands such that they are of the same resolution before proceeding with the formula, or can I go ahead with this formula without resampling?
I am using Google Earth Engine for the analysis.


texts:
I am trying to upload shapefile polygons onto my Garmin 60Cx GPS.  I use the Trace tool and load shapefiles via Minnesota DNR's application.  There are almost always gaps in the shapefile on the GPS, making lines and the shapefile hard to read.  Oddly enough, these gaps are often where I have not lifted the mouse/clicked anything to offset the trace tool.

Is there something that I am doing wrong?


texts: Spatio Temporal Interpolation with gstat in R

texts:
I am very new to spatial evaluation and come from psychology.
I am using the software R and the packages "gstat" and "spacetime".
I would like to do a spatio temporal interpolation. For this I follow the paper of Gräler et al. (https://cran.r-project.org/web/packages/gstat/vignettes/spatio-temporal-kriging.pdf)
Unfortunately I can't find/fit the right variogram model.  I can create the empirical variogram and this is also conclusive to me, but then I do not get any further. I do not understand how to define the individual parameters such as "sill" or "nugget" or what they stand for.
Here are my previous approaches:
My ST-Dataframe:
> df.stf
An object of class "STFDF"
Slot "data":
# A tibble: 2,406 x 6
     KRS month DeprIndex MMW10 MMW25 NewInfect
   <dbl> <int>     <dbl> <dbl> <dbl>     <dbl>
 1  1001     1    -4.08   NA   NA           NA
 2  1002     1    -2.28   13.3 NA           NA
 3  1003     1    -3.29   17.8 11.7         NA
 4  1004     1    -4.31   NA   NA           NA
 5  1051     1    -2.51   17.7 NA           NA
 6  1053     1    -1.07   13.0  9.93        NA
 7  1054     1    -0.863  NA   NA           NA
 8  1055     1    -1.63   NA   NA           NA
 9  1056     1     0.887  NA   NA           NA
10  1057     1    -1.21   14.7  8.84        NA
# ... with 2,396 more rows

Slot "sp":
class       : SpatialPolygonsDataFrame
features    : 401
extent      : 3280359, 3921536, 5237511, 6103443  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +units=m +no_defs
variables   : 1
names       :   AGS
min values  :  1001
max values  : 16077

Slot "time":
           timeIndex
0001-01-30         1
0002-01-30         2
0003-01-30         3
0004-01-30         4
0005-01-30         5
0006-01-30         6

Slot "endTime":
[1] "0002-01-30 CET" "0003-01-30 CET" "0004-01-30 CET" "0005-01-30 CET" "0006-01-30 CET" "0007-01-30 CET"

So, I have 401 german counties and the mean monthly PM10 value over six month.
My empirical variogram:
eVgmPm10<-variogramST(MMW10~1,df.stf,tlags = 0:5)

eVgmPm10$dist<-eVgmPm10$dist/1000
eVgmPm10$avgDist  <- eVgmPm10$avgDist/1000
eVgmPm10$spacelag <- eVgmPm10$spacelag/1000

plot(eVgmPm10, map=F)

So far everything works as it should (I guess)
Now I wanted to make as in the paper explained a seperable model:
separableModel <- vgmST("separable",
                        space = vgm(0.9, "Exp", 200, 0.1),
                        time = vgm(0.9, "Sph", 3.5, 0.1),
                        sill = 124)

I understand, that I have to create a variogram for space and time. But I do not know how to define the parameters correctly (psill, model, range, nugget) for the variograms. If I use the same parameters as in the paper I get folowing plot:

As you can see, in this model, I just have two lags and the variogram looks realy strange. So, I think it is because of the wrong parameter for my model. I also tested another approach like in the paper.
sumMetricModel <- vgmST("sumMetric",
                        space = vgm(20, "Sph", 150, 1),
                        time = vgm(10, "Exp", 2, 0.5),
                        joint = vgm(80, "Sph", 1500, 2.5),
                        stAni = 120)

Again, I do not have any idea and clue how to set the parameters so I took the same as in the paper and plot all three.

I am pretty sure it depend on the parameters but I do really not know how to find the correct variogram model to fit.


texts: QGIS Attribute tables and adding fields data to all attribute tables through a lookup from an imported database attribute table

texts:
I understand you can add data to a single attribute table in QGIS via the field calculator from a database CSV layer. This allows a lot to be done with a single attribute table. Now this is good if you have added all your features to a single layer. However I would like to add my data before I merge so that if the data changes I don't have to merge and repeat this process.

Is there any way to set the field calculator on more than one attribute table and alter data on those tables depending on the data in them?

I am talking about a power tool that allows you to alter data in all attribute tables based on a look up to a existing table. Has anyone found anything to do this?

I have a database in QGIS which was prepared from an external source . It sits in QGIS nicely. I wanted to access the data from each each feature that has its own fields eventually to put them in one layer, so I can fill in information from the database. I can do this on each layer or feature but wondered if there was a global way to set it on all features without visiting each feature. I have seen a YouTube aid which shows you how to do this on a single layer with multiple features in it but not many single features.


texts: Postgresql table won't display in QGIS 2.2

texts:
I just downloaded QGIS 2.2 to check it out and now one (just the one) of the tables from one of my PostgreSQL database won't display.

It loads just fine in 2.0.  The records are all there (there's only 3 in this case), I can see them in the attribute table, I can query them, but they don't show up on the map.

When I try to zoom to a particular feature it pans the map down to the origin of the datum (coordinates 0,0) - but when I zoom to the extent of the layer it sets the map back to where it is supposed to be.

What am I missing in the jump from 2.0 to 2.2?


texts: Render all from react component into .setHTML() Mapbox-gl-js with React

texts:
How can I Render react component inside .setHTML() whe in create popup?
import PopUpClick from "../../components/Map/PopUpClick";

map.on('click', 'mapPins', function(e) {
            map.getCanvas().style.cursor = 'pointer';
            const coordinates = e.features[0].geometry.coordinates.slice();
            const name = e.features[0].properties.title;
            const img = e.features[0].properties.img_url;
            const id = e.features[0].properties.pin_id;
            // const popUpHTML = ReactDOMServer.renderToString(
            //     <PopUpClick
            //         name={name}
            //         img={img}
            //         id={id}
            //     />
            // );

            const popUpHTML = (<PopUpClick
                    name={name}
                    img={img}
                    id={id}
                />
            );

            new mapboxgl.Popup({
                closeButton: false,
                offset: 16,
                anchor: 'top'
            }).setLngLat(coordinates)
            .setHTML(popUpHTML)
            .addTo(map);
        });

When I click on symbol marker I got [object Object] in my div
I need all functionality from react component, and with ReactDOMServer.renderToString i got only pure HTML.
I need to have all react stuff ass Redux and click listeners.


texts: Why learn/use Python Toolboxes over Python Script Tools?

texts:
I've written a few Python Toolboxes (which are new at ArcGIS 10.1), but am yet to decide whether/when I should write them rather than Python Script Tools in a standard toolbox.

I thought the Online Help might enlighten me when it prefaces some dot points with:

  Once created, tools in a Python toolbox provide many advantages

However, the five advantages listed all seem to be over not being able to use Python to write tools, and none seem to specify an advantage of Python Toolboxes over Python Script Tools.

The two advantages that I can think of are:

I can now write a "pure" Python tool in a single Python script without having to hook it up to a separately authored dialog with its Tool Validation looking like it was tacked on but I'm happy to be pragmatic rather than pure in this regard
I could now use code (Python or any language capable of writing text files) to automate the writing of Python toolboxes but I am yet to come across a requirement to do this

Am I overlooking the compelling case that led Esri to provide the Python Toolbox capability and, if so, what is it?


texts: Writing PyQGIS code to select items in line layer by length?

texts:
I have been trying to write a short code to select items in a line layer that are longer than 2000m, but I don't have much luck at the moment now.

layer=iface.activeLayer()
selected_fid=[]
lenth=[]
for f in layer.getFeatures():
    selected_fid.append(f.id())
    lenth.append( f.geometry().length())

for i in range(0,len(lenth)):
    if lenth[i]>2000:
        print i

This is what I have, and the first thing I notice is that the len function does not work in QGIS python console, I am not sure why. It gives me an error as below:

len(lenth)
Traceback (most recent call last):
  File "<input>", line 1, in <module>
TypeError: 'list' object is not callable

Additionally, I am not sure how to set those items using setSelectedFeatures function.

Any suggestions?


texts: Calculating daily max and min temperature data from GLDAS over several years using Google Earth Engine

texts:
I am still new to the Google Earth Engine and try to calculate the daily max and min temperature based on the GLDAS.
Could you advise me on how I correct this code to calculate daily max and min temperature?
var startDate = ee.Date('2018-01-01');
var endDate = ee.Date('2021-01-01');

// Specify Country names
var worldcountries = ee.FeatureCollection('FAO/GAUL/2015/level2');

var filterCountry = ee.Filter.eq('ADM2_NAME', 'Kurgantepa district');

var country = worldcountries.filter(filterCountry);

Map.addLayer(country);
Map.centerObject(country, 10);

var dataset = ee.ImageCollection("NASA/GLDAS/V021/NOAH/G025/T3H")
  .select('Tair_f_inst')
  .filterDate(startDate, endDate);

var list_dataset = dataset.toList(dataset.size());

print(list_dataset);

var getTemperature = function(image) {

  var value_temperat = ee.Image(image)
    .reduceRegion(ee.Reducer.first(), country)
    .get('Tair_f_inst');

var temperat_K = ee.Number(value_temperat);

  return temperat_K;
};

var count = dataset.size();

var temperat_list = dataset.toList(count).map(getTemperature);

print("Tair_f_inst list", temperat_list);

var allDates = ee.List(dataset.aggregate_array('system:time_start'));

var allDatesSimple = allDates.map(function(date){
  return ee.Date(date).format().slice(0,10);
  });

print(allDatesSimple);

var len = temperat_list.size();

print(len);

var paired = allDatesSimple.zip(temperat_list);

print (paired);

var title = {
  title: 'Daily Temperature',
  hAxis: {title: 'Time'},
  vAxis: {title: 'Temperature (K)'},
};

var chartDaily = ui.Chart.image.seriesByRegion({
  imageCollection: dataset,
  regions: country,
  reducer: ee.Reducer.mean(),
  band: 'Tair_f_inst',
  scale: 5000,
  xProperty: 'system:time_start',
  seriesProperty: 'SITE'
}).setOptions(title)
  .setChartType('ColumnChart');

print(chartDaily);

var myFeatures = ee.FeatureCollection(paired.map(function(el){
  el = ee.List(el); // cast every element of the list
  var geom = country;
  return ee.Feature(null, {
    'date': ee.String(el.get(0)),
    'value':ee.Number(el.get(1))
  });
}));

//print(myFeatures);

// Export features, specifying corresponding names.
Export.table.toDrive(myFeatures,
"Temperature", //my task
"GEE_Folder", //my export folder
"daily_temperat",  //file name
"CSV");


texts: GeoServer: Get styles from PostgreSQL/PostGIS table

texts:
We're running GeoServer with all layers coming from a PostGIS data source. Most people work directly on the PostGIS database, but we have one person sitting 10.000km away, and for him, connecting to PostGIS via a VPN is almost impossible because the latency/connection quality is so bad.

So I created an account for him in GeoServer, under which he works on the PostGIS base via QGIS+WFS. This works perfectly for him (less overhead) basically.

This system has one little flaw: With QGIS, we created a table containing all the SLD data for the layers in the PostGIS database. QGIS seems to manage this in its own way, but it works perfectly when working with QGIS directly on the PostGIS database.
But: That one user working on the data via WFS/GeoServer cannot see these styles. The style management in GeoServer, as far as I'm concerned, is a "closed" system.

Yet, is there any trick whatsoever to make the styles from our PostGIS styles table available in GeoServer, most appreciatedly respecting the relations between the styles and the layers?

The table definition is as follows:

CREATE TABLE public.layer_styles
(
  id integer NOT NULL DEFAULT nextval('layer_styles_id_seq'::regclass),
  f_table_catalog character varying(256),
  f_table_schema character varying(256),
  f_table_name character varying(256),
  f_geometry_column character varying(256),
  stylename character varying(30),
  styleqml xml,
  stylesld xml,
  useasdefault boolean,
  description text,
  owner character varying(30),
  ui xml,
  update_time timestamp without time zone DEFAULT now(),
  CONSTRAINT layer_styles_pkey PRIMARY KEY (id)
)

The f_table_name column is a reference to the PostGIS table/layer to which the style applies. stylename is just a unique, freely chosen name. styleqml contains some QGIS-specific information, and stylesld contains the raw SLD code, again including the layers the style applies to.

So if I'm taken correctly, I could, as well, copy all the SLD fields and create them as styles in GeoServer. But since there are some changes sometimes, it would be absolutely great if there was a way to have GeoServer's styles synchronized with this table.

Any ideas?


texts: leaflet on android app with GPS localisation

texts:
I put my leaflet html code on android webview component and it works awsome like native app of maps ,but the GPS location doesn't work .
mycode

<!DOCTYPE html>

<title>GIS</title>
<meta http-equiv='content-Type' content='text/html; charset=UTF-8' />
<link rel="stylesheet" href="https://unpkg.com/leaflet@1.2.0/dist/leaflet.css" />
<link rel="stylesheet" href="https://domoritz.github.io/leaflet-locatecontrol/dist/L.Control.Locate.min.css"/>


texts: Select point with highest value within a buffer

texts:
I work with ArcMap 10.2.1 and I created buffers around points. Some buffers contain many points and others just one. I would like to select one point for each buffer, but the points can represent different attributes having different values.

So when I have for example one buffer with two points, one point could be organism one and the second point could be organism two. Organism one has a value 30% and organism two has value 70%. I would like to select organism two and delete organism one from the buffer (file).
Only thing is I have many buffers, with varying amount of points (organisms) and different percentages..So I would like to do it for all the buffers..
I am using modelbuilder as I do not know python.
Does anybody have a suggestion? That would be great!


texts: Getting attribute from intersecting polylines using ArcGIS Desktop?

texts:
I am using ArcMap 10.2.2

I've got 1 polyline feature class, which includes aprox. 500 lines with attributes (streets).
To fill out certain documentation, I need to know which streets cross each of the streets. Those. I need the names of the streets crossing it to appear in the attributes of each line. Maybe a comma. Or if there is some other way.


texts: SQL within Python Script to get only part of the expression by GetParameterAsText

texts:
I was wondering if anybody could help me.
Im trying to create a tool within ArcMap using a python script.
Within the script i've been trying figure out how to change it so that I  only need to get the numerical value from GetParameterAsText. So that for the tool i do not need to include the whole SQL expression just have the user input a value.

This the last method that I have used

wf_value = arcpy.GetParameterAsText(4)

wf_expression = ' "AGE_18_64"' + ">=" + "'%'" %wf_value

arcpy.Select_analysis(counties, workforce, wf_expression)

ERROR MESSAGE:
ExecuteError: ERROR 000358: Invalid expression 3000
Failed to execute (Select)

As I'm guessing this way doesn't work as the wf_value is a float not a string but I can't seem to get it to work whichever way I've tried.
Any help or suggestions would be great!
Thanks!


texts: Generating DEM from isolines and spot heights using GRASS?

texts:
I have two shape files containing height information of a certain area and I need to create a raster DEM with GRASS using them.

file 1 consists of height information in contours (isolines)
file 2 consists of height information with spot heights

What would be the best technique to generate raster DEM with these two types of data?

Do i need to generate some spot heights from contours and then integerate with existing spot heights for DEM generation or any other method?


texts: Location of Index file for geoserver on Openshift

texts: I've configured an openshift application to run Geoserver with a postgis database.  Where in the open shift directory should i place my index.html?


texts: OpenLayers 3 support for WMS 1.3.0 CRS:1

texts:
I am trying to get a TileWMS layer working with the Coordinate System as CRS:1 defined in WMS 1.3.0. By defining the right projection stated in CRS:1, axisorientation as 'es', I could get it working for ImageWMS. However, when using TileWMS, I found all the tileGrids defined in OpenLayers 3 are always aligned with 'en' direction. tileindex x increments as x goes right, tileindex y increments as y goes up. So it is different from the CRS:1 definition. This causes the tiles are put into wrong direction in y axis.

My question is that is there any way that I can change tiling direction to follow the axixorientation defined in my projection, in this case, is 'es'?

Or is there any advice or examples for 'es' projection like CRS:1 in tiled API of OpenLayers 3?


texts: Further compression of LiDAR data besides LASZip

texts:
I am looking for a way to compress files other than LASZip that may make the files a bit more smaller, is there an efficient way to do this?

Would zipping up a LAZ file further compress the file?  Or is it unnecessarily redundant.

I am working with a large amount of data that I have zipped into LAZ format, but was wondering if there were further steps I could take to reduce the size even more.


texts: How to use search and update cursor on two tables to update values?

texts:
I would like to read the values from the company field of linefc, if the values from this field in table then a the WORSTOFFENDER field should be updated with the value YES, if it is not, then NO.  What needs to be added to complete this?

with arcpy.da.SearchCursor(linefc,["Company"]) as cursor:
    for row in cursor:
        linecur =row[0]

with arcpy.da.SearchCursor(table,["DAMAGING_PARTY"]) as tablecursor:
    for tablerow in tablecursor:
        tablecur =  tablerow[0]

print "DONE"


texts: I'm using ArcMap, I have a shapefile which is my wastewater pipelines, now I need to merge some rasters which are about ground movement and groundwater table with this shape file to know what is the ground situation of each pipe. How can I merge all these together? Should I change pipelines to points and then merge them to other rasters?


texts: ArcGIS Modelbuilder model clogs up memory?

texts:
I have a model in which a zonal histogram is computed of each line in a feature class, showing the amount of cells each line intersects in a raster file. After this, this column is joined to a table that has values of the first line in it.

This works. The problem is, that as can be read here, intermediate results are stored on the RAM (if I understand correctly). In my case, the first ~100 iterations are very fast (<0.1 seconds). However, after that it becomes increasingly slow, and since I have 50k lines in total, I need some way to speed this up.

What is the easiest way of running this model elsewhere, where it does not fill up the RAM?


texts: Go through QGIS Python for loop through features in reverse order?

texts:
I am trying to go through a for loop in QGIS's Python module in reverse order for my features.

from qgis.utils import iface
from PyQt4.QtCore import *
layers = iface.legendInterface().layers()
for layer in layers:
    name = layer.name()
    myLayer = QgsMapLayerRegistry.instance().mapLayersByName( name )[0]
    if name.startswith('Test'):
        for f in layer.getFeatures():

I would like to be able to go through that second for loop in reverse order. I know python has the reversed() function but that is only for lists. I can't seem to use it for layer.getFeatures().

Is there a way to do this?


texts: Installing gdal to geoserver running with apache tomcat

texts:
I have been following Mastering GeoServer pdf from Colin Henderson to install gdal but didn't succeed. Here is what I did:

install Java, I've installed JRE windowsx64 (because my pc is x64)
install apache tomcat, I've decided to go with the version 8.5.13 instead of version 9 (because it would be more reliable) then I've installed apache-tomcat 8.5.13 via 32-bit/64-bit Windows Service Installer. While installing I've ticked service startup and native and unticked documentation, host manager and example. I've configured tomcat, in Java section I've added following parameters ;

-Dserver
-Xmx1024m
-Xms512m
-XX:MaxPermSize=128m

locahost:8080 is working now.

deploying geoserver to apache tomcat, I've downloaded GeoServer 2.11.0, war file from Geoserver's website.

localhost:8080/geoserserver/web Works now.

Installing the GDAL binary libraries;
From gisinternals, I have tried many of them over and over, none of them worked. For example I've downloaded complete versions of
(I have c++2012 and c++2013 os my system)
then set the environment variables as needed. Then applied this "To perform the installation, we just need to extract the .jar files contained within the archive into the WEB-INF/lib folder of our GeoServer installation."

it didn't work.

What's wrong ?


texts:
I have a website where an OpenLayers 3 map fills part of a page. The site is intended for mobile viewing, and my problem is that I'm not able to scroll the page by pushing and dragging on the map.

I've tried disabling all interactions, and found that setting interactions: OpenLayers.interaction.defaults({mouseWheelZoom: false}) allowed me to use mouse scroll when in developer mode in Chrome.

I've adjusted the controls and interactions for the map according to my needs, but this integration with the rest of the page is still a problem.

I know that it's possible to get this effect with OpenLayers 2.

Do any of you know a way to enable this? I guess the event is called touchmove.


texts: Summing column values in table using arcgisscripting?

texts:
I need to create a python script that sums a column's values in a table. Thereafter, the script is supposed to get the percentage of each individual row in that column and populate this percentage values in a another column. I been struggling for quite some time now, I checked the net, but there is nothing of this sort. The nearest is an ArcGIS Discussion Forum thread, but its for an older version of ArcMap, or I am using it wrong.

import arcgisscripting, sys
gp = arcgisscripting.create()

intable = sys.argv[1]
field = sys.argv[2]

# Create search cursor
rows = gp.SearchCursor(intable)
row = rows.Next()
x = 0.0
# Enter while loop for each feature/row
while row:
    x += row.getvalue(field)
    print x
    row = rows.next()
#note value can be rounded depending on field type
gp.calculatefield(intable,sys.argv[3],float(x),"PYTHON")


texts: What are some examples of model-based upscaling of raster data?

texts:
Remotely-sensed raster data may be too fine-grained for many purposes.

A typical approach to coarsen the grain is to upscale by reducing the resolution.  If the data are continuous (e.g elevation, surface temperature) cells on the raster can be aggregated (and therefore reduced in number) by using a mean function.  Or if they are discrete (e.g. vegetation) they may be aggregated using the mode of adjacent cells.

Moving windows are a second example of upscaling.  Here a "window" is slid across the raster and the cells under the window are replaced using some function of their values (e.g. using a mean or mode).

I believe that model-based approaches to upscaling will simplify a raster by adopting a more complex model of the raster (i.e. the cells and their adjacencies).

Can anyone clarify my definition of model-based approaches, and provide an example?
Do the resulting "cells" even need to be regular, could they vary in shape and size?

EDIT:  To be more specific about what I mean by model-based upscaling.  Here is an example of what I think may be considered a model-based approach.

Consider using a mathematical graph to represent discrete features on a raster.  Then use some method to identify components of these features (i.e. groups of vertices).  A Voronoi tessellation, for example, using the components as generators, could then be used to produce an irregular tiling that represents an upscaling of the original raster.

... I could give more detail here, but opt not to, as I am interested in mostly whether or not this would be considered a model-based approach to upscaling ... or alternatively how such operations could be described in generic terms.  Is this merely a way of doing interpolation (as defined in an answer below)?


texts: Is it possible to create Thematic Maps using two columns in Mapinfo?

texts:
I have created Thematic map with one column in Mapinfo. Is it possible to create Thematic Map using two columns in Mapinfo?

For Example:

  I have a column with name of place and a column with date. Is it
  possible to combine name of place and date to generate thematic map.


texts:
I have an excel file with about data from 600 locations, each location having data for 1200 months(12 months* 100 years). I need to interpolate all point locations into rasters for each month of each year separately. I understand I need a script to automate the process because it would be a huge number of rasters. Can anyone recommend me a similar script?
I am absolutely new to GIS.


texts: Disable editing after commit in Leaflet.Editable

texts:
I'm using Leaflet.Editable to draw polygon on my map, but I don't want to edit the polygon that I'm drawing. I want it to be fixed on my map.
How can I do that?


texts: Infobox not fixed to marker when Google Map is wrapped (on zoom 0,1,2)

texts:
I use a custom infobox that is fixed to the marker's position. Everything is good (see screenshot 1) until the map reaches the zoom 3 or lower, precisely, when the map is wrapped. The infobox is no more fixed to the marker and it's moved some pixels to the left of the marker (see screenshot 2).

Screenshot 1

Screenshot 2

I use the following code to convert the marker's position (LatLng) to "XY" points (See the code) and to fix the custom infobox to the marker.

var scale = Math.pow(2, map.getZoom());

var nw = new google.maps.LatLng(
    map.getBounds().getNorthEast().lat(),
    map.getBounds().getSouthWest().lng()
);

var worldCoordinateNW = map.getProjection().fromLatLngToPoint(nw);

/**
 * [@marker_position] - The position (LatLng) of the marker */

var worldCoordinate = map.getProjection().fromLatLngToPoint(marker_position);

var pixelOffset = new google.maps.Point(
    Math.floor((worldCoordinate.x - worldCoordinateNW.x) * scale),
    Math.floor((worldCoordinate.y - worldCoordinateNW.y) * scale)
);

/**
 * [@mapPosition] - This will get the position of the map on the screen */

var mapPosition = jQuery('div.mapContainer').position();

var infobox_half_width = jQuery('div.infobox').width() / 2;
var margin_top = marker_icon_height + jQuery('div.infobox').height();

/**
 * Set the position of the infobox on the map.
 * [@left] - The position of the marker (Horizontaly) ACCORING to the map's left position on the screen.
 * [@top] - The position of the marker (Verticaly) ACCORING to the map's top position on the screen.
 * [@margin-left] - Move the infobox 50% of its size to the left.
 * [@margin-top] - Move the map to the top of the marker (According to the marker image size). */

jQuery('div.infobox').css({'left':(pixelOffset.x + mapPosition.left + 'px'),
                 'top':(pixelOffset.y + mapPosition.top + 'px'),
                 'margin-left':('-' + infobox_half_width + 'px'),
                 'margin-top':('-'+margin_top+'px')
               }).fadeIn('slow');

What i want is to keep the infobx fixed to its marker even when the map is wrapped.

I think the problem is in this part of the code (see code bellow) but i don't know how to adapt it to maintain the infobox fixed to its marker when the map is wrapped. The code needs something (in the calculations) to tell the map to calculate the infobox's position related to the bounds of the completely visible part of the world. Does anyone have an idea on how to resolve this issue?

var scale = Math.pow(2, map.getZoom());

var nw = new google.maps.LatLng(
    map.getBounds().getNorthEast().lat(),
    map.getBounds().getSouthWest().lng()
);

var worldCoordinateNW = map.getProjection().fromLatLngToPoint(nw);

/**
 * [@marker_position] - The position (LatLng) of the marker */

var worldCoordinate = map.getProjection().fromLatLngToPoint(marker_position);

var pixelOffset = new google.maps.Point(
    Math.floor((worldCoordinate.x - worldCoordinateNW.x) * scale),
    Math.floor((worldCoordinate.y - worldCoordinateNW.y) * scale)
);

Note: I didn't want to use the GMaps overlay to create the infobox just to prevent the slow loading of the map when there's too many markers. I've tried it already with 600+ markers (markers + overlays displayed by default) and it takes too much time to load. Creating the custom infoboxes outside of the map and using that code to fix each infobox to its marker had a positive impact on the loading speed.


texts: Updating Polygon only on multiple selection, not single, using ArcPy?

texts:
I have a script that selects all addresses within a polygon (also multi polygons), and then sums a column "Anzahl" from those selected addresses. The Sum is then written to the Polygon in the column Anzahl.

Problem is, it only works on a single selected polygon. If I select multiple polygons it breaks line 33:

line 33, in <module>
    edit.stopEditing(True)
RuntimeError: start edit session

How can I change this start/stop edition to let it work on multiple selected polygons.

The goals is for each polygon select the addresses inside, sum the column "Anzahl" of the selected addresses. And write the sum, to that polygon.

import arcpy

# Get Parameters
polygon_layer = arcpy.GetParameterAsText(0)
address_layer = arcpy.GetParameterAsText(1)

selected_oids = []
with arcpy.da.SearchCursor(polygon_layer, "OID@") as cursor:
    for row in  cursor:
        selected_oids.append(row[0])

workspace = arcpy.Describe(polygon_layer).path

anzahl_total = 0
try:
    edit = arcpy.da.Editor(workspace)
    edit.startEditing(False, True)
    for oid in selected_oids:
        arcpy.SelectLayerByAttribute_management(polygon_layer, "NEW_SELECTION", "OBJECTID = {0}".format(oid))
        arcpy.SelectLayerByLocation_management(address_layer, "WITHIN", polygon_layer)
        #address_count_within = int(arcpy.GetCount_management(address_layer).getOutput(0))
        #arcpy.AddMessage('Number of address points within Polygon (OBJECTID={0}): {1}'.format(oid, address_count_within))
        with arcpy.da.SearchCursor(address_layer, "Anzahl") as cursor:
            for row in cursor:
                anzahl_total = anzahl_total + row[0]
        #print anzahl_total
            #edit.startOperation()
            with arcpy.da.UpdateCursor(polygon_layer, ["anzahl"]) as cursor:
                for row in cursor:
                    row[0] = anzahl_total
                    cursor.updateRow(row)
            #edit.stopOperation()
    edit.stopEditing(True)

finally:
    pass
    # qry = "OBJECTID IN ({0})".format(','.join([str(o) for o in selected_oids]))
    # arcpy.SelectLayerByAttribute_management(polygon_layer, "NEW_SELECTION", qry)


texts: Where can I find a shapefile of the UK electricity grid network?

texts: Does anyone know where can I download a shapefile with the UK electricity grid network?


texts: Converting QGIS X,Y info to lat, long using Python

texts:
I am trying to convert QGIS x,y points into the lat, long coordinates. I have (335687.7471962841, 4487822.535352517) supposedly point to somewhere in the Chicago metropolitan area. I tried the accepted answer here (without inits) and got (37.34869549513416 3.0155343399807206) which is somewhere in the Mediterranean Sea. Where am I doing wrong?

Here is what I have done:

from pyproj import Proj, transform

inProj = Proj('epsg:3857')
outProj = Proj('epsg:4326')
x1,y1 = 335687.7471962841,4487822.535352517
x2,y2 = transform(inProj,outProj,x1,y1)
print (x2,y2)

Output:
37.34869549513416 3.0155343399807206


texts: ogr2ogr difference between -a_srs and -s_srs

texts:
What's the difference between the two flags -a_srs and -s_srs.
From what I read they both define the EPSG the Input-File has. When using -s_srs ogr2ogr says that -t_srs must be used as well.
But when setting -a_srs instead of -s_srs it runs without -t_srs.


texts: Extract paleoclimate variables from wordclim

texts:
I1m relative new in to RStudio and I'm trying to download the paleoclimate variables (http://ww.worldclim.com/paleo-climate1) from worldclim, but I can't, the unique solution that I found was to download all the rasters on a folder and to tell rstudio the folder that contains the files. I tried with the common comands for getdata from raster package:
getData('worldclim', var='bio', res=10)

But I need specific scenarios and models (eg. LGM from CCSM4) but the details of that comand explains that not have more arguments apart of variable and resolution:
if name='worldclim' you must also provide arguments var, and a resolution res. Valid variables names are 'tmin', 'tmax', 'prec' and 'bio'. Valid resolutions are 0.5, 2.5, 5, and 10 (minutes of a degree). In the case of res=0.5, you must also provide a lon and lat argument for a tile; for the lower resolutions global data will be downloaded. In all cases there are 12 (monthly) files for each variable except for 'bio' which contains 19 files.
getData('worldclim', var='tmin', res=0.5, lon=5, lat=45)
getData('worldclim', var='bio', res=10)


texts: QGIS transparent polygon returned after digitizing

texts:
After the digitalization of my raster map as per below:

I saved it in .geojson format and .kml format.

The geoJSON works fine:

but once I export as a kml and open (in Google MyMaps), then I have borders only:

In the Google Earth it looks the same.

Is there a step, which I forgot during my digitalization process?


texts: Cokriging using categorical covariables in ArcMap?

texts:
I am using cokriging to model a continuous variable across my study unit (soil carbon). I have one continuous covariable (soil moisture) but also a categorical covariable (vegetation community). I am using the Geostatistical Analyst Wizard to do this. (I ran linear models in R to establish the significant covariables and any transformations required).

Am I right in thinking that I can't include the categorical covariable directly into the model? If so, can anyone please suggest a way that I can include this information within a cokriging model?

I am using cokriging as it is enabling me to compare the cross validation statistics with and without the covariables. But open to suggestions of other methods which allow inclusion of covariables and provide some idea of model fit, that I can apply in ArcMap.


texts: Using CQL for a poor man's temporal filter on a server offering only WFS 1.1.0

texts:
There's a WFS I'm accessing via script to pull in features that will match a particular day.  The attribute I care about in the output is DATE_TIME.  However because it's WFS 1.1.0 it does not support temporal CQL queries (per the GeoServer documentation anyway).

The DATE_TIME values are stored in ISO format, YYYY-MM-DDTHH:mm:ss.  I would like to issue a CQL query that will return only rows matching a wild-carded day or month, but have been unable to find examples of this kind of thing that will work with 1.1.0:

/ogc/wfsServer?service=WFS&request=GetFeature&typename=[name]&version=1.1.0&propertyname=DATETIME,EVENT,DESCRIPTION,ETC&[put CQL here to do ONLY RETURN ROWS WHERE DATETIME=2017-01-01%]&maxFeatures=1000&outputFormat=application/json

Not looking for my homework to be done for me - If there is a reference page with examples on this, would be happy for a link, as I haven't turned up anything by googling.


texts: How to do Kernel density analysis in QGIS 1.7?

texts: Do someone have some suggestion? We would like to tray Kernel Density analysis as we did with ArcGis using QGIS?


texts: Using ArcPy to fetch Information from File Geodatabase

texts:
I am new to GIS.
I have a .gdb folder with data and I need to process that data to find educational institutions which are within say X Km of the Roads and within Y Km of any other educational institution.
I need to generate a shapefile with the resulting data.
Do I need to install ArcGIS to do this or is it possible to do this in ArcPy alone?
I am pretty OK with Python and I seek a Python Snippet.  I am not looking for full code but just an idea to get me started.


texts:
I need to reinstall GRASS 6.4 on this Ubuntu computer due to compatibility issues. I had installed GRASS 6.4 a few weeks ago to work with a simulation program. I used the command

apt-get install grass

as explained at the website https://grass.osgeo.org/download/software/linux/. Later I installed GRASS 7.1 to try to start making compatibility changes to the simulation. After doing this, the simulation no longer worked. I decided to try uninstalling all GRASS and reinstalling GRASS 6.4.

However, when I use the same command at the terminal, it only installs GRASS 7.0. I tried using a repository for GRASS 6, but that installed GRASS 7.1.

What is the simplest method I can use to reinstall GRASS 6.4?


texts:
I have looked and do not see this answered in any other discussions.

Purpose: Site suitability analysis for brewers spent grain processing facility

Method:  Locate intersecting buffers with greatest amount of spent grain - need field to sum barrels based on intersections (ex. Mountain Brewing Co. produces 3,000 barrels a year and Sierra Nevada produces 800,000 barrels a year - if the radii from the two breweries intersects the intersection would have a sum of 803,000 a year)

Using ArcMap 10.1, I have one shapefile that I have a 15 mile radius around each point (brewery).  I am trying to identify all intersecting buffers by summing a specific field.  I would eventually like for the intersections with the greatest number to be a 3D representation of barrels per year based on the sum field (I don't think I can do this in ArcMap 10.1 but need to transfer everything over to ArcScene, but I'm hoping not).


texts: Postgis: Is there a way to do an intersection query based on lines while taking into account the dateline?

texts:
This is an exact copy of a question I asked on Stackoverflow a few weeks back and still haven't found a good answer for it. The original question is from here: original question.

I've been searching about this for a while and can't seem to find anything useful.

In my database, I have an outline of the world coastline saved as lines. I want to take another line and determine where on the coastline that line intersects. My current query is like this:

SELECT st_intersection(makeline(p1.latlong, p2.latlong)::geography, the_geom::geography)
FROM worldcoastline, points p1, points p2
WHERE p1.id = 1 AND p2.id = 2

All is well and good, and everything seems to line up when I test it in KML and such, except when the dateline is involved. Instead of going across the dateline, it will go all the way around the world, and show me all these collisions it couldn't have had (across Russia, Europe, America). The most interesting thing I find though is if I throw in a st_intersects on the makeline as a geography, it returns false. I'm considering adding it to the where clause, but the problem is lines that cross the dateline and intersect land will return really whacky results anyways!

Is there anything I can do in PostGIS to fix this problem?

EDIT: I just found this in the reference:

  ST_Intersection — (T) Returns a geometry that represents the shared
  portion of geomA and geomB. The geography implementation does a
  transform to geometry to do the intersection and then transform back
  to WGS84.

So I guess that's why it's having dateline problems. Any alternatives or suggestions?

EDIT: I just wanted to point out that I did add the st_instersects WHERE clause, but of course if a line does intersect the dateline as well as a coastline I will end up with a messed up ST_Intersection because it will still try to find the intersections going the long way around the world instead of across the dateline


texts:
I have a similar problem to Getting back to side bar after detaching layer panel in QGIS?, but with toolbars.

So I accidentally detached digitizing toolbar (and other toolbar as well), and cannot put it back. I tried double click, untick and tick it again in Menu View - Toolbar, hovering it for several seconds over the empty toolbar place, none works.

Anybody know how?


texts: Kriging values greater than the most extreme value of the input layer?

texts:
I'll admit I don't understand kriging to the fullest. Yet, I've read that the weights should add up to 1 as in IDW and Natural Neighbors.

However, surrounding the biggest elevation points of my vector layer, I get some raster cells with greater values.


texts: Can osm2po produce a travel distance buffer for a given time?

texts: Given a starting point and a travel time, I would like to produce a travel distance polygon. Can I request this from osm2po via http?


texts: Using arcpy.MakeFeatureLayer_management on Shapefile and Feature Class

texts:
The larger goal with this code is to create a temporary polygon based on the calculated coordinates (which it does), then to see if the polygon intersects with the "unobscured_FC". After I run the select by location the polygon will be deleted, and all sorts of other hilarity will ensue. I only need the temporary polygon for the select.

Now, I've read that you can't run arcpy.SelectLayerByLocation_management() on shapefiles or feature classes, which is why I tried to make feature layers from the two. However, despite that, I'm getting this error:

  ExecuteError: Failed to execute. Parameters are not valid.
  The value cannot be a feature class
  ERROR 000840: The value is not a Raster Layer.
  ERROR 000840: The value is not a Mosaic Layer.
  Failed to execute (SelectLayerByAttribute).

I've searched thoroughly for an answer to this problem, and they all appear to be telling me that I need to do exactly what I'm already doing.

How am I failing here to make both of these into feature layers?

Code below...

unobscured_FC = r"C:\example.gdb\example_FC"
arcpy.MakeFeatureLayer_management(unobscured_FC, "unobscured_lyr")

#These coordinates are calculated above in a code block I've excluded...
coordList = [[[Xmax,Ymax], [Xmin,Ymax], [Xmin,Ymin], [Xmax,Ymin]]]
point = arcpy.Point()
array = arcpy.Array()
featureList = []

for feature in coordList:
    for coordPair in feature:
        point.X = coordPair[0]
        point.Y = coordPair[1]
        array.add(point)
    array.add(array.getObject(0))
    polygon = arcpy.Polygon(array)
    array.removeAll()
    featureList.append(polygon)

tempShape = os.path.join(sys.path[0], "temp.shp")
arcpy.CopyFeatures_management(featureList, tempShape)
arcpy.MakeFeatureLayer_management(tempShape, "tempShape_lyr")

arcpy.SelectLayerByLocation_management("unobscured_lyr", "INTERSECT",   "tempShape_lyr", "", "NEW_SELECTION")


texts: Cleaning GPS tracks from large jump / outliers

texts:
Can someone please explain the steps to take in order to eliminate the parts of these linestrings which jump across large distances (not along with the road network)?

They were created from GPS coordinates with this PostGIS code:

create view tour_bus_position as
select tour.id,
       st_makeline(bus_position."position" order by bus_position.record_time) as tour_bus_position_shape
from tour
join bus_position on tour.id = bus_position.tour_id
where tour.id = bus_position.tour_id
and bus_position.accuracy_meters < 100::double precision
and tour.state = 6
group by tour.id;

With this result:

Update
The next thing I tried is setting a maximum distance between the points ordered by record time:

with point_order as(
   select tour.id,
   bus_position."position",
   bus_position.record_time,
   row_number() over(partition by tour.id order by
   max(bus_position.record_time)) as point_order
from tour
join bus_position on tour.id = bus_position.tour_id
where  tour.id = bus_position.tour_id
and bus_position.accuracy_meters <= 10::double precision
and tour.state = 6
and tour.arrival_at_first_stop::date > '2018-04-25'
group by tour.id, bus_position."position", bus_position.record_time
order by bus_position.record_time asc
)
select first_point.id,
   st_makeline(first_point."position" order by first_point.record_time) as
   tour_bus_position_shape
from point_order as first_point
inner join point_order as second_point on first_point.id = second_point.id
and first_point.point_order = second_point.point_order-1
where ST_Distance(first_point."position", second_point."position") < 100
group by first_point.id;

With this improved but still not great result:

My goal is to get an accurate measurement of the total distance travelled. With the above result, it takes all of the 'noisy' points and overestimates the distances travelled.

I've tried reducing the bus_position.accuracy_meters <= 10 which helped quite a lot, but it doesn't help me deal with the points where there are large gaps in between.

I imagine there is some kind of pre-processing I can do to the points so that the linestrings turn out more realistically, or I can clean up the linestrings once they are created.

What is the best way of using PostGIS or QGIS methods?

Maybe adding a maximum distance parameter between the points to be considered in the same linestring?


texts: What is required for post-processing differential correction? (And can Garmin units provide it?)

texts:
2 very closely related questions about differential correction with a Garmin unit:

1). Some data was collected using a Garmin GPSMAP 76Cx unit and the data collectors are interested in whether post-processing differential correction is possible. Does anyone know if this model of unit supports this kind of functionality?

If I'm understanding this right (and please correct me if I'm wrong), various brands of GPS units collect a whole slew of information behind the scenes as they record points and save this in a proprietary format, such as SSF for Trimble. This is then convertible to a universal RINEX format using software (sometimes open-source?), which is then used to do post-processing differential correction. Does the Garmin GPSMAP 76Cx produce this kind of file, equivalent to, say, an SSF file on a Trimble Juno SB? I can't seem to find any such information in its detailed specs.

2). If there is no such background info collected by the unit, is any kind of post-processing differential correction possible to data collected by this unit? It may even be the case that the data collectors didn't create any points on the unit but simply wrote down the time and coordinates it displayed. As I understand it, RINEX / SSF / whatever files contain a lot of other info -- the exact "GPS time" (as opposed to local or UTC time) when each point was taken, satellite configuration and pseudo-range, a Doppler shift based on whether a satellite is approaching or retreating, number of satellites detected at time of collection, etc. etc.

Presumably if you don't have a unit that's doing all this for you, or you don't even bother to log the point on the unit but just write down the coordinate & time, you're out of luck if you want to do post-processiong correction?

Would be great if you all could confirm! Thanks so much.


texts: Calculate z score into new field in ModelBuilder?

texts:
I'm trying to calculate z scores for a given column. I am working with a shapefile, that contains fields like percent poverty, percent tree canopy and temperature. My ultimate goal is to add new columns and then calculate z scores for each these fields.

I can calculate a single z score using the procedure described in
Deriving Z-score values for given field using ArcGIS for Desktop?

However, this is something that I'll be doing a lot of times for different geographies and fields. I'd like to be able to run the process in ModelBuilder. I'm not sure if I should create a new tool, or if I can write a script within the "Calculate Field" tool. My Python is not so great, so I'm struggling to make either approach work.

I recognize that this question is similar to Calculating Z-scores using ModelBuilder in ArcGIS Desktop?, but I'm not sure how to convert the responses in that post, which are for a raster dataset, to what I'm working on.

So, does anyone know of a pre-existing tool to do this, or is someone willing to offer up some code to make this work?


texts: Is it possible to use FCIS: Fully Convolutional Instance-aware Semantic Segmentation without using CUDA?

texts:
I was reading a work on github (https://github.com/chrieke/InstanceSegmentation_Sentinel2) where it requires the use of FCIS and mxnet for deep learning, following the tutorial I saw that the installation of CUDA is required an extension of NVIDIA while I have available a video card AMD Radeon HD Series.

I wanted to know if I have trouble proceeding because I don't have nvidia or can I use OpenCL as a better alternative?

I started with running the command sh ./init.sh but it returns:

skipping 'bbox.c' Cython extension (up-to-date)
running build_ext
Traceback (most recent call last):
  File "setup_linux.py", line 56, in <module>
  CUDA = locate_cuda()
File "setup_linux.py", line 51, in locate_cuda
  for k, v in cudaconfig.iteritems():
AttributeError: 'dict' object has no attribute 'iteritems'
Traceback (most recent call last):
File "setup_linux.py", line 56, in <module>
   CUDA = locate_cuda()
File "setup_linux.py", line 51, in locate_cuda
  for k, v in cudaconfig.iteritems():
AttributeError: 'dict' object has no attribute 'iteritems'

Secondly, running the command demo.py from FCIS says it doesn't find the command even though the path is correct:

./FCS/fcis/demo.py
Traceback (most recent call last):
File "./fcis/demo.py", line 21, in
from utils.image import resize, transform
File "/home/..../..../lib/utils/image.py", line 6, in
from bbox.bbox_transform import clip_boxes
File "/home/..../..../lib/bbox/bbox_transform.py", line 3, in
from bbox import bbox_overlaps_cython
ImportError: cannot import name bbox_overlaps_cython

I changed it to from bbox import bbox_overlaps_cython to from bbox.bbox import bbox_overlaps_cython,and used python 2.7 as suggested but did not solve my problem.


texts:
I am new in the field of geospatial analysis and I'm trying to learn how to process satellites images for remote sensing purposes. I think that my problem is a very basic one.

I'm following some tutorials that suggest to download Landsat CDR from Earth explorer provided by USGS. Unfortunately, despite I am correctly logged in, the option to download what I want do not appear to me (see attached images taken from one previously cited tutorial).

How can I do to resolve this problem?


texts: Postgis is Broken after installing Nominatim

texts:
I followed Nominatim Tutorial to install PostGIS and after that I can't access PostGIS at all
I already trying reinstalling and installing everything over and over again with no success.
Everytime I try the function

SELECT postgis_lib_version()

I get the error:

ERROR:  could not load library "/usr/lib/postgresql/9.3/lib/postgis-2.1.so": /usr/lib/postgresql/9.3/lib/postgis-2.1.so: undefined symbol: HeapTupleHeaderGetDatum

My Ubuntu is

  Distributor ID:   Ubuntu
  Description:    Ubuntu 12.04.4 LTS
  Release:    12.04
  Codename:   precise

My Postgresql is 9.3 with PostGIS 9.1

The Libs I was instructed to install is:

  apt-get install build-essential libxml2-dev libgeos-dev libpq-dev libbz2-dev libtool automake libproj-dev
  apt-get install gcc proj-bin libgeos-c1 git osmosis
  apt-get install php5 php-pear php5-pgsql php5-json php-db
  apt-get install postgresql postgis postgresql-contrib postgresql-9.1-postgis postgresql-server-dev-9.1

In the tutorial here


texts: Hansen loss and lossyear bands do not always match in Earth Engine

texts:
I'm looking at the Hansen forest loss v1.5 data in Earth Engine. According to the Hansen documentation,

  In addition, to reduce confusion, beginning with version 1.4 we are no longer releasing loss as a separate layer from lossyear. Loss as previously releaed corresponds to nonzero values of loss year.
  (link)

If I understand correctly, this means that the loss band is deprecated, and only the lossyear band should be used. However, the Hansen v1.5 data in Earth Engine still has a loss band (link). Furthermore, there are pixels where the loss band is equal to 1 (indicating a loss has occurred), but the lossyear band is equal to 0 (indicating that no loss has occurred) (see example script in Earth Engine).

What is the cause of this discrepancy? In the pixels where loss == 1, but lossyear == 0, are those pixels that were classified as loss in the last year that Hansen used the loss pixel, but in subsequent versions have not been classified as loss? It seems like this might be the case, as version of Hansen prior to 1.4 do not seem to have this discrepancy (see example script).

Possibly related: there appears to be a typo in the Earth Engine documentation for the Hansen dataset. According to the EE docs, forest loss is represented by the lossyear band, and is "Encoded as either 0 (no loss) or else a value in the range 1–16, representing loss detected primarily in the year 2001–2017, respectively" (link). I don't think I'm totally losing my mind: there are sixteen numbers in the range 1-16, but 17 numbers in the range 2001-2017. (Right?)

I'm inclined to assume that the loss band is legacy and only lossyear should be used, except that the Hansen forest loss Earth Engine app and an example in the Earth Engine code playground both treat a 0 lossyear value as a loss in the year 2000, instead of no loss, as described by the Hansen documentation. I'm pretty sure the EE app was produced by the UMD team, making things more unclear. (EE App, EE example script).


texts: Winkel II fails while reprojecting on the fly in QGIS 3.12

texts:
Currently, in QGIS 3.12 the Winkel Tripel projection working on the fly but Winkel II, which is also implemented, is not working.

This is the fact or there is a bug?

If I reproject dataset into ESRI:54019 (in QGIS 3.10-12 CRS are named ESRI instead of EPSG), it is working.

When reprojecting dataset to The Winkel Tripel (ESRI:54042) it is working well.

Comparisson from here Compare Projections Winkel II vs. Mercator does not really answer my question.


texts:
I am using this code

BANDS=['B1', 'B2', 'B3','B4', 'B5', 'B6', 'B7']
area = ee.FeatureCollection('ft:1GxT2Q22KIGDjLFY1J11IX8lTAQPIoNyvBsG2fX0')
collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR').filterBounds(area).filterDate('2019-10-15', '2020-02-17')\
.select(BANDS).median().clip(area)
map=folium.Map(location=([-10.892409779059864,-46.763759689976126]),zoom_start=10)
map.addLayer(collection.select(['B7','B5','B4']),{'min':0,'max':0.5,'gamma':0.7})
map

But i got this error assent is not found

Collection.loadTable: Collection asset 'ft:1GxT2Q22KIGDjLFY1J11IX8lTAQPIoNyvBsG2fX0' not found.


texts: Plotting great circle lines in Cesium (just getting started)?

texts:
How does one plot Great circle lines in Cesium (3D, not 2D) using the WGS84 ellipsoid? I found something like this in another thread but I'm having trouble applying it to my application (simply plotting great circle PolyLines of a given width and colour.

Also, I noticed that Google Maps have a light blue line within a dark blue line. Is it possible to draw this on Cesium in one go, or should I use two separate lines? And how would I place markers at the beginning and end?

As you can see, I'm just getting started with Cesium, but I would like some insight.

I am trying to plot great circle lines between airports. The data for each airport is in a JSON object, and the sub-objects look like this:

var airportData = {
    ...
    "EDDM": {
        "icao": "EDDM",
        "iata": "MUC",
        "name": "Munich International Airport",
        "city": "Munich",
        "state": "Bavaria",
        "country": "DE",
        "elevation": 1487,
        "lat": 48.3538017273,
        "lon": 11.7861003876,
        "tz": "Europe\/Berlin"
    },
    ...
    "KSFO": {
        "icao": "KSFO",
        "iata": "SFO",
        "name": "San Francisco International Airport",
        "city": "San Francisco",
        "state": "California",
        "country": "US",
        "elevation": 13,
        "lat": 37.6189994812,
        "lon": -122.375,
        "tz": "America\/Los_Angeles"
    }
    ...
}

And the function to actually plot the lines:

function plotLine(start, end)
{
    var viewer = new Cesium.Viewer("map");

    var origin_lat = start.lat;
    var origin_lon = start.lon;
    var dest_lat = end.lat;
    var dest_lon = end.lon;

    // What goes here?
}

And if someone could put an example that works in the Cesium Sandcastle that would be great.


texts: Getting IndexError while classifying .tif using GDAL python

texts:
In the following code (derived from code written by dmh126) I want to classify .tif file based on it's pixel values after converting it into an array. I think there is some mistake in for loop which is giving me an error "IndexError: index 101 is out of bounds for axis 0 with size 101". I am posting my code here.

from osgeo import gdal
import numpy as np
driver = gdal.GetDriverByName( 'GTiff')
file = gdal.Open( 'Gadchiroli.tif')
band = file.GetRasterBand(1)
lista = band.ReadAsArray(0,0,file.RasterXSize,file.RasterYSize).astype(np.int)
print(lista)
print(range(file.RasterXSize))
print(range(file.RasterYSize))

# reclassification
for j in range(file.RasterXSize-1):
    s=lista[j]
    print(s)
    i=0
    #print(j,end='')
    for i in  range(file.RasterYSize-1):
        print(i,j,end='')
        if s[i] < 280:
           print(s[i])
           lista[i,j] = 1
        elif 280 < s[i]< 290:
            lista[i,j] = 2
        elif 290 < s[i] < 300:
            lista[i,j] = 3
        elif 300 < s[i] < 310:
            lista[i,j] = 4
        elif 310 < s[i] < 320:
          #print(s[i])
            lista[i,j] = 5
        elif 320 < s[i] < 330:
            lista[i,j] = 6
       # print(s[i])
        else:
            lista[i,j] = 7

# create new file
file2 = driver.Create( 'class.tif', file.RasterXSize , file.RasterYSize , 1)
file2.GetRasterBand(1).WriteArray(lista)

# spatial ref system
proj = file.GetProjection()
georef = file.GetGeoTransform()
file2.SetProjection(proj)
file2.SetGeoTransform(georef)
file2.FlushCache()


texts: Spatial distribution based on point data

texts: I have a specified geometric boundry with three data sets that correspond to rainfall measurements at three stations of which two are outside of the geometry and one inside. I would like to find a rainfall distribution (no matter how approximate it is) for the whole geometry according to these three data sets. How could I do this?


texts: Saving Starting / Ending points of connecting line in attribute table

texts:
I have 3 layers ( 1 line and 2 points) as shown in image. The layer id,s are mixture of integers and text such as cable4389175, node1 , service_point5316. I connect the two point layers (node and service_point) with line by using the IDs. I got the connections but the attribute table is not showing the the data of starting(fromNode) and ending points (toNode) of connecting line. For instance, node1 is connected with service_point53168 by forming a line (cable4389175).
I have also attached a demo plot to demonstrate what I am trying to achieve.

Is there a plugin or solution that will connect the points with line and also save the starting and ending points in the attribute table?


texts:
i download data from here https://nsidc.org/data/SPL2SMP/versions/3#,  but i can not open it in QGIS.

Do you know how to open this data in QGIS?


texts:
Using: QGIS 2.18 with a plugin: ECOPLAN Scenario evaluator, inserting self-made shapefiles made in Arcmap.

Received this error:

Traceback (most recent call last):
              File "C:/Users/PAU-EXT-LOVE/.qgis2/python/plugins\ECOPLAN_SE\ecoplan_se.py", line 7123, in runScenario
                CreateScenarioMask(ShpFile, RasterFile, GeoDbLandcover)
              File "C:/Users/PAU-EXT-LOVE/.qgis2/python/plugins\ECOPLAN_SE\ecoplan_se.py", line 266, in CreateScenarioMask
                SourceLayer = SourceShapefile.GetLayer()
            AttributeError: 'NoneType' object has no attribute 'GetLayer'

I have no knowledge of coding. Meaning I am not helped with previous answers.
What can I do to prevent this error? What does this mean?


texts: Labeling only a part of the attribute value of a string column in GeoServer

texts:
I am looking for a solution for labeling points based on an string attribute value in my table in a way that it shows only a portion of the attribute value in GeoServer. I am using ysld (sld) styling. The idea is that wherever there is an 'AND' in the string, the label should show only the portion before the 'AND'. For example in the table below there are few rows with values that have AND in them.

When the style is applied to this table it should show only 'ABC' for row 1, 'GHIJ' for row 2, and so on.
Below is the style I have without the split in the attribute
  - filter: ${name <> 'OWNER INFO NA'}
    scale: [min, 1067]
    symbolizers:
    - text:
        label: ${name}
        fill-color: '#FFFFFF'
        font-family: DejaVu Sans
        font-size: 10
        font-weight: bold
        anchor: [0.5,0.5]
        rotation: 0
        placement: point
        x-autoWrap: 50
        x-spaceAround: 10
        halo:
          radius: 1.25
          fill-color: '#000000'
          placement: point

What should I add in the 'label:' function in the style to achieve that?


texts:
I looking a way to integrate Esri solutions with SAP solutions.

I have seen APOS http://www.apos.com/ESRIGIS/ but don't have any experience with integration betweend Esri and SAP. So then maybe someone who had similar problem can give any hint.


texts: Enabling SQL Server 2008 database to access it as Esri geodatabase via ArcSDE using ArcGIS 10.0?

texts:
I have 1 SQL Server 2008 and 1 Oracle database (did not check which version it is).

I'd like to convert them to geodatabases, first that 'MS SQL Server 2008', then Oracle. So, leave Oracle for a while. My question is about SQL Server.

Database has about 200 tables, each under dbo schema. I used sdeupdate, got:
SDE Error(-51) Underlying DBMS error

-51 was since 2002 so I did not even try to dig more.

Let it be underlying.

So, is it possible at all with ArcGIS 10.0 to take a simple database and upgrade to geodatabase?

c:\somewhere> sdesetup -o upgrade -d sqlserver -s . -u sa -p SA_PASSWORD -D DATABASE_NAME -N -l "arcsdeserver,100,************,31-dec-2012,******************"

ESRI ArcSDE Server Setup Utility Tue Feb 28 17:51:49 2012
----------------------------------------------------------------
Upgrading ArcSDE schema.....
Successfully upgraded ArcSDE schema.

Upgrading geodatabase schema:
Error: Underlying DBMS error (-51).
       Error upgrading GDB schema.
Geodatabase schema object install not completed.

Seems, it's gonna take some while. Will prepare batch scripts to restore DBs, install profiler (seems it doesn't come bundled with MSSQLSRV2008)...

[Update]
Here I tried to play with empty databases.

On a clean, just created SQL database, using sde schema :

sdesetup -o install -d sqlserver -s . -u sde ...

OK

On a clean, just created SQL database, using dbo schema :

sdesetup -o install -d sqlserver -s . -u dbo ...

All fine again. Esri systables are named now like dbo.GDB_ITEMRELATIONSHIPS, etc.

On a clean, just created SQL database + new empty table dbo.Table_1, using dbo schema :

sdesetup -o install -d sqlserver -s . -u dbo ...

Fine again.

Now I Connect from ArcCatalog Add Spatial DB connection with:

service string= sde:sqlserver:.
user=sa
version: dbo.DEFAULT

OK - great success!

So it's possible to upgrade existing SQL database to geodatabase.

(Seems someone had fk-uped customer DB version, this I'll need to check)...

EDIT, This part is about -51 error.

Found there were sde and arcims users in DB. This probably was left from older esri software installed long ago. Removed both.
In sde_setup.log I see this:

In sde_setup.log

Successfully created  GDB_Tables_Last_Modified table...
Error creating  GDB_Items table...
ERROR Creating Geodatabase tables
, Error = -51
,EXT_Error = 1934
,EXT_ERROR1 = Microsoft SQL Server Native Client 10.0: CREATE XML INDEX failed because the following SET options have incorrect settings: 'ARITHABORT'. Verify that SET options are correct for use with indexed views and/or indexes on computed columns and/or filtered indexes and/or query notifications and/or XML data type methods and/or spatial index operations.
,EXT_ERROR2 =
,

It seems for me it may be that I missed SET BLOBSIZE or something options...
UPDATE:

Repeated all steps and additionally ran this SQL script:

-- SQLServer2008R2_SetRequiredMagicForXML.sql
SET ANSI_NULLS ON
SET ANSI_PADDING ON
SET ANSI_WARNINGS ON
SET ARITHABORT ON
SET CONCAT_NULL_YIELDS_NULL ON
SET NUMERIC_ROUNDABORT OFF
SET QUOTED_IDENTIFIER ON

Still that ARITHABORT error. Ok there are only 2 options ON and OFF, so will invert that setting and try again. Just too long and useless error message.

Update:
SET ARITHABORT OFF did not work. Played with empty new database, once set it ON and then OFF. And in both cases result was OK. This means ARITHABORT doesn't matter. That's incorrect error message. Too many if-s in code.

Update:
Replaced our code with strategy OLEDb using best practices in the world.

IMPORTANT: I'm playing with backup of customer's database on my machine. Never try to do anything like this quickly on production server without trying it in sandbox first. Use DB backups.

UPDATE: finally proved myself with your help that it's possible when playing with small databases. So final answer is YES.

Always - NO, sometimes you must fail fast.

Later I will check FAQ how to correctly "close a case" here.