texts: Spatially Analyzing two layers within same File Geodatabase using Python / Shapely?

texts:
I'm using a nested for-loop to do Shapely analysis between two Python lists created via two pre-reads of the same GDB file.  One list contains all the features of a particular layer while the other list contains all the features of another layer.

I'm using an ~ 80MB test GDB to try and figure things out in an attempt to do an all-Python approach.

My code takes about 44 minutes to fully analyze all permutations on a thin-client running Windows 7 connected to a Linux VM via Remote Desktop where my code runs.  It's what I'm limited to.

By comparison, a co-worker of mine is using Arcpy's MakeFeatureLayer_management and SelectLayerByLocation_management calls to go through the same GDB in 24 seconds directly on a Windows box and 8 minutes on a thin-client.  Basically the thin-client is 20x slower.

Roughly speaking, I'm speculating that my code SHOULD have a run time of about 2.2 minutes running directly on a Linux box based on my co-worker's observation of his own code.  2.2 minutes is still way too long.

What approach should I take to reduce my processing time?  I'm very new to GIS, so I'm not sure if there is a more obvious faster way of achieving what I'm trying to do.

NEW:

with fiona.drivers():
    try:
        if os.path.isdir(sourceDatabaseFile):
            numBLDGwithinPerennialWaters = 0
            shapelyEvaulationsAvoided = 0
            discoveredBuildings = [] # List of discovered buildings to avoid re-evaluating in Shapely
            listOfPerennialWaters = []
            listOfAllBuildings = []

           with fiona.open(sourceDatabaseFile, 'r', layer="HydrographySrf") as gdb: # Pointer to HydrographySrf layer of source file

                for hydroFeature in gdb:
                    if (len(gdb) > 0):
                        if (('ZI024_HYP' in hydroFeature['properties']) and (hydroFeature['properties']['ZI024_HYP'] > 0)):
                            listOfPerennialWaters.append(hydroFeature)
                    else:
                        print "No bodies of water to parse in hydrographySrf!!!\n"
                        break

                print "Finished pre-parsing HydrographySrf layer.  ", str(len(listOfPerennialWaters)) + " perennial waters discovered!!!\n"

            with fiona.open(sourceDatabaseFile, 'r', layer="StructurePnt") as gdbStructures: # Pointer to StructurePnt layer of source file
                for building in gdbStructures:
                    if (len(gdbStructures) > 0):
                        listOfAllBuildings.append(building)
                    else:
                        print "No buildings to parse in structurePnt!!!\n"
                        break
                print "Finished pre-parsing StructurePnt layer.  ", str(len(listOfAllBuildings)) + " buildings discovered!!!\n"

            outerLoop = 0

            for perennialWater in listOfPerennialWaters:
                outerLoop = outerLoop + 1
                for building in listOfAllBuildings:
                    if (building['id'] in discoveredBuildings):
                        print "Building already discovered so skipping to next building for evaluation! \n"
                        shapelyEvaulationsAvoided = shapelyEvaulationsAvoided + 1
                        continue
                    elif shape(perennialWater['geometry']).contains(shape(building['geometry'])):
                        discoveredBuildings.append(building['id'])
                        print "Found building ID: ", str(building['id']), "contained in perennial waters!!!\n"
                        numBLDGwithinPerennialWaters = numBLDGwithinPerennialWaters + 1


texts: Definitions and limitations of MapInfo TAB file objects

texts:
I am working on a MapInfo plugin to convert MapInfo files to GeoJSON files.

I am having trouble matching up the object types of the two formats. Each GeoJSON object type is described in the documentation. The MapInfo object types are a bit more ambiguous. The best I can find is in the MapBasic reference file which says, for example:

OBJ_INFO_TYPE value: OBJ_TYPE_REGION
ID: 7
Object type: Region Object

This is obvious enough but what exactly is a region? Is it analogous to a Polygon? Or is it any closed object?

Can anyone point me towards a decent definition of each object type?

Update 2010-02-03:

I sent a message to the MapInfo support people who pointed me towards a section in the MapInfo User Guide headed "Map Objects as Part of Layers" which spells things out perfectly. I still can't find this page in the guide that I have (version 10.5) so I'll paste the relevant section below for anyone else in my position.

Quote starts here:

MapInfo Professional are made up of layers of map objects.
There are five basic types of objects:

Regions: closed objects that cover a given area. These include polygons, ellipses, and
rectangles. For example, country boundaries, postal code boundaries, sales territories.
Point objects: represent single
locations of data. For example,
customer locations, restaurants,
parking meters. Points can also be
combined into multipoint objects.
Line objects: open objects that cover
a given distance. These include
lines, polylines, and arcs. Examples
are streets, rivers, power lines.
Text objects: text that describes a map or another object, such as labels and titles.
Collection objects: combination of region, line, and multipoint objects.

You can have each type of object in a separate layer (most common), or you can combine
objects in the same layer.


texts: How to add label to point without marker displaying in QGIS with python

texts:
I'm using a python script to manipulate the display of shapefiles in QGis.

My code basically selects an active layer and sets a field as the label of the point. Like this:

layer = qgis.utils.iface.activeLayer()
label = QgsPalLayerSettings()
label.readFromLayer(layer)
label.enabled = True
label.fieldName = 'objnam'
label.placement= QgsPalLayerSettings.OverPoint
label.setDataDefinedProperty(QgsPalLayerSettings.Size,True,True,'8','')
label.writeToLayer(layer)
QgsMapLayerRegistry.instance().addMapLayers([layer])

This give's me the following result:

I would like to add code so that the coloured marker will not be visible.
I've tried adding the following above my code in attempt to make the marker transparent:

renderer = layer.rendererV2()
symbol = QgsFillSymbolV2.createSimple({'style':'solid', 'color': '255,255,255,0', 'style_border':'no'})
renderer.setSymbol(symbol)

But this seems to do nothing.
I have a feeling that I misunderstand how the symbol layer works.
Any help or guidance is welcomed.


texts: "dataset does not exist or is not supported" using large GRID raster

texts:
I am trying my hand at arcpy for the first time to create a conditional raster from a set of threshold values. I have 3 large GRID format rasters (mos_zb3, mos_zb4, pca_finb2) I am trying to feed into arcpy using the following code:

import arcpy

from arcpy import env

from arcpy.sa import *

arcpy.CheckOutExtension = "Spatial"

env.workspace = "C:\final_mosaic_files\for_python"

Z3 = Raster("mos_zb3")

Z4 = Raster("mos_zb4")

PC2 = Raster("pca_finb2")

outCon = Con(Z4 < -0.9274,2, Con(Z3 <-0.6223,2, Con(PC2 > 1.13,2, Con(Z3 <-0.05342,2,1))))

outCon.save("C:\classifications\actual_real_deal\CART")

I keep getting "dataset does not exist or is not supported" errors in the beginning when python tries to read the rasters in as objects. I know the rasters are in the workspace folder, and I have used them in several other analyses so am confident about their integrity. I have searched the other questions about this error message on the exchange and don't see a direct answer for my problem.


texts: Test case for geo-distance implementation?

texts:
I wrote a method, based on this presentation, that filters points based on their proximity to a target point. I would like to write a few test case. What would be good test cases for this method?

// The method returns all points in the database that are
// close enough to the target point
List<Point> findPointsNear(double longitude,
                           double latitude,
                           double distanceInKilometers);

Edit - I'm more interested in a sanity test or two, to make sure I didn't screw anything while adapting/implementing the algorithm. Do you know of existing test code I could use? I'm ok with the method failing near the actual poles, as long it works in most countries.


texts:
I am using the NEDSnapshot of the NHDPlus Dataset as my raster files. They can be found at  https://www.epa.gov/waterdata/nhdplus-california-data-vector-processing-unit-18. I have included an imgur album of my various files if that is helpful http://imgur.com/a/oVT3g. The album contains images of the individual rasters, the vrt, and the mosaic file which is faulty. Any thoughts on how to resolve this issue? I'm new to both python and gdal, so I'm kind of lost as to what I can do to resolve this issue.

My code:

import os
f1 = r'C:\Users\User\Data\NHDPlusCA\NHDPlus18\NEDSnapshot\Ned18a\elev_cm'
f2 = r'C:\Users\User\Data\NHDPlusCA\NHDPlus18\NEDSnapshot\Ned18b\elev_cm'

os.system('gdalbuildvrt mosaic_CA.vrt {} {}'.format(f1, f2))
os.system('gdalwarp mosaic_CA.vrt -co BIGTIFF=YES --config GDAL_CACHEMAX 6000 -wm 500 mosaic_CA.tif')


texts: QGIS spatial join, calculating distance from polygon to point layer

texts:
I am trying to perform a kind of 'spatial join' in Arcgis, joining to a polygon layer atributes of the closest feature of a point layer, as well as distance. But 'spatial query' doesn't offer an option 'closest feature on target layer'.

I guess I could create a distance matrix, I've tried but the distances are something like 4695923.10207, can it be because I am working with WGS 84 non projected layers?

And, is it possible to perform the spatial join I mentioned, resulting on a layer, instead of creating a distance matrix?


texts: Problems with Zonal Statistic in graphical modeler QGIS 3

texts:
I want to make an average of two rasters previously generated in a model but it always gives error. I have used both SAGA and GDAL to calculate the raster out of graphical modeler without problems, but using the algorithm in the model gives the following error:

File "C: /PROGRA~1/QGIS3~1.8/apps/qgis /./ python / plugins \ processing \ algs \ saga \ SagaAlgorithm.py", line 187, in processAlgorithm
if parameters [param.name ()]. source (). lower (). endswith ('sdat'):
AttributeError: 'QgsProcessingOutputLayerDefinition' object has no attribute 'source'


texts: How to control field properties when importing csv to Arc gdb

texts:
I am writing Python code to import csv files with UK postcodes to ArcGIS file gdb.
The import works fine, however, text fields in newly created table in gdb have length of 8000 characters.

Is there a way to control field properties in a table before importing from csv?


texts: Measuring distance to two Points in ArcPy?

texts:
I would like to know how to measure the distance from two Points in ArcPy?

I know only Geometry has the distanceTo method but how do I could make a typecast in ArcPy?

I got the following error message:

Runtime error
Traceback (most recent call last):
  File "<string>", line 10, in <module>
AttributeError: 'Point' object has no attribute 'distanceTo'

Here is my code:

'''

... ptIni = {}
... ptEnd = {}
... for row in arcpy.SearchCursor("polilinha"):
...     ptIni = row.Shape.firstPoint
...     ptEnd = row.Shape.lastPoint
... print ptIni.distanceTo(ptFim)


texts: Conducting geographically weighted regression on a raster

texts:
I have a relatively big raster file on which I want to conduct geographically weighted regression. Size of raster ~ 2 GB and number of rows is around 100,000

One way might be to convert from raster to polygon and then use the GWR in ArcGIS. However, the raster to polygon conversion is taking too long.

Is there any function in R or python which can conduct GWR on big rasters? Alternatively, any approach on how to sample a huge raster so that I can run GWR on the smaller sample?


texts: Locking annotations to base map in Pisa Composer?

texts:
In Pisa Composer: I am creating a complex geologic map.

I created geologic base map and then added annotations (notes, arrows) in Composer.

I need to lock the annotations to the base map in Composer to be able to create different figures by changing scales & panning.

Group and lock?

I cannot find the documentation for this.


texts: Viewing large set of polygons in PostGIS using OpenLayers?

texts:
I am prety new to GIS. I have set up a postgis database with al large dataset and a openlayers front-end to display it. But since the dataset is getting prety big I am looking for a way to only load the polygons that are currently visible. Currently I am looking at PHP mapscript. Is it posible to achieve wat I want with this or do i need mapserver or something else for this?

I hope someone can point me in the right direction.


texts: Calculating average raster value for area in ArcGIS Desktop?

texts:
I have Raster data for a country which has Evaporation Data in its pixels. I also have Vector layer which has the 30 cities of the country.

What I need is the average Evaporation value of each city from the raster data using ArcGIS?

I have access to the spatial analyst extension.


texts: Tool to retrieve coordinates from SVG file

texts:
I have an .SVG file of a map with known extreme North/South/East/West coordinates of the area, and I'm required to retrieve coordinates of certain points on it.

What free/freeware software tools, capable of running on Windows 10 x64 (LTSB) I could use for that purpose?


texts:
I have a CSV file with climate data, latitude, and longitude coordinates.
RH      PREC    TEMP      WS      lat     lon
65.93   1.51    287.43    4.67    63.31   177.62
70.43   1.89    285.32    1.54    63.55   177.98

When I load it into QGIS I get a point shapefile (as expected).

Right now, everthing is more or less 1D. I know the shape should be something like (76,109), see image above. And then I would add my different data variables to these in the NetCDF file with a resolution of 0.3.
I tried this but got the same scattered plot but now in a NetCDF file. I actually want the points to be rasterised at a 0.3 degree resolution with all the climate data included and then made into a NetCDF file, like this.

So, how do I convert the points (first image) into a NetCDF file with gridcells of 0.3 degree? I am flexible with program to use, as long as it is not ArcGIS and Java.
1:


texts:
I'm new to FME. I have a process that intersects a line buffer from Oracle with a shp stored in local. In the end there is a dissolver to reduce the number of polygons and write it in a shp.

The problem is that from Oracle come 11 million records ... and in the dissolver I get an error.

My question would be if there would be a way to execute the process in the following way: that I read 500,000 records, process them, write them in the shp, then another 500,000, then others and so on until the end. It could be done?

This is an image of my workspace. As you can see, lines buffers in workspace, not in oracle. The attribute to Group_by in the dissolver is inherited from oracle, and all the lines have same value on this attribute. I want the buffer as clean as possible, without overlaps, but I don't know if this is the best way to do that. If I only have a a single value for the exposed attribute, can I execute Parallel Processing? Or should I create artificial grouping to be able to run the dissolver through parallel processing?

To get the intersection I use Clipper and works fine.


texts: OpenLayers 3.11.1 - Replace default attribution by an image

texts:
I am looking for a way to replace the default attribution of my map by an image.
I am trying to add a static legend to my map this way.

Any idea how this could be done?


texts: I'm using QGIS 2.18.16 and notice that every time I open a project I sit and wait while it pre-loads composers.  Is there any way to turn off the pre-loading of composers?  There's quite often when all I need to do is open the model to check something quickly, it takes longer to load the composers than it does for me to find the information I needed.  It would be nice to simply load composers as-needed, not by bulk pre-load.


texts:
Possible Duplicate:
  Can QGIS open KMZ files?

Is it possible to import Google Earth generated KML or KMZ files into QGIS and then use these files as background data for a project?

I intend to add non-KML hydrology data to these base layers.

Can you provide some easy tutorials?


texts: How to compute the distance between centroids and polygons?

texts:
Suppose I have an Excel spreadsheet with the latitude and longitude of the centroid of every county in the U.S.  I want to know the shortest distance from each county to any part of Oklahoma (presumably, somewhere along the border).  I then want this distance in this Excel spreadsheet as a new column.

Is this possible to do and how?  I know next to nothing about GIS


texts: Dissolve/merge polygons using ST_Union in CARTO

texts:
I am attempting to dissolve the boundaries between several polygons using a common attribute within a column in CARTO and the ST_Union function.

In my scenario, I have created a table for all US Zip Codes, & have given each Zip Code a value within a column called area_zip_r. I would like to dissolve/merge the polygon boundaries for all zip codes that share a common "area_zip_r" value, creating a single polygon for each value within "area_zip_r".

After some initial research and the help of this Stackexchange article I have come up with the following command which does dissolve/merge the polygons, but also leaves gaps of data for an unknown reason.

SELECT
  cartodb_id,
  area_zip_r,
  ST_union(the_geom_webmercator) as the_geom_webmercator
FROM
  {table_name}
GROUP BY
  area_zip_r,
  cartodb_id

Results: Many of the area_zip_r values are seen in the output with the desired effects, but I only see a small percent of the area_zip_r values I expect to see. I am returned with approximately 500 unique area_zip_r values when I expect 2000. This is verified when mapping the results. I've attached photographs to show the area my original table covers, and the results of my ST_Union command.

What am I doing wrong?


texts: Histogram matching in Google Earth Engine?

texts:
Does anyone know if there is a more GEE-friendly way of achieving this? I'm experimenting with S1 SAR image histogram matching in GEE and was wondering if foremost the code can be improved. Read here that the use of for-loops and getInfo() calls should be discouraged due to the overhead it brings with transferring data downstream to the browser and upstream to the server.

Here's the code:

/**
* Function that returns a normalized image using unitScale()
* @param image the image object
* @param band the band to select
* @param aoi the area of interest for scoping the reducer
*/
function getNormalizedImage(image, band, aoi) {
  var redux = getReducedImage(image, ee.Reducer.minMax(), aoi);
  return image.unitScale(
    ee.Number(redux.get(band + '_min')),
    ee.Number(redux.get(band + '_max'))
  );
}

/**
 * Function wrapper for a reducer
 * @param image the image object
 * @param reducer the reducer to apply, see API docs
 * @param aoi the area of interest for scoping the reducer
 */
function getReducedImage(image, reducer, aoi) {
  return image.reduceRegion({
    reducer: reducer,
    geometry: aoi,
    scale: 10,
    bestEffort: true
  });
}

/**
 * Function that returns a histogram array object
 * @param image the image object
 * @param band the band to select
 * @param aoi the area of interest for scoping the reducer
 */
function getHistogram(image, band, aoi) {
  var redux = image.reduceRegion({
    reducer: ee.Reducer.histogram(256, 1),
    geometry: aoi,
    scale: 10,
    bestEffort: true
  });
  var histoDict = ee.Dictionary(redux.get(band)).get('histogram');
  return histoDict;
}

/**
 * Function that returns a cumulative distribution function
 * as a list object
 * @param image the image object
 * @param band the band to select
 * @param aoi the area of interest for scoping the reducer
 */
function getCDF(image, band, aoi) {
  var histoArr = ee.List(getHistogram(image.multiply(255), band, aoi));
  var len = histoArr.length().subtract(1);
  var nPixels = ee.Array(histoArr).accum(0).toList().get(len);
  var CDF = ee.Array(histoArr).accum(0).divide(nPixels).multiply(len);
  return [CDF.round().toList(), len];
}

/**
 * Function that returns a source image that has been
 * matched with the target's image histogram using
 * the remap function from Earth Engine
 * @param imageR the source image object
 * @param imageZ the target image object
 * @param bandR the source band to select
 * @param bandZ the target band to select
 * @param aoi the area of interest for scoping the reducer
 */
function matchHistogram(imageR, imageZ, bandR, bandZ, aoi) {
  var cdfR = getCDF(imageR, bandR, aoi);
  var cdfZ = getCDF(imageZ, bandZ, aoi);
  var lenR = cdfR[1];
  var lenZ = cdfZ[1];
  var map = [];

  var cdfR0 = cdfR[0].getInfo();
  var cdfZ0 = cdfZ[0].getInfo();

  for (var r = 0; r < cdfR0.length; r++) {
    var max = Infinity;
    var min = 0;
    for (var z = 0; z < cdfZ0.length; z++) {
      var diff = Math.pow(cdfZ0[z] - cdfR0[r], 2);
      if (diff < max) {
        max = diff;
        min = z;
      }
    }
    map.push(min);
  }

  var seq = ee.List.sequence(0, lenR);
  var imageM = imageR.multiply(lenR).round().remap(seq, map, 0).divide(lenR);
  return imageM.select(['remapped'], [bandR]);
}

getInfo() is being called in the functions matchHistogram and getCDF.

UPDATE: I've managed to improve getCDF by removing the getInfo() call there, matchHistogram seems a bit more tricky to solve.


texts:
I'm having a debate with a colleague of mine over whether generated heat maps are quantitative or qualitative in nature.

The underlying data source is GPS points. The GPS devices that we use are rigged to save power based on speed so they do not produce a uniform sampling rate.

In my mind generating heat maps for these data over say 2 weeks is a quantitative representation of space use. My argument being GPS points cluster around the areas of space most used. So to be able to say over 2 weeks this area see’s 50,000 points vs another's 20,000 points I view as a quantitative representation of area use.

My colleague's view is that the heat maps are qualitative because the sampling rate isn't steady. Some GPS devices are napping for 4 hrs, others less so, so have no idea how many points are missing from your counts for a true quantitative display.

Thoughts?


texts: Qgis 3.10 Raster Calculator error in standalone script

texts:
I want to do simple raster calculations on approx. 1'000 small rasters in a standalone script in QGIS. I verified the code several times with various SE answers and youtube tutorials. The same calculation in the raster calculator GUI works properly. In the QGIS python console however, it always throws the result "1" (CreateOutputError) to me and does not create any output raster.
I tried to change input file and output directory in order to check if there was an issue with my tif-file or with the drive I wanted to save the output to.
I don't know what else I could test? In some forums, there was an issue with raster calculator on QGIS 3.10 reported. Any hint?
My code:
    raster_layer = QgsRasterLayer('C:\Users\baer\Desktop\GIS_temp\area_0.tif', 'ras')
    output = r"C:\Users\baer\Desktop\GIS_temp\area_0_red.tif"

    entries = []
    ras = QgsRasterCalculatorEntry()
    ras.ref = raster_layer.name() + '@1'
    ras.raster = raster_layer
    ras.bandNumber = 1
    entries.append(ras)

    calc = QgsRasterCalculator('"ras@1" + 100', output, 'GTiff', raster_layer.extent(),raster_layer.width(), raster_layer.height(), entries)

    calc.processCalculation()

Result was always:
    >>> calc.processCalculation()
    1

I also tried to add
    context = QgsCoordinateTransformContext()

which didn't help. Also varying the calculation expression didn't change anything.
Enabling OpenCL acceleration, resulted in the output to be "6" (BandError) instead of "1" (CreateOutputError), but still no output raster.


texts: The stability of randomForest models after increasing predictor variables

texts:
With reference to the posts: Incorporating terrain data to predict canopy cover using randomForest in R and Random-Forest Classification of 10cm Imagery for species-distribution in R (no point-shapes)

I would like to know how and why would one's model to map canopy vs. non-canopy areas, improve if several predictor variables are added such as:

Vegetation indices
slope
elevation
aspect
multiple bands

If I do add several predictor variables to the training set, I have to predict the model over a raster stack consisting of the same layers right?

Therefore, can I incorporate bands from Landsat too (after down-scaling to 5.8 m) with my LISS IV bands, as part of the training data and the raster stack on which the model will predict?


texts: How to combine ST_Intersection() and St_Area() into a single query?

texts:
I have a queried layer and a buildings layer and I am currently running the following query to find areas of intersection using the St_Intersection() function inside the QGIS DB Manager plugin:

SELECT ST_Intersection(buildings.geom, queried.geom) as inter from buildings, queried

(I'm hoping the function is correct as it's been running for over an hour now...)

But I want to repeat this and combine it using ST_Area() so that the output layer constains all intersection features with another column in the table showing the area.

How could I combine these two functions into one query?


texts: Creating Python code to pin and rotate fixed length segment using shapely

texts:
I am very new to Python GIS.
I am creating code in Python using shapely geometry library. One of the goals is to create a fixed length segment, pin one end of it to a known coordinate and rotate the segment, keeping, of course, the length constant. I have created the line with
line = LineString([Point(px0, py0), Point(px1, py1)])

where px0 and py0 is where I want the segment to be anchored. I can now look for px1, and px2 using trigonometry. However, it would be very nice if I could simply describe this segment by its length and angle. I wonder if there is a library already built that does the trick. I have searched the shapely documentation without success.


texts: About zonal statistic function in ARCGIS

texts:
I calculated zonal statistic in Arcgis to get max/min value of each LC types. Broadly, LC image include forest, grassland, waterbody etc. But output max value of zonal statistic function give me 0.972 on the waterbody, 0.9908 on grassland, 0.9976 on forest. I think it is not reliable. What is the reason on this? Actually NDVI image has been masked on waterbody. Original value of the NDVI image range from -0.03 to 0.9976. Do you guys have another suggestion for this?

Thanks a lot


texts:
I've added the GDAL extensions to GeoServer and it worked. However, I followed the instructions in order to add GDAL ECW & MrSID Extensions for GeoServer documented in the installation guide (Adding support for ECW and MrSID on Windows) but I got the error in the screenshot below:

Considering that I've created a GDAL_DRIVER_PATH variable in the environments variables as shown below

And not appearing in the GeoServer. I am using GeoServer 2.8.1

What might be the issue here?


texts: Working with geographic coordinates in shapely?

texts:
I have a very basic question about working with the shapely module in Python. I have some geographic coordinates received from a GPS. I would like to use shapely to calculate the great cicle distance in meters between two points.  I started with:

>>> from shapely.geometry import Point
>>> p1 = Point(43.374880, -78.119956)
>>> p2 = Point(43.374868, -78.119666)

I believe that this gives me two points in the cartesian coordinate system, which isn't going to be very useful.  These coordinates come from a GPS which presumably is using the WGS84 CRS.  I have seen some example using pyproj and shapely.ops.transform to transform between coordinate systems, but these all seem to involve data points that are already in some valid cartographic coordinate system with an EPSG identifier. I'm not sure what do with naive points like the above that don't have an associated geographic CRS.

How do associate these points with the WGS84 CRS such that p1.distance(p2) will yield the distance in meters? Am I using the wrong tools?


texts: Mapping Russia - Int'l dateline splitting polygon

texts:
I am having the issue where the International Dateline is splitting a polygon I am working with.  I have a layer of the "Federal Districts" of Russia, and the eastern most federal district is being "cut" by the international dateline.

I currently have the data frame's coordinate system set to Asia Lambert Conformal Conic and this places the two "cut" polygons next to each other (how they should be).  I really need the polygon to not be cut and show as a continuous polygon.

I've researched up and down, tried changing the projection/coordinate system, but with no luck.

I've attached a screenshot of the one multi-part polygon that is being split by the international dateline.  I would like it so that there is no cutting of the polygon, and therefore no boundary line going through the polygon.

I've even tried transforming this multi-part polygon into separate, single-part polygons, then merge them back together, but anytime any part of any polygon crosses this line, it gets cut.

Layout View:

Data View:


texts: Comparing multiple polyline shapefiles to single key file

texts: I am working on a project where we are comparing participants' spatial memories. We take them through a circuit, then later in the lab they draw the route they traveled. I want to be able to compare this drawing against the actual route of travel. There are 40 in total, and I have them digitized and georeferenced, along with the "key" polyline file. there are 3 segments in each polyline layer. I am hoping to arrive at some sort of score - perhaps a "degree of difference". I had hoped to find a Levenstein dissimilarity tool or the like, but haven't found anything. Feature Compare seems to not be able to do this.


texts:
I am trying to understand the meaning of this error. I created spatial data using the sf package in R and wrote out the data as a GeoPackage (i.e. a .gpkg). However, whenever I add the data into ArcMap, I receive the following error.

Unfortunately, my reproducible R example works and my real example still does not, thus something must be wrong with my code. And I cannot figure out what it is. The real code is too long and cannot be copied here.

The weird thing is that when I write the data in R out as a shapefile, there is no error and the data is correctly spatially referenced, etc. What does this error mean? And why would it produce as a .gpkg and not as a .shp?

A similar error message was seen here, but it is a little different.


texts:
I have two tables in psql that I am trying to join based on certain attributes - they are very large - 17 M rows and 2.7 M rows respectively

mydb=> SELECT reltuples::bigint AS estimate FROM pg_class where relname='foo';
 estimate
----------
17087196
(1 row)

mydb=> SELECT reltuples::bigint AS estimate FROM pg_class where relname='bar';
 estimate
----------
  2763829
(1 row)

On both, I have created spatial indices with

CREATE INDEX foo_gix ON foo USING GIST (the_geom);

The query I am running is a point in polygon analysis based on data collected at historical time intervals - counting points in polygons when the timestamps on both tables match.  Points are locations of mobile phone connections (bar_history and bar) and polygons are buffers around certain areas (foo_history, foo)

It looks like this:

select s.id, s.place_id, s.time, count(l.location) as total
FROM foo_history as s LEFT JOIN foo as p ON s.place_id = p.id
LEFT JOIN bar_history l ON ST_Contains(ST_Buffer(ST_Transform(ST_SetSRID(ST_Centroid(p.polygon),
4326), 32615), 100), ST_Transform(l.location, 32615)) LEFT JOIN bar ON bar.phone_id = l.id
WHERE bar.network_id = 2
group by s.id LIMIT 5

This query returns results successfully in less than a few seconds.  However, when I add join by timestamps:

select s.id, s.place_id, s.time, count(l.location) as total
FROM foo_history as s LEFT JOIN foo as p ON s.place_id = p.id
LEFT JOIN bar_history l ON ST_Contains(ST_Buffer(ST_Transform(ST_SetSRID(ST_Centroid(p.polygon),
4326), 32615), 100), ST_Transform(l.location, 32615)) LEFT JOIN bar ON bar.phone_id = l.id
WHERE bar.network_id = 2
AND
to_timestamp(floor((extract('epoch' from l.time::timestamp) / 600 )) * 600) =
to_timestamp(floor((extract('epoch' from s.time::timestamp) / 600 )) * 600)
    group by s.id LIMIT 5

the query runs endlessly, even though nearly all of the timestamps will match each other from both tables.  I am not a database expert.

How can I speed up this query or add some kind of index that might speed it up?


texts:
Possible Duplicate:
  Calculating Lat/Lng X miles from point

I would like to know how could I translate a geographical coordinate to an specific distance. So for example let's suppose I'm facing North at (-78.4437,-0.31399). I need to know the coordinate that's 100m farther and has +30 degrees offset from North.
Thanks a lot.


texts: Storing raster data with Java into Oracle Spatial

texts:
I have a Java application that processes raster data.
At the end of the process I plan to store the raster in oracle spatial.

What is the best way to do this?

I had a look at geotools but I did not find a solution there.

Is 'wrapping GDAL' the best way and, if so, then how?

Or is there a workaround or a strategy available?


texts: Planimetric Area Measurement in Geographic Projected Vector Data in QGIS3

texts:
I am trying to calculate the area of agricultural fields in QGIS. If my vector data is in Default WGS 84 - EPSG:4326 geographic projection, the difference between Planimetric and WGS 84 Area Calculation becomes so large. Here is a screenshot of the comparison of results (WGS 84 & Planimetric Area Calculation in EPSG:4326 and EPSG:32635 - WGS 84 / UTM zone 35N Projected Coordinate System) The difference is much smaller in EPSG:32635 Projected Coordinate System.
You can replicate the result with any vector file. Is it all wrong to calculate the area as planar in EPSG:4236 projection?

area_wgs = WGS 84 Area Calculation
area_pla = None/Planimetric Area Calculation
in Project settings as:

OR


texts: OpenDroneMap returns error on merging depthmaps?

texts:
I am doing a reconstruction on 424 Parrot Sequoia RGB images using OpenDroneMap.  The reconstruction runs for about 2 hours then throws the  error shown below. I have successfully done several reconstructions using ODM on my computer but I have never seen this problem.

I am running ODM on Ubuntu 16.04.

Load and prepare mesh:
PLY Loader: comment PCL generated
Reading PLY: 52843 verts... 105614 faces... done.
Warning: Zero-length normals detected: 12 face normals, 0 vertex     normals
Generating texture views:
NVM: Loading file...
NVM: Number of views: 424
NVM: Number of features: 0
        Loading 53%...libpng error: Write Error
Aborted (core dumped)
Traceback (most recent call last):
  File "/home/gerry/OpenDroneMap-v0.4/run.py", line 47, in <module>
    plasm.execute(niter=1)
  File "/home/gerry/OpenDroneMap-v0.4/scripts/mvstex.py", line 139,     in process
        '{keepUnseenFaces}'.format(**kwargs))
  File "/home/gerry/OpenDroneMap-v0.4/opendm/system.py", line 34, in run
raise Exception("Child returned {}".format(retcode))
Exception: Child returned 134


texts:
I am trying to do a simple format of a text box in ArcMap. I went into the 'About Formating Text' and copied a tag to change color and when I apply the setting, the text doesn't change color and the tags appear in the text like this:

Can anyone suggest what I am doing wrong?


texts: convert a gpx file to a feature class loses ident field

texts:
Using the ArcGIS 10.3 gpx to feature tool (GPX To Features (Conversion))  the "ident" field is lost in the translation.

DNR Garmin is an option but has limits to the number of points it can process and I'm running into memory issues.

Is there another tried and tested method to translate from .gpx to .shp  whilst keeping the ident field intact or can the ESRI tool do this and I am missing somthing?

I have hundreds I am trying to batch.

I am using ArcGIS 10.3 for Desktop but am open to Python solutions.


texts: Populating polygon with most frequent point attribute using QGIS?

texts:
What I have is two features, a polygon showing wards, and a point shapefile showing properties which has an attribute detailing the type of consumer within the household, such as country living, urban cohesion, family basics etc.

The point data has a lot of features, as you'd expect, and what I am trying to do is workout what type of consumers live within a particular ward. Obviously within a single ward there will be a variety of consumers, so I am after finding the most common consumer. So ward x may be classified as mainly country living consumers, ward y family basics etc. The end result will likely be a new field in the ward shapefile with attribute for each ward with the most frequent consumer type.

I have tried a few things, in both ArcGIS and QGIS, such as spatial joins but I am not having much luck.


texts: How to map a range of values to a single value in FME

texts:
A dataset with a set of points has an attribute called HEIGHTS which is filled with numbers like, 10.8, 0.4 and 150.6 etc... I want to transform this dataset to a shapefile. In the transformation I want to add a new column to the shapefile which will need to contain a color code.

How can I map the height values to this color code?

For example all heights between 5-10 should have the color dark blue,  all heights between 10 and 20 should have the color light blue etc......
I tried to use the valuemapper transformer, but it only seems to work for 1:1 mappings.

I am using FME 2009.  This means that I cannot use the AttributeRangeMapper and AttributerRangeFilter transformers. Can I conclude that what I want is not possible in FME 2009 and I should look at a newer version on FME?


texts: Dissolving ploygons in QGIS or Grass to new polygons of defined population sizes?

texts:
I have a shapefile with 3800 polygons representing the smallest administrative areas in Ireland. Many of these, in rural areas, have tiny populations, under 100 people. In urban areas, all would have adequate populations.

For each area I have a population figure, and various other attributes, all of which can be intelligently summed to give valid results for areas produced by combining the smaller areas.

I need to merge the areas to get say, 500 or 600 areas of reasonable population. I need to leave the larger (in population) areas alone, where possible. I would greatly prefer to do this automatically. This is necessary for model fitting purposes, primarily.

I've looked at a number of the dissolve questions here:

'Batch Dissolve Polygons based on Attributes in QGIS or FWTools' I think would require me to know which areas to merge, when this is what I need to calculate.
'Practical way of managing polygons and unions of them' and 'How do dissolve polygons from shapefile using open source tools?' have the same issue.

Any suggestions?


texts: Zooming to specific loaded layer using PyQGIS

texts:
I'm trying to zoom to a layer that was loaded with a script, just like manually would be done by right clicking the layer and pressing the "Zoom to Layer" button.

I've tried several ways and the closest I got was using canvas but it wouldn't zoom in but instead it would open a secondary screen with the layer extended.


texts: UPDATE a value field for area of the geom and convert sq m to acres using cartodb

texts:
I'm trying to use the SQL commands in CartoDB to UPDATE a value field that shows total acres of the area of the_geom. My SQL is pretty poor/entry level and has some issues:

UPDATE table_name SET shape_area=(SELECT ST_Area(the_geom::geography) area_sqm
FROM table_name)

I'm missing something... any help?


texts: FME Querying an Microsoft SQL Server Non-Spatial to perform an attribute join

texts:
I am trying to join a spatial table with some values coming from a MS SQL table.

I have it working properly with an excel but when I try to use SQLExecutor or SQLCreator I get stuck where the DatabaseJoiner only allows me to choose SQL tables and not the result of the SQLCreator.

But in the DBJoiner I can't get the reader to be the result from the SQL

Do I have to write to a csv/excel and then read use this for the db joiner? it's best to avoid the creation of files even if temporary.

UPDATE:
I need to run the query to get the fields I need to join from the MSSQL -they aren't in a table but a result of a query.

Basically the process is.
1. Get current cadastre and clean it to be in the format we require (done as part of FME Merge two files that only have some attributes the same)
2. Perform 3 sets of queries from the db that *a -match lot plan to the council's land_no and get some address fields, *b -use the land_no and get a property no from a second query, *c -use the property number to get the owner info and put this into the table.

3- the queries seem to work fine but when I use DB Joiner of Feature Joiner it's overwriting the fields that were previously populated


texts: Selecting adjacent polygon in vector layer in OpenLayers 3?

texts:
I am trying to change feature of all polygon's with in 100 radians of given lat long.

initially i used buffer but instead of highlight polygon's it creating as new layer. below is the code i have tried

var poitnExtent=features[i].getGeometry().getExtent();
        var bufferedExtent = new ol.extent.buffer(poitnExtent,100);
        //console.log(bufferedExtent);
        var bufferPolygon = new ol.geom.Polygon(
        [
            [[bufferedExtent[0], bufferedExtent[1]],
            [bufferedExtent[0], bufferedExtent[3]],
            [bufferedExtent[2], bufferedExtent[3]],
            [bufferedExtent[2], bufferedExtent[1]],
            [bufferedExtent[0], bufferedExtent[1]]]
        ]
        );
        var bufferedFeature = new ol.Feature(bufferPolygon);//console.log(vectorBuffers.getSource());

        map.addLayer(vectorBuffers);
        vectorBuffers.getSource().addFeature(bufferedFeature);


texts: Keeping TOC layers unchecked on first load in ArcGIS JavaScript API

texts:
I am new to ArcGIS and it's API.
I am using NLiu's TOC widget for displaying the layers. It is working perfectly, however, I wanted to keep all layers unchecked when it is loading on the browser (for performance reason).

I used the visible property of dynamic layer as:

var layer= new esri.layers.ArcGISDynamicMapServiceLayer(mapService,{"visible": false});

This hides the layer, but when I am expanding its childs, all of them are selected (checked), when I am clicking to uncheck one, the rest others are displaying back.

So, what I want, to keep all its child nodes un-selected, when I am checking a single checkbox, so it should display its relevant data on the map.

I tried to find a way on its javascript files but not succeed.


texts:
I have a Layer:

          pointLayer = new OpenLayers.Layer.Vector("Point Layer", {
            strategies: [strategy],
            styleMap: new OpenLayers.StyleMap({
                "default": {
                    label : "${text}",
                    pointRadius: 10,
                    fillColor: "#ffcc66",
                    fillOpacity: 0.8,
                    strokeColor: "#cc6633",
                    strokeWidth: 2,
                    strokeOpacity: 0.8
                },
                "select": {
                    fillColor: "#8aeeef",
                    strokeColor: "#32a8a9"
                }
            })
        });

I create this layer without any feature. I periodically call the server and get json objects from it. The json objects have coordinate attributes and a text attribute. After i get the objects, i put points on the layer. I set the text attribute in the point features:

pointFeature.attributes = {
    text: jsonObj.text
};

The problem is that the labels don't appear on the layer (but the points appear). I noticed that the labels appeared only when i added the point features to the layer before i added the layer to the map.

I found the redraw() method, but it doesn't work because after the server call i remove all the points from the layer and add the new points to it.

To summarize the problem: i want to display labels of point features which were added to the layer after the layer was added to the map.


texts: Creating outline around layer which has separated features using ArcGIS Desktop?

texts:
I have a layer representing a neighborhood. And the roads are separating groups of features(blocks) from each other. I want an outline around the whole neighborhood and not around each feature or each block. I want a line around the whole layer. I have tried dissolving, masking and merging. None have worked.

How could I do this using ArcGIS Desktop?


texts: Exporting metadata to HTML format via ArcGIS Desktop?

texts:
In ArcGIS 9.3, it was easy to export metadata into an HTML format. In ArcGIS 10, it looks like everything is in XML format and that there's no easy way to convert the XML to HTML.

Has anyone successfully exported layer metadata into an HTML format using ArcGIS 10?


texts: ENVI: Postclassification, remove pixels of certain size

texts:
I did a classification over multiple Landsat imagery in ENVI 5.1, which turned out not too bad (ice and snow in Cyan in the picture below). Now I am about to do a postclassification to filter out pixel values which I don't want in my classification (in my case small ice bergs and snow patches) to smoothen the results.

Is there a tool in ENVI where I could define the number of pixels which should be eliminated? What I have in mind is to get rid of all areas which are smaller than let's say for example 1 km2.

How would you do this? Is there an easy way in ENVI? I've tried the Sieve Classes tool and also Majority/Minority Analysis tools up until now but I would be interested what you normally use to do this kind of thing, when area definition is important?


texts: Accessing initial point data inside of point metrics functions with lidR

texts:
I am using the R package lidR to do some analysis on a large forestry point cloud. I am trying to write my own point_metrics() function for identifying outliers. How can I access the X, Y, and Z coordinates of the point that the metric is being run on inside of the function? For example with the following code:
is.outlier <- function(x, y, z) {
  return(x)
}
point_metrics(las, ~is.outlier(X, Y, Z), k = 5)

returns the x values of the point's k nearest neighbors. I want to find the x value of the point in addition to the k nearest neighbors in order to find the distance away each knn is from the initial point.


texts: How to get the centroid of a set of points in PostGIS?

texts:
I use PostgreSQL with the PostGIS extension.

I have a set of points in the_geom column from a table myschema.myobjects.  I want to create a select statement to get the centroid of this cluster, so from a select statement like this:

SELECT the_geom FROM myschema.myobjects

I need to find the correct syntax for a statement like:

SELECT ST_AsText(ST_Centroid( (SELECT the_geom FROM myshema.myobjects) ));


texts: Why isn't my layer loading via pyqgis despite the file name and path being correct?

texts:
I am writing up a python script on PyQGIS to load a point shapefile.

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

src = "C:/Users/Aaron/Documents/pygis/Points.shp"
layer = QgsVectorLayer(src,"Localities", "ogr")
if not layer.isValid():
    print "Layer failed to load!"

else:
    print "Layer was loaded successfully!"

QgsMapLayerRegistry.instance().addMapLayer(layer)

Everytime I run this the "Layer failed to load!" message appears.

But checking the file name and path, well it appears to be correct.

What is the error that is causing the layer to fail to load?


texts: Is there a way to show a layers feature count in the TOC?

texts:
I'm using ArcGIS 10.  In ArcMap's table of contents (TOC), is there a way to have the "Layer Name" automatically show a count of the total number of features in each layer?

I was thinking the TOC would look something like this:

Roads (27)
Streams (100)
Parcels (12)

I found this option for Unique Value renders, but:

I am not an ArcObjects guy, and
I want to work with just the Single Value renderer.

The "List By Selection" tab sort of has this capability, but only when there are selected features.


texts: Single Sign In authentication to PHP app and GeoServer

texts:
We have an application written in PHP / Javascript to manage PostGis-data that is displayed using OpenLayers and GeoServer.
Since access to the data is limited to authenticated users they actually need to login twice: one in our application and a second time in GeoServer.

The main login-page is located at

https://DOMAIN/login.php
and the GeoServer login is at

https://DOMAIN/geoserver/web/

What would be the best way to realize a single sign-in, so that the users have only to log in once at the main page and the username and password are passed to the GeoServer-login?


texts: Changing Map Tip popup direction dynamically

texts:
Is there a way to dynamically (or semi-automatic) change direction of a popup?

Default direction is 'right-down', which in some cases, where hovered feature is at the bottom of a canvas, Map Tip is cut (no. 3 on the image). So It would be better to have popup placed like no. 2.

I described only one case, but it applies to all situations, where features are close to border of a canvas.


texts: While delineating catchment in QGIS using Upslope Area, error:could not create index insufficient memory prompted and software got hang again and again. What should I do?


texts: Extracting Coordinates from Oracle Spatial Object Table using FME?

texts:
I've searched all available online sources but none appear to have this issue.

I've added an Oracle Spatial Reader, connected to desired table and now want to extract x/y but I get this error when running.

  Error - Cannot have NULL value for 'z coordinate'

I've looked at the table in Oracle developer and indeed the Z values are Null.

I'm aware that this table has multiple Spatial columns so I've checked:
Handle Multiple Spatial Columns: Yes

I've tried a SELECT statement to run before translation:

SELECT c.GEOM.SDO_POINT.X AS X, c.GEOM.SDO_POINT.Y AS Y,
FROM A.ENTITIES c
WHERE ACCESS_CODE = 'O'

I've tried to expose the attributes:
oracle_sdo_point.x
oracle_sdo_point.y
oracle_sdo_point.z

Is there a way to force this to 2D rather than looking at the Z values?


texts: How to include a single quote / apostrophe in an ArcGIS Server REST API query string

texts:
This is similar to How do I query a string in Arcobjects that has a single quote in it? but in my case the apostrophe is actually in the word I'm searching for: &where=locality = 'O'CONNELL'

This fails in an ArcGIS Server query:


texts: How to delete empty feature classes in FME at the end of executing a workbench?

texts:
I have an FME workbench which uses a python parameter to copy an existing file geodatabase from one location to another and uses the copy as the writer.  This geodatabase has about 50 empty feature classes with fields defined.
The input of the workbench is a folder of shapefiles which each get mapped to the appropriate feature classes using a SchemaMapper transformer.
What I want to do is delete the empty feature classes from the geodatabase after the workbench completes.  I have ArcGIS installed on the same machine and I have been successful at deleting the empty feature classes using arcpy in an FME Shutdown Script.
The script I'm using is:
for dirpath, dirnames, filenames in arcpy.da.Walk(env.workspace, datatype = 'FeatureClass'):
    for filename in filenames:
        fc = os.path.join(dirpath, filename)
        featureCount = int(arcpy.GetCount_management(fc).getOutput(0))
        if featureCount == 0:
            arcpy.Delete_management(fc)
            print "deleted {0}".format(fc)

However, I was wondering how I could achieve this without the use of arcpy?  What would I need to do to delete the empty feature classes using only the tools available in FME?
I'm using FME 2016.0.1.


texts: Way to "rip" bicycle directions from Google Maps Classic?

texts:
It's very unfortunate that Google Maps Classic no longer exports routes as a KML file, and Google Maps Engine Lite isn't as easy to use.

So I was wondering if it were possible to somehow "rip" a route drawn in Google Maps Classic, possibly through some local JavaScript.


texts: Calculate coverage percentage of an area covered by shapefiles using GDAL？

texts:
I have thousands of shapefiles covering around the world. I need to calculate coverage percentage. For example, a pixel located in Hawaii was covered by 100 of all shapefiles, then this pixel value of 100 will be output.

I can convert them to rasters, but I cannot stack all rasters in a file for calculation because of their large volume.

So, how can I  calculate coverage percentage by using shapefiles or separated rasters?


texts:
I'm querying the ArcGIS REST API and I'm trying to write a request that includes a wildcard. The normal request (without wildcard) looks like this:

https://serviciosgis.catastrobogota.gov.co/arcgis/rest/services/ordenamientoterritorial/catastro/MapServer/6/query?f=json&where=ESOCLOTE=%22005611022010%22&spatialRel=esriSpatialRelIntersects&geometry={%22xmin%22:-8248284.5749014811590314,%22ymin%22:523903.3499047987861559,%22xmax%22:-8248034.5749014811590314,%22ymax%22:524153.3499047987861559,%22spatialReference%22:{%22wkid%22:102100,%22latestWkid%22:3857}}&outFields=*

This one works, and returns a valid json object with one feature. However, if I include a wildcard with virtually the same parameters it doesn't return anything. Here's the request with the wildcard:

https://serviciosgis.catastrobogota.gov.co/arcgis/rest/services/ordenamientoterritorial/catastro/MapServer/6/query?f=json&where=ESOCLOTE+like+%22005611022010%22&spatialRel=esriSpatialRelIntersects&geometry={%22xmin%22:-8248284.5749014811590314,%22ymin%22:523903.3499047987861559,%22xmax%22:-8248034.5749014811590314,%22ymax%22:524153.3499047987861559,%22spatialReference%22:{%22wkid%22:102100,%22latestWkid%22:3857}}&outFields=*

Am I writing the request correctly?

The first one says where=ESOCLOTE=%22005611022010%22 (this works) and the second one says where=ESOCLOTE+like+%22005611022010%22 (this doesn't work).

Any ideas?


texts: Extracting exact land coverage raster from TIFF format raster using R?

texts:
I have land coverage raster data that downloaded from JRC data catalogue (JRC website and JRC land coverage data download link), where data file like shapefiles at a very high degree of resolution, with information of land/soil coverage, for example, information on agricultural land coverage, information about mountain, city coverage and so on. I want to extract all information on agricultural and city land coverage only for Germany and render them in my final output as csv or xlsx format.

However, I tried to read those data in R by using raster::stack function and I intend to crop the raster of only Germany. I got some idea from SO community (useful post), Using Jeffrey Evans' solution, here is runnable R code down below:

url = "https://cidportal.jrc.ec.europa.eu/ftp/jrc-opendata/LUISA/PrimaryOutput/Europe/REF-2014/JRC_LUISA_Input_Corine_land_cover_2006_r_ref_2014.zip"
download.file(url, basename(url))
gunzip(basename(url))    # got some problem to unzip the file in R; better to use decompress the file in local site
tifDat <- list.files(getwd(), "tif$")     # tif file was located in two nested folder
land_cover = raster::stack(tifDat[1])

I plotted this tif format raster in R down below:

plot(land_cover)

My goal is I want to crop this raster only for land coverage of Germany, so I did approximate projection in R, here is the possible R code I could give it try down below:

germany_territory <- readOGR("germany_landCover",layer="germany_landCover")     # read Germany NUTS3 level polygon (by district level)
proj <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84  +towgs84=0,0,0"
proj_raster_landcover <- projectRaster(land_cover, crs =  proj)

germany_landCover_city <- extract(x = proj_raster_landcover,                         y = germany_territory)

but these don't work for me, so I also tried something like this down below:

data(wrld_simpl)
germany <- wrld_simpl[wrld_simpl@data$NAME == "Germany",]
ex_deu <- extent(raster::area(land_cover))
germany_land <- crop(land_cover, ex_deu)

still can't achieve what I want.

desired output:

in cropped raster, I want a plot something like this down below:

I want to extract all information on agricultural and city land coverage only for Germany and render them in my final output as csv or xlsx format.

How can I easily deal with TIFF format raster in R?

How can I correctly crop the land coverage for Germany?

Any way to make this happen easily?

How can I correct my implementation on above?

Any more thoughts?


texts: SLD not changing style when publishing layer

texts: I've tried several times creating an SLD in QGIS and then bringing into GeoServer. However, after selecting the style in the Edit Layer-> Publish tab, it remains in the original default color.


texts: How to create point symbols with smooth color gradient fill and associated legend?

texts:
I like the feature heatmap as a layer style, but it doesn't have any symbolism in the legend at the compose manager. So I took a print screen from a characteristic part of the heatmap and now I am pasting it in the legend near the Text of the heatmap reference. But for every new project I have to make the same job, or move it if I move the legend...is there a simpler way for it?

I am using heatmap for stand alone trees instead of svg symbols or simple points. I think it makes the map prettier. So the heatmap has small radius circles and white area. A green square is not what I am looking for. Mostly I made this question, if there is somewhere a parameter for editing the point layer with the heatmap style and I haven't seen it.


texts: Difference between EPSG:4236 and EPSG:4326?

texts:
I know that the EPSG 4326 is the WGS84 (standard) projection, and is the projection that I am using in my code, but by a mistake, I have wrote 4236 and this worked in the code.

The question is what espg=4236 work? And which is the difference between 4326 and 4236?

The Python code:

geodf = gpd.read_file('./zone.shp')
geodf = geodf.to_crs(epsg='4236')
geodf.to_file('./zone.geojson', driver='GeoJSON')


texts: which spatial server goes better with raster images stored in DB Spatialite?

texts:
I would like to serve images from Spatialite as WMS services.

As far as I know:

Geoserver is compatible with vector features only -if you install the extension-, but not with spatialite BLOBs (in this case raster images).

Arcgis desktop since 10.2 supports vector layers in SQLITE, not raster images, what about ArcGIS Server, does it?

So, what other options do I have having RASTER spatialite support?


texts: Why getting wrong NDVI values resulted from scaled reflectance in ArcMap?

texts:
I have made atmospheric corrections to my Landsat 8-OLI raster using 'ENVI FLAASH' having surface reflectance values scaled between 1-100.

In ArcGIS I am using simple expression to compute NDVI using 'Raster Calculator' tool as illustrated below.

However, the results are not respectable. The same results are faced using 'ENVI Band Math' tool. Interestingly, this is not the case with the outcome from  'Spectral Indices' or 'NDVI' tools built within ENVI.

Why there is great discrepancy in results in all these cases?


texts: Forcing label under point shapefile in QGIS 2.18?

texts:
I need to force the labels of the street names to be under the main point shapefile of a map. As I understood it, this is not possible, but perhaps somebody has a workaround. For now it would be enough if it could work in the print composer. Here is a sample situation:

The red diamonds are the main information of the map and the street names should also be seen, but if conflicting with the point shapefile they should be under those diamonds.


texts: How to calculate past coordinate pairs from current coordinate pairs, a heading, and speed?

texts: I am very new to this side of the world. I am a classical embedded software developer. We want to incorporate GPS in our project. If I have the current Latitude, Longitude, Height,Heading (I hope this is what I think it is, the direction of the vehicle) and speed of the vehicle. Can I calculate the coordinates 20 meter before (I can also convert this to time scale). I am not sure how to proceed with this or what other information. I am not sure if this is totally possible at all.


texts:
What is the equivalent in QGIS (or other FOSS) to this:

rows = arcpy.SearchCursor(table,
                          fields="ORIG_FID; field1; field2; field3; field4",
                          sort_fields="ORIG_FID A")

for row in rows:


texts: Geometry Data in GeoDataFrame will not plot

texts:
I want to import a CSV to use as a geodataframe. When I plot it, nothing shows, which makes me think I'm doing something wrong. Here's my code.

import pandas as pd
from geopandas import GeoDataFrame
myfile = 'data.csv'
df = pd.read_csv(myfile)
crs = {'init': 'epsg:4269'}
gdf = GeoDataFrame(df, crs=crs)
gdf.plot()

The 'geometry' column in my CSV contains two rows, which contain LineStrings e.g.

LINESTRING(-117.674411 33.561825,-117.61164 33.589388)
LINESTRING(-117.687375 33.602688,-117.690752 33.604706)


texts: Openlayers 2: Check the map content and tell if the feature exists in the map bound

texts:
In my map there are a number of layers. I want my custom legends to show when the related layer feature is visible in the current map view. I am able to turn on/off the display of the legend according to display scale but I need to hide the legend if the feature is not in the map view even if it is in the range of display scale.

I go through the layer collection of the map. I tried to debug when the vector layer feature is visible in the map view but I get layer.features.length as always zero.

I need it total in terms of openlayers script. I do not want to utilize GeoExt.


texts:
I have a raster file on precipitation and I want to create zonal statistics as a table using a shapefile with boundaries of countries across the world. However it appears that I need to take into account that in different latitudes the  corresponding area is different and thus I need to weight my statistic (mean etc.) for area.

Any suggestions on how I can do that?


texts: QGIS 3.16.4 crashes after opening 3D map view

texts:
I have a problem with QGIS 3.16.4 3D Map Viewer.
When I try to open the viewer, QGIS crashes immediately. This happens both when I use only OSM as basemap, without data layer, but also when I visualize a simple point layer with 6 points. Does anyone have similar problems or a fix for this?
QGIS Info
QGIS Version: 3.16.4-Hannover
QGIS code revision: 654e76b
Compiled against Qt: 5.11.2
Running against Qt: 5.11.2
Compiled against GDAL: 3.1.4
Running against GDAL: 3.1.4
System Info
CPU Type: x86_64
Kernel Type: winnt
Kernel Version: 10.0.19041


texts: GeoNames cities list missing many cities that appear on website

texts:
I am trying to set up a table with cities.
I have downloaded the cities500, cities1000, cities5000 lists, but strangely enough they are missing rather big and popular cities/towns that ARE present on GeoNames website.

A couple of examples for these :

Ronda, Spain - https://www.geonames.org/2511730/ronda.html
Ribe, Denmark - https://www.geonames.org/2614813/ribe.html

Is there something I'm missing here? How come even though GeoNames does seem to contain those cities, the lists on the download server don't have them?

UPDATE :

I was using a bad CSV editor tool, both Ronda and Ribe do exist in the dump file.


texts: How to add labels using featurelayer and sceneview in ArcGIS API for JavaScript 4.3

texts:
I am using using featurelayer and sceneview in ArcGIS API for JavaScript 4.3 and want to add labels.  I have followed code examples from the API documentation but the labels do not show.  How do I get the labels to appear?

this.map = new Map({
  basemap: 'streets',
});

this.map.layers.add(new FeatureLayer({
  url: 'http://slcipuapp2.slcgov.com/arcgis/rest/services/waterhardness/MapServer',
  id: '0',
  visible: true,
  labelsVisible: true,
  outFields: ["*"],
  popupTemplate: this.locationsTemplate,
  labelingInfo: this.labelTemplate
}));

this.labelTemplate = new LabelClass();
this.labelTemplate.symbol = {
  type: "text",
  symbol: {
    type: "text",  // autocasts as new TextSymbol()
    color: "black",
    haloSize: 1,
    haloColor: "white"
  }
};
this.labelTemplate.labelExpressionInfo = { expression: "$feature.reading" };

this.view = new SceneView({
  container: "viewDiv",
  map: this.map,
  camera: {
    position: [
      -111.7496032288646, // lon
      40.478774371255994, // lat
      7671.145575070754  // elevation in meters
    ],
    tilt: 72.34558755016394,
    heading: -26.290004118126618,
    fov: 55
  }
});


texts: QGIS on MAC cannot be opened because of an error

texts:
I know this question has been answered before, 3 months ago, but cannot work out how to join that thread.

I have been using Mt Lion with QGIS successfully for some weeks but now get the  'Cannot be opened..' message on attempting to start.  I have reinstalled the GDAL framework Complete but without success.

Any suggestions most welcome!

TW

The core of the error message is:

Application Specific Information:
dyld: launch, loading dependent libraries

Dyld Error Message:
  Library not loaded: /Library/Frameworks/PROJ.framework/Versions/4/PROJ
  Referenced from: /Applications/QGIS.app/Contents/MacOS/QGIS
  Reason: Incompatible library version: QGIS requires version 8.0.0 or later, but PROJ provides version 7.0.0


texts: How create legend border same as map background like in the image? I'm just can get the border white like the rest of legend frame. I'm trying in Arcmap and ArcGIS Pro.


texts: OpenLayers: add a listener to a single polygon

texts:
I want to add a listener to a single polygon. Currently, I am trying to add one that makes a polygon opaque when you click on it, but I should be able to add any general listener. OpenLayers.Layer has an events property, but there is no such property for OpenLayers.Feature or OpenLayers.Geometry.Polygon.

This is my setup:

An OpenLayers.Layer.Vector layer, pLayer, called "polygon layer" that contains all of the polygons.
A function handleClick that sets this.style.fillOpacity = 1. this is supposed to refer to the polygon clicked on
For each polygon, I call:

pLayer.events.register('click', polygon, handleClick);

The issue with this setup is that when I click on any of the polygons, it fires the event for all polygons on that layer and every polygon becomes opaque. How do I make it so that the event only fires for individual polygons?

I need this to be general. Basically, I want to be able to add any type of listener to a single polygon. My overall goal is to create an interface between Google Maps and OpenLayers, and I Google Maps you can add an event to a single polygon.

Will I need to put every polygon on a separate layer? If that is the case, how do I set up the renderer for each of those polygon layers? Or is there a way to keep all polygons in the same layer and have the event handle each of them individually?


texts: QGIS in Qt Widgets (Cross Platform Application)

texts:
I am currently developing a widget Application in Qt, and I need to use GIS in this software.

Naturally I first found the ArcGIS API for Qt, but like ArcGIS, it only runs on Windows and Linux, and not MacOS.

I then thought of QGIS, as it runs on all three operating systems. Yet, I cannot find any online resource that tells me how to integrate QGIS in my application for all three platforms. I only found a few questions here about windows and linux, but nothing about MacOS, nor how to integrate it for all operating systems.

If this is possible, are there resources online that explain how to do this?


texts: How to output file with all overlaps between two Shapefiles - Neighborhood/Zip Codes

texts:
I have a pretty specific problem I would like help with.

I have two shapefiles:

neighborhoods from https://www.zillow.com/howto/api/neighborhood-boundaries.htm

zip codes from the census bureau: https://www.census.gov/geo/maps-data/data/cbf/cbf_zcta.html

The end result I would ideally like to achieve is this: a csv with each neighborhood in a specific file from Zillow (let's say we use the NY file) along with all zip codes from the second shapefile that overlap with that neighborhood.

For example:

Essentially what I would like to achieve is https://www.health.ny.gov/statistics/cancer/registry/appendix/neighborhoods.htm

but for all other major cities in the USA, using the Zillow mapping.

Seems like the way to approach this is using Python with Shapely/Fiona, but wondering if there is an elegant solution.


texts: Getting Node.appendChild error with leaflet

texts:
I am trying to add a popup function to some polygons. I've been consulting a few different tutorials trying to work it out and I'm concerned that I am getting conflicting instructions.

The map I have now is displaying perfectly. The next step for me is to add the popup when the polygon is clicked on. I keep getting the error "Argument 1 of Node.appendChild is not an object." I have no idea what this means. My only guess is that it has something to do with using a var within landscape.js to create the landscape.geojson, since that might not be creating it in the right way to use it for the popup.

<link rel="stylesheet" href="./leaflet/leaflet.css" />


texts: Filtering Geoserver WFS by using a reverse bounding box

texts:
So I have a WFS being served through Geoserver that links to a Postgis table.

It's a very large table so up until now I've been checking the box that says "Only request features overlapping current view extent" but the problem is that the data sometimes keeps going on a loop and it's impossible to work in some areas

The workaround I thought about was to use a bounding box for the area with very high frequency of events and another for everything outside of it. How to do this?


texts:
I was wondering what are the necessary files for hosting a webgis server with OpenLayers3 without using a CDN?

I am testing on a localhost without internet access. It  doesn't work to merely copy the direct js (and css) file pointed to by a CDN such as

http://openlayers.org/en/v3.0.0/build/ol.js

to the local webserver as <script src="js/ol.js"></script>.

Also, the same question holds for jquery and bootstrap. I had

<script src="http://code.jquery.com/jquery-2.1.1.min.js"></script>
<script src="https://maxcdn.bootstrapcdn.com/bootstrap/3.2.0/js/bootstrap.min.js"></script>

But I wasn't sure how to get a complete local version for those.


texts: How to zoom in/out selected elements of MapControl?

texts:
I have a mapcontrol that includes many elements: text, line.
I used this code to add elements into mapcontrol:

 pGraphicsContainer.AddElement(pElement, 0);

I can use zoom button to zoom in/out the mapcontrol, but how can I zoom in/out only seleted elements ?

Is there any buttons for this purpose in Arcgis ?

Thanks and regards,

Tai


texts: QGIS Process Modeller Multiple Input (3.4.4)

texts:
I am trying to build a process with the process modeler.
The process is working with one shapefile via Input vectorlayer.

But in the final Version the model should do a Folder or Muliple Input. But if I Chose Folder or Multiple Input, I can't connect them to an other part of the model.

Could please anyone help me?


texts:
I was looking into PostGIS and its defined Spatial Reference Systems, and I found that some have no proj4.
Why don't they?

SELECT auth_name || ':' || auth_srid, proj4text FROM spatial_ref_sys ORDER BY proj4text

  "EPSG:22300";""

  "EPSG:29701";""

  "EPSG:3577";"+proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132
  +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs "


texts: ArcMap suddenly started exporting massive .png files?

texts:
I'm exporting a very simple map. It has a single polygon shapefile on top of ESRI base imagery. There are several text boxes as well, but just simple Arial text.

I'm exporting to .png at 220 DPI and the file size is huge - 19MB. I've exported many similar maps and they are usually between 1 and 4MB. Any idea what might be going on?


texts: Adding KmlLayer using local file (File API)

texts:
I am trying to load a local KML file to the map for this I am using FileReader object (File API). To use google.maps.KmlLayer I need to use readAsDataURL to transform to a URL the file, the problem is that KML File is not loading to the map and I don't get any error. I think may be url attribute in KmlLayer doesn't accept data url.

This is the code:

gp.addEventListeners = function(){
  $("#file-upload").change(function (evt){
      var file  = evt.target.files[0]; // FileList object
      var reader = new FileReader();

      reader.onloadend = function () {
         var ctaLayer = new google.maps.KmlLayer({
          url: reader.result
        });
        ctaLayer.setMap(gp.map);
      }

      if (file) {
        reader.readAsDataURL(file);
      }
  });
}

In OpenLayes the same code is working.

What I can do? I need help :(

Sorry for my English.


texts: Defining a color ramp for temperature range & applying it to different rasters?

texts:
I'm using QGIS 2.18 (Las Palmas).

I have downloaded Climate-Data from WorldClim and I would like to display all of the monthly maps within the same color scale (and I don't mean within Styles always choosing the same color ramp because if I choose a specific color ramp for January where the Temperatures vary between -40 and -10 and in February the temperatures vary between -20 and +5, the maps will look similar). This has the effect, that when clicking through the maps, the monthly changes become evident, even without having to look at the temperature legend).

Is there a way I can create a color scale for my temperatures between -40 and +20°C although Not all my maps vary between these temperatures?


texts: How to spatially join a grid with a landuse layer in QGIS?

texts:
I'm new to QGIS and I'd like to know if exist in QGIS the equivalent command for "spatial join" of ArcGIS. I have tried with "join attributes by location" but seems to be impossible setting the match options. I have to find the correspondence between each element of a grid and a vector layer which represent the land use.

I'm dealing with two vector layers. What I have to do is associate to each polygon of the grid,the value of the class land use which cover the specific area. With "match option" I mean for example, the possibility to associate to each element of the grid only the class that fall in the centre of the area. (The land use is a v polygonnal vector layer)


texts: Why does ArcPy plot just one point when there are many?

texts:
Why does following ArcPy code, when run in Python window in ArcGIS Pro, plot only one point, the last one in for loop, that is (538973.76, 6478442.80)?
My aim was to get all the points (x,y) plotted.
import arcpy
from arcpy import env
env.outputCoordinateSystem = r"C:\temp\coordsys.prj"
env.workspace = r"C:\test.gdb"

x = [538993.76, 538983.76, 538973.76]
y = [6478462.80, 6478452.80, 6478442.80]

for i in range(len(x)):
   point = arcpy.Point(x[i],y[i])
   pointGeometry = arcpy.PointGeometry(point)
   arcpy.CopyFeatures_management(pointGeometry, "Hydrant")


texts: How to create non-overlapping intervals in the legend?

texts:
So I am trying to create a legend where the entries are, for example, like this:

1.1 - 2.0
2.1 - 3.0
3.1 - 4.0

as opposed to like this:

1.1 - 2.1
2.1 - 3.1
3.1 - 4.1

I essentially don't want the upperbound of one interval to be the same as the lowerbound of the next interval. I know I can do this by setting the Legend Format to "%11 - %20" but in that case, I would end up with two digits after the decimal, where as I only want one digit after the decimal.


texts: Invoking QGIS style dialog from Python standalone script

texts:
I can load a layer, open a canvas, and draw the layer in the canvas using the Qgis python API.

Can I bring up the properties dialog, or even just the "Style" tab of the properties dialog, from such a standalone application? I have a bad feeling I can't because all the code is in the main Qgis Application rather than present in qgis.core and qgis.gui, and doesn't have a python interface.

Unless I've missed something, the solution might be to rebuild the dialog and its behaviour in Python, re-using the relevant .ui file?


texts: Reading a shapefile as an array using Python?

texts:
I am trying to read a shapefile as an array with Python but I keep getting the same error:

Warning 1: Failed to find field ID on layer shapefile_maido_tipe, skipping.

ERROR 4: `temp.tif' not recognized as a supported file format.

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Anaconda2\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", line 714, in runfile
    execfile(filename, namespace)
  File "C:\Anaconda2\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", line 74, in execfile
    exec(compile(scripttext, filename, 'exec'), glob, loc)
  File "C:/Users/python/Desktop/programme_python_tipe/test.py", line 38, in <module>
    print layer('shapefile_maido_tipe.shx')
  File "C:/Users/python/Desktop/programme_python_tipe/test.py", line 36, in layer
    return gdal.Open('temp.tif').ReadAsArray()
AttributeError: 'NoneType' object has no attribute 'ReadAsArray'

Here is the code that I got from Convert polygons in shapefile to a NumPy array? :

import gdal
from osgeo import osr
from osgeo import ogr

def layer(shapefile):

    # 1) opening the shapefile
    source_ds = ogr.Open(shapefile)
    source_layer = source_ds.GetLayer()

    # 2) Creating the destination raster data source

    pixelWidth = pixelHeight = 1 # depending how fine you want your raster
    x_min, x_max, y_min, y_max = source_layer.GetExtent()
    cols = int((x_max - x_min) / pixelHeight)
    rows = int((y_max - y_min) / pixelWidth)
    target_ds = gdal.GetDriverByName('GTiff').Create('temp.tif', cols, rows, 1, gdal.GDT_Byte)
    target_ds.SetGeoTransform((x_min, pixelWidth, 0, y_min, 0, pixelHeight))
    band = target_ds.GetRasterBand(1)
    NoData_value = 255
    band.SetNoDataValue(NoData_value)
    band.FlushCache()

    # 4) Instead of setting a general burn_value, use optionsand set it to the attribute that contains the relevant unique value ["ATTRIBUTE=ID"]

    gdal.RasterizeLayer(target_ds, [1], source_layer, options = ['ATTRIBUTE=ID'])

    # 5) Adding a spatial reference

    target_dsSRS = osr.SpatialReference()
    target_dsSRS.ImportFromEPSG(4326)
    target_ds.SetProjection(target_dsSRS.ExportToWkt())

    return gdal.Open('temp.tif').ReadAsArray()

print layer('shapefile_maido_tipe.shx')

# maido is the name of a mountain
# tipe is the name of a french school project

Can you tell what am I doing wrong please ?

In this link you will find the shapefile I use: https://drive.google.com/drive/folders/0B7P95aWmH4DUYUZfcTVvbzF4bG8?usp=sharing

PS: The Reunion Island National Park is the owner and also the creator of the shapefile. I have to mention them every time I communicate this file.


texts: Memory efficient shortest path algorithm for offline Android app

texts:
Finding the shortest path in a road network requires the road network to be loaded up in memory. Using Dijkstra's algorithm, for example, it searches through all the immediate road coordinates and puts a weight on that road based on the distance and time to pass through it. This keeps on going until it reaches the destination coordinate which is when Dijkstra is able to find the shortest path.

This approach is very memory intensive since the whole road network needs to be in memory.

I would like to do some more research and look at the literature when it comes to finding the shortest path on an offline mobile app with limited memory.

Since I am quite new to this finding the appropriate literature to get acquainted with is a problem in itself.

I'm just looking for pointers.


texts: How can I ignore suggested username for new members

texts:
I'm trying to create a couple of new members for our organization. The usernames must match a series the usernames that exist in another app. The problem is when I enter the username ArcGIS Online automatically overwrites the entered username with another suggestion, usually by appending numbers to the end.

So for example I want to create 'RoadCrew1' but ArcGIS Online overwrites it to 'RoadCrew14' as soon as I click out of the username text box. I cannot turn off these suggestions. I suspect that these usernames existed in our organization before but if so have been deleted. Is there any way to get these usernames back?


texts: How to display shapefiles on an OpenLayers web mapping application that are provided by users as an input without making any changes in the database

texts: I have created a web mapping application using OpenLayers 3 which is displaying the data that is stored in a database. Now I want to display additional data that is going to be provided by the users in the form of shapefiles. The data should only be displayed on the main application and not stored in the database as I do not want to make any changes in the database. I know it can be done through WMS but how it can be done without altering the database?


texts: Leaflet Marker Cluster: Popups won't open after first

texts:
I built a map using Leaflet and markercluster.js, and I'm having a strange issue where after I click on the first popup, I'm unable to open any other popup. This happens whether you close the first popup manually or not.

When I don't use markercluster.js, the popups appear as normal, so I know it has something to do with that library. This is the second time I have used this library; the first time I had no issues. I'm getting no errors in the console.

Also, I had this issue regardless of if the .openPopup() method was attached to the marker variable.

Code below:

  var puetrasList = new L.MarkerClusterGroup();

          for (i=0 ; i < dataStore.length ; i++){
              var a, point;
              a = dataStore[i];
              lat = parseFloat(dataStore[i][0]);
              lng = parseFloat(dataStore[i][1]);
              title = "
" + dataStore[i][2] + "
" +
                                  "
" + >dataStore[i][3] + "
" +
                                  "
" + >dataStore[i][4] + ", " + dataStore[i][5] + "
" +
                                  "
" + >dataStore[i][6] + "
";
              marker = new L.Marker(new L.LatLng(lat, lng));
              marker.bindPopup(title).openPopup();
              console.log(marker)
              puetrasList.addLayer(marker);
          }

          //----- Summon Map ---->
          map.addLayer(puetrasList);
          map.addLayer(layer);


texts: add-in to open automatically "data frame properties

texts: How can  created an ArcMap Python add-in using ArcGIS 10.1 to open  automatically data frame properties And change the background color of another?


texts: Getting elevation heights from folder of thousands of shapefiles that apparently come from LiDAR?

texts:
I obtained a few folders purporting to come from LiDAR data taken of a U.S. county. The directories are set up like this.

Data-folder
- be_rasters
-- 692500_725000
--- dblbnd.adf
--- hdr.adf
--- metadata.xml
--- prj.adf
--- sta.adf
--- w001001.adf
--- w001001x.adf
-- (thousands more folders with this naming convention format and files within them named the same thing)
- contour
-- shp
--- Contour_692500_735000.cpg
--- Contour_692500_735000.dbf
--- Contour_692500_735000.prj
--- Contour_692500_735000.sbn
--- Contour_692500_735000.sbx
--- Contour_692500_735000.shp
--- Contour_692500_735000.shp.xml
--- Contour_692500_735000.shx
--- (Thousands more files with similar naming conventions, with the same extensions)

It appears that each shapefile contains data for a small area of the county. I want to obtain the average elevation above/below sea level for each section of the county.

How do I obtain that information from these files?

I have QGIS. I tried opening a few .shp files in it. But trying to open all the .shp files at once freezes the program.


texts: Updating datefield from a string with UpdateCursor

texts:
I have a string I'm slicing and then using UpdateCursor to update rows with the outputs. The string is formatted like this:
Foo_YYYYmmddHM_Xy_Ab_Bar
I'm struggling in getting only the YYYYmmdd part of the string in a date field. My current approach is:
with arcpy.da.UpdateCursor(dataset, ["String", "Year_short","Date_date"]) as updateRows:
    for row in updateRows:
        string_split = row[0].split("_")
        row[1] = datetime.datetime.strptime(string_split[1], '%Y%m%d%H%M').year
        row[2] = (datetime.datetime.strptime(string_split[1], '%Y%m%d%H%M')).strftime("%Y%m%d")
        updateRows.updateRow(row)

row[1] updates OK because the field is a short type, however no matter what I've tried row[2] will fail unless I include the HM component part of the date. Above is my last attempt but I have also tried the following with no success:
   row[2] = datetime.datetime.strptime(string_split[1], '%Y%m%d%H%M').date
   row[2] = datetime.datetime.strptime(string_split[1], '%Y%m%d%)


texts: tmap::tm_facets() in interactive mode in html R Markdown file, small multiples variable title renders in wrong place

texts:
If you look at this reproducible snippet, when knitting, the first small multiples set renders both its and the next small multiples' title variable in the first set of maps.

Here's the reproducible example:

---
title: "TEST"
author: "author"
output:
  html_document
---

```{r}

library("tmap")
library("tmaptools")

tmap_mode("view")
data(World)

tm_shape(World) +
  tm_fill("pop_est", thres.poly = 0) +
  tm_facets("name", free.coords=TRUE, drop.units=TRUE) +
  tm_layout(panel.label.size = 2, legend.show = FALSE)

tm_shape(World) +
  tm_fill("pop_est", thres.poly = 0, palette = "BuGn") +
  tm_facets("name", free.coords=TRUE, drop.units=TRUE) +
  tm_layout(panel.label.size = 2, legend.show = FALSE)

```

As you can see, the first set of multiples show the country names twice, whereas they're not rendered in the second set at all. This is showing itself in a more complicated report I'm working on elsewhere, but this is reproducible using nothing but tmap data and functions.

Is this a bug or am I doing something wrong?


texts: Can a PostGIS table contain features with both single-point and multi-point geometries?

texts:
I have a geopackage containing point features. Most features are single points but some are multipoints. When I attempt to import to PostGIS using ogr2ogr, I received the following error:

ERROR:  Geometry type (MultiPoint) does not match column type (Point)

I then attempted to create a table beforehand where the geometry column is defined as multipoint. When attempting to load the geopackage into this table, I received this error:

ERROR 1: ERROR:  Geometry type (Point) does not match column type (MultiPoint)

Is it thus impossible to have a PostGIS table containing both single points and multi-points? Can I convert the singular point features to multi-points?


texts: How to get the click coordinates for the use in a custom form

texts:
I'm having a Spatialite layer to store photo locations from objects (stored in another Spatialite layer) taken in the field. Therefor I'm using a custom photo form like explained here.

Since now I type the corresponding object id in the "Objekt" field as an attribute of the photo layer. As an improvement I want to have a combo box (see picture of the photo form) were I can choose the object ids from a selection of the 5 nearest objects.

What I need is the coordinates from the point that will be created after I close the form. Or in other the words: the coordinates from the canvas click before the custom form appears. To search the 5 nearest objects (which is simple to achieve) I need the coordinates before the custom form opens, then a function can search for the neighbors and the results can be added to the combo box.

By now I'm using the following map tool (for the use in the Python console) to get click coordinates from the canvas:

from qgis.gui import QgsMapTool
from PyQt4.QtCore import pyqtSignal

class PointTool(QgsMapTool):
    def __init__(self, canvas):
        QgsMapTool.__init__(self, canvas)
        self.canvas = canvas

    canvasClicked = pyqtSignal(['QgsPoint'])
    def canvasReleaseEvent(self, event):
        point = self.toMapCoordinates(event.pos())
        self.canvasClicked.emit(point)

def getCanvasCoordinates(point):
    QgsMessageLog.logMessage(str(point))

mapTool = PointTool(iface.mapCanvas())
iface.mapCanvas().setMapTool(mapTool)
mapTool.canvasClicked.connect(getCanvasCoordinates)

How can I implement this in the custom form to get the coordinates, that I need to make the next steps explained above?

Is there a way to connect the getCanvasCoordinates function to the QAction like:

iface.actionAddFeature().triggered.connect(getRasterValue)


texts: Equations to convert from global cartesian coordinates to geographic coordinates

texts:
First, I needed to convert geographic coordinates (like those received from GPS) to global cartesian coordinates. I found an earlier question (How do I calculate a xyz-position of a gps-position relative to an other gps-position?) which worked perfectly. Now I need to go the other way around.

In Wikipedia, there are equations to convert between those, but the equations from geographic to cartesian are different, and not indicating which is latitude and which is longitude. I know it's because it's not geographic per se, and I also know that those angles have different constraints than those of latitude/longitude. This makes me not trust the reverse equations there. I searched elsewhere but I couldn't find the answer.

Precision is not important. Consider a perfect sphere.


texts: In OpenLayers zoom to extent of all overlays

texts:
In Openlayers I can zoom to the bounding box of all layers like this:

map.zoomToExtent();

I can also zoom to the bounding box of a single layer:

map.zoomToExtent(layer_name.getDataExtent());

Is there an easy way of zooming to the extent of all layers except the base layer? In other words showing the extent of all overlays.


texts: Change the symbol for a single feature in ArcObjects?

texts:
Is it possible to change the symbol for a single feature in ArcObjects? If so how?

I know it's possible to setup the map symbology for different field values of a feature class, but I'd like to achieve this without changing any of the underlying field data if possible.

P.S. I'm using C# and ArcGIS 10.0


texts: How to create a new value with information from multiple columns?

texts:
I have dataset that has columns that need to be sorted based on multiple columns.

The first column contains a name and the second a number, the names are not unique. What I want is a NewID from 1 to 10 per name based on the descending value of the number.
Like this:

Any ideas how this can be done? Preferably in ArcMap.


texts: Debugging CalculateField_management ERROR 000623: Invalid value type for parameter field?

texts:
I am trying to select records in my table that have certain words in the "Name" field and then code them as 0 (not pharmacy) or 1 (pharmacy)  This is just a test searching for two words to select non-pharmacies.

I am getting a:

  ExecuteError: ERROR 000622: Failed to execute (Calculate Field).
  Parameters are not valid. ERROR 000623: Invalid value type for
  parameter field.

I don't know which variable is invalid.  Any thoughts?  I also do not know if the Calculate Field will work on only the records selected.  I cannot get far enough to check that.

# import system modules
import arcpy
import random

# set workspace environment
myWorkspace = arcpy.env.workspace = "C:/PharmacyProtocol/tables"
arcpy.env.overwriteOutput = True

# set variables
out_table = "VT_TRL_DLC.dbf"

# set field variables
fieldList = arcpy.ListFields(out_table)

for field in fieldList:
    if field.name == "Name":
        nameF = field
    elif field.name == "Pharmacy":
        pharmF = field

#create list of words that are NOT pharmacies
noPharmWords = ['AMERICAN LEGION','BAR']

expression = 0

#find words in Name field, code Pharmacy as "0" if not already coded
for codeWord in noPharmWords:

#create SQL text clause
fldName = '"Name"'
endClause = "%'"
SQLstr = fldName + " LIKE '%" + codeWord + endClause

#select records from table with codeWord
cursor = arcpy.da.SearchCursor(out_table, ["Name"], SQLstr)

arcpy.CalculateField_management(out_table, pharmF, expression)

del cursor


texts: Ocean flow models - predict the path of an object dropped in the ocean

texts:
This question is from a user who wants to calculate the predicted path of a bottle dropped from a ship. He has a spreadsheet of hundreds of these bottles and wishes to display the predicted vs actual landing points of the bottles.

Is there a publicly-accessible model which will show the predicted path of an object, anywhere in the earth's oceans?

For example, NOAA published a series of maps during the Gulf of Mexico oil crisis, showing the predicted path of the oil spill.

I'm presuming that this is a fairly specialised branch of GI-science and the answer will be that such a model depends highly on local variations, the date, wind and other weather conditions, etc.

As a next-best option, is information on the general patterns of the ocean currents publicly accessible in a GIS-ready format? It may be possible to create a geoprocessing model from this?


texts: Merging or dissolving points while maintaining attributes using ArcGIS Desktop?

texts:
I have a point shapefile with points that have identical geometry but different attributes. I would like to merge the points in the same location, but keep the highest value of the attribute from them.

I tried using dissolve and statistics fields, it dissolved the points but it kept the highest value of each attribute from the entire shapefile, while i only want the point to keep the value of each cluster.

For example, in a position i have a pole, with 3 fiber optics lines coming, and 2 leaving. I would like for the position of the pole to merge the points into one, and keep the highes value, in this case, 3.


texts: Jagged line and markers when publish map on ArcGIS Server

texts:
We developed a native app with Esri iOS SDK and Esri Android SDK using ArcGIS Server, when we publish our map using Server, road lines become jagged. So when using this road service in our app development, app visual become so terrible.

How can we smooth the lines or disable jagged display of roads?

Not: You can see real image at this link : http://i.stack.imgur.com/sEPuw.png


texts: Is it possible to use symbols instead of colours to visualise raster cell values?

texts:
I have a raster image (GeoTIFF) with integer values. I know how to visualise raster values in QGIS using colormaps or similar, but I would like to assign symbols to each colour value instead - just like a colormap, but with symbols.

Example:

Question:
Is this possible and in that case how?


texts: Error converting vector shapefile to raster (.TIF) in QGIS?

texts:
I am working with a shapefile of 10km grid regions of the UK obtained from digimap.

I need to convert this to a raster in order to do some analysis in R regarding agricultural intensity of 10km UK sites.

Everytime I load the vector file in QGIS and do raster->conversion -> rasterize I get an error which reads:

The following layers were not correctly generated./private/var/folders/q3/vldvgpt970g5mn1y_46n4jjr0000gn/T/processing_cf35adc2256045b9b004c26c6ce46245/94b6409973714574a8226988cee87064/OUTPUT.tifYou can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.

I know there are a lot of other questions about this but none of the answers have helped me to solve what I am trying to do.

Perhaps I am not inputting the correct parameters in the dialog box, could anyone help me with this? I've tried using layer extent and canvas extent in the output extent box and both have not worked.


texts:
I am looking for a way to create a multi distance buffer and have each buffer a different colour. For example, a buffer around a point has 5 rings. I would like to have each ring a different colour fill.

Is there a simple way of doing this?


texts:
I have a point shapefile that is symbolized with an svg.  The shapefile has a attribute with Y values.  I would like the points with lower Y values to draw on top of points with higher Y values.  Below is a SE link with instructions for an earlier version of QGIS that has that feaure called Define Order.  I do not see that as an option in v2.8.  How to change feature draw order in QGIS?.

Is that feature still an option in QGIS?


texts: Does a national(USA) interstate rest stop/rest area dataset exist?

texts:
And I don't care about format--worry about that later!  I found some maps in Google maps, usually per state.  I also have one map made my someone else, that no longer appears in the Google Maps search for user content.  Maybe that user made it private, but I still can see it.

The reason I ask is related to my google maps navigation question on android.se.  I think the easiest way to get this info on the droid is to create or use a shared map in Google Maps.  I can find Florida, but no others show up.

Maybe it's a vocabulary issue?  What else are are rest stops called?