texts: SAGA Inverse Distance Weighting in QGIS - How to Set Minimum Number of Points?

texts:
I am analyzing a postGIS point dataset in QGIS version 2.8.2, about 200k features total. To calculate weighted average values I am using the Inverse Distance Weighted interpolation technique. However, I am concerned that the "high" values I'm seeing in the outputs are the result of a handful of outlier values, not the basis of any meaningful spatial pattern.

Is there a way in the SAGA IDW plugin - or another IDW tool - to set a minimum number of points required for the raster surface to calculated the weighted averages, the same way it considers a maximum number of points?


texts: ArcObjects Load features that don't exist

texts:
I'm trying to write a method in ArcObjects .NET that will load features from a file GDB into another geodatabase but only if they don't otherwise exist.  Records are considered existing if there is a matching GlobalID.

Without looping through record by record and setting up a query for each record, is there a more efficient way to do this?


texts:
I have a WFS server that I need to filter information out of with filter encoding. If the attributes are on the surface, I don't have any problems getting their information. But I don't know how to get the information if it is a subfeature.

In the example below, I would use:

<fes:PropertyIsEqualTo>
  <fes:ValueReference>app:dørtype</fes:ValueReference>
  <fes:Literal>Slagdør ut</fes:Literal>
</fes:PropertyIsEqualTo>

to get all with that value, but I don't know how I would get "app:rampeLengde" under "app:Rampe".


texts: Using GIS to determine Average Land Distance From Equator?

texts:
How far away is the average piece of land on Earth from the equator? Put more precisely, what is the median absolute value of the latitude of landmass on Earth?

From spherical geometry, the median absolute value of latitude on the Earth's surface is 30 degrees. However, land might tend to be further or closer to the equator than this.

Note: I originally posted this in the Earth Science Stack Exchange, but this question is a better fit here instead.


texts: Change stroke width from Milimeters to map units

texts:
I would like to change the stroke width from mm to map units. The code below is mainly taken from ThomasG77's answer from a previous question. I only found a way to change the units of Labels unfortunately.
Use colums in attribute table to set colors and line thickness directly automatic
layer = iface.activeLayer()

# Get colors in rgb codes from Attribute table
layer.renderer().symbol().symbolLayer(0).dataDefinedProperties().property(QgsSymbolLayer.PropertyStrokeColor).setExpressionString("\"Farben\" || ',255'")
layer.renderer().symbol().symbolLayer(0).dataDefinedProperties().property(QgsSymbolLayer.PropertyStrokeColor).setActive(True)

# Get Stroke width from Attribute table
layer.renderer().symbol().symbolLayer(0).dataDefinedProperties().property(QgsSymbolLayer.PropertyStrokeWidth).setExpressionString("\"Linien\"")
layer.renderer().symbol().symbolLayer(0).dataDefinedProperties().property(QgsSymbolLayer.PropertyStrokeWidth).setActive(True)

# Refresh rendering of the layer
layer.triggerRepaint()


texts: Create line segment at point perpendicular to line

texts:
I have a long road with 50+ points on it to mark the drainage pipes that go under the road. However, the pipes are not points, they are long cylinders. They are (almost) always perpendicular to the road. Is there a way to turn each point in to a 20 meter line which is perpendicular to the road? I would at the same time like to preserve some of the attribute fields in the point, where the diameter of the pipe is specified. This can of course be done using a spatial join later.

I have tried creating the line features myself in Editor mode, using the Constrain Perpendicular function on the Feature Construction toolbar. But doing this manually takes too long (there are more roads with even more drainage pipe points).

The pink lines are how I imagine the output to be:

I am using ArcGIS Desktop 10.6 and ArcGIS Pro 2.4. I am not yet well-versed in Python.


texts: QGIS buffer polygon with different widths on each side

texts:
I create maps for wildfire risk assessments and am constantly asked to plot "defendable space" buffer zones around building envelopes that are different distances on different sides of the buildings. (See image)

I have been doing this in QGIS by creating a buffer for each distance required and then manually drawing one that matches as close as I can.

This is obviously less than ideal. Is there a way I can create these buffers with QGIS? I also have access to ArcGIS and MapInfo.


texts: Resample function doesn't work for Landsat-8 NDVI 8-days product in Google Earth Engine

texts:
I am using Landsat 8 Collection 1 Tier 1 8-Day NDVI Composite in Google Earth Engine and I want to export it with the projection same as code editor (EPSG:3857). However, when I add the reprojected image in the map, it appeared very weird.
The original image looks like:

but after reprojection the image looks like:

I also tried to use Landsat 8 Collection 1 Tier 1 8-Day TOA Reflectance Composite to directly calculate NDVI and reproject it, but it still failed.

Another problem is the NDVI calculated by image is different from products provided by Google, is there any differences in calculating NDVI?

The code is here:

https://code.earthengine.google.com/a88b1b6da7529e3a10e642e28002f4f7


texts: Shifting a set of points and get new coordinates?

texts:
I have a transect that I've added as a delimited text layer that has UTM X and Y coordinates for about 300 points all within a river meander.

 They currently do not line up with the channel shown in the raster, so I would like to use the selection tool to highlight and shift them, just like the first ~1 min of this video for QGIS 2.8.2. https://www.youtube.com/watch?v=yiihopTnewg

When I select this layer, I am not able to toggle editing. I cannot seem to find the reason why. Once I do move these points, will I be able to export the new coordinates? Is there a better way to do this, considering I have about 20 transects to shift from 0.4-3m?
Example of coordinates:

(hi I'm very new, apologies if this should be obvious)

edit: I'm using QGIS 3.6.


texts: Get the point coordinates of images of an interactive web map?

texts:
I want to parse this web map for the point coordinates of the displayed icons:
https://www.mercedes-benz.de/passengercars/mercedes-benz-cars/dealer-locator.html

The HTML class of the items is "leaflet-marker-icon leaflet-zoom-animated leaflet-clickable"

Is there a way to find out the lat/lon-coordinates of every displayed item in the map?

I have not found any helpful workarounds.


texts: Unable to read TIFF file with raster2pgsql

texts:
I am having some trouble with the raster2pgsql function (on a Windows 7 in the OSGeo4W shell). I am receiving the error "Unable to read raster file". I can import shapefiles, ok using org2org so I know that the connection works, it just does not work for the raster2pgsql function.

To retrace my steps, I first changed the directory from m:/> to the folder containing my raster:

cd M:\Users\userID\Documents\R\Koala_Proj\koala_data\koala_information_base\habitat_suitability

I then check the files within the directory using >dir to see what the raster is called. From this, I receive a list of the one raster in this folder and the affiliated files, such as:
koala_raster.tif

I then try to use the raster2pgsql to import this raster to the above directory in postgresql:

raster2pgsql -s 4283 -I -C -t 100x100 koala_raster.tif public.khsm │
psql -h localhost -p 1234 -d koala -U postgres

I have also tried:

raster2pgsql -s 4283 -I -C -t 100x100 *.tif public.khsm │
psql -h localhost -p 1234 -d koala -U postgres

From these queries, I receive the error "Unable to read raster file".

However, rather than the file that this error refers bring named "koala_raster.tif" it prints "._koala_raster.tif". I am not sure what the "._" at the beginning of my raster refers to, or how to fix this problem I am having.

I have looked for fixes in

Uploading raster format(*.tif) to postgis through raster2pgsql
Cannot upload an IMG file into PostGIS using raster2pgsql
and the help files

However, I still have not found a solution that fixes my problem.


texts: Incorrect values of NDSI index after TOA correction

texts:
I would like to calculate NDSI index (Normalized-Difference Snow Index) which basically is the same as NDVI but presents snow not vegetation :). As in NDVI the range of values for NDSI should be <-1, 1>. However, after application TOA correction values of NDSI extend range <-1,1> and it's difficult to find treshold for snow because most of the NDSI values are cummulated in range <0.9,1>. In publications I've found the treshold value as ~0.4 which is almost impossible to find on my images.

I calculated NDSI also for DN values and everything seems to be fine - range of NDSI values is <-1,1>, treshold ~0.4. But I would like to apply TOA correction because in next stage I'll compare NDSI images from different acquistition time.

I've tried different software to do the analysis: PCI Geomatica, ArcMap, RStudio... Results are the same...


texts: How to compute drive access to a station based on vehicle origin points?

texts:
I have a shapefile that contains station locations and another shapefile that is the geocoded origins of those who drove to the station. I would like to make a buffer around the station that shows the area where 90% of the vehicles came from so I can show the general drive access area for each station without the extreme outliers.

Again, from my current files I would like to create a raster or shapefile that shows the general catchment area of each station.


texts:
I am currently working with QGIS 2.4 on Windows 8.1 pro and with QGIS-server on Debian 7. I want to publish a map with diagram but the render is not good.
At right is the desktop map and at left the server version. The circular diagram stays black on the server version.

Is anyone already uses and success to display diagrams and histograms with qgis-server ?


texts:
I have a problem with IDW interpolation. I need to use canvas extent and other attributes (I created an option with only attribute 0, other attributes aren't working). The documentation is not very helpful.

My code:

layer = iface.activeLayer()

layer_data = QgsInterpolator.LayerData()
layer_data.source = xlayer
layer_data.zCoordInterpolation=False
layer_data.interpolationAttribute =0
layer_data.SourceType =SourcePoints
layer_data.distanceCoefficient=2

tin_interpolator = QgsIDWInterpolator([layer_data])

export_path = r"C:\Users\Maria\Documents\miau.tif"

cos= canvas.extent()

output = QgsGridFileWriter(tin_interpolator,export_path,cos,2000,2000)
output.writeFile()

rlayer=iface.addRasterLayer(export_path, "interpolation_output")


texts: How to view or export WMS data in tabular format?

texts:
Is there any way to export or view data from a WMS in tabular format?

The info tool brings up a table showing the attributes for the selected feature, but I want to view attributes for all features.


texts: Convert coordinates from readShapePoly in R to long-lat coordinates

texts:
I'm trying to find some centroids of a shapefile for Danish municipalities and afterwards find the driving time between them. I use R's readShapePoly function from maptools combined with the gCentroidfunction from rgeos, and everything works. However, I get spatialpoints such as

SpatialPoints:
     x       y
1 571860.7 6225016
Coordinate Reference System (CRS) arguments: NA

Which is clearly not something I can use in Google to grab travel times. I'm looking for a way to convert these numbers to longitude-latitude, but have no idea how.

When I read the data using readOGR from the rgdal library I get the same coordinates but it tells me the following about what I assume is the projection (but the coordinates are the same)

Slot "proj4string":
CRS arguments:
+proj=utm +zone=32 +ellps=intl +units=m +no_defs

Reproducible example:
I've put the data for the example here:
https://github.com/sebastianbarfort/shapefiles

This should reproduce the problem:

library(maptools)
library(rgdal)
library(rgeos)

map = readShapePoly("~/Downloads/shapefiles-master/kommuner1983.shp")
centroid = gCentroid(map)
centroid


texts:
I'm trying to clip a csv-file containing xyz-data (points) by an irregular shaped polygon in ogr2ogr. My code runs without any errors or warnings, but it returns an exact copy of the input-csv. What am I doing wrong?

My code is:
ogr2ogr -f CSV output.csv -clipsrc clip_polygon.shp -clipsrclayer clip_polygon input.vrt

The vrt-file looks like:

<OGRVRTDataSource>
    <OGRVRTLayer name="test_layer">
       <SrcDataSource>input.csv</SrcDataSource>
       <GeometryType>wkbPoint</GeometryType>
       <LayerSRS>EPSG:25831</LayerSRS>
       <GeometryField encoding="PointFromColumns" x="field_1" y="field_2" z="field_3"/>
   </OGRVRTLayer>
</OGRVRTDataSource>


texts: Clipper: Split polygon layer with line layer - generates empty layer

texts:
I am using QGIS 2.14.3 Essen.

I wanted to split several polygons with several lines. Both layer are in GK4 coordinate system (used layers generated by: QGIS right-click on layer -> Save as... ) and they optically intersect. But when I use the clipper function to split the polygons by the line layer, I get an empty layer. Why?? Is there any possibility to get log data of the tool? Using the clipper function the other way round, works (getting the lines inside of polygon).

You can download my data here:
https://kommsafe.de/#/public/shares-downloads/Ek2uYFrP9Yo9DFiHPPSGi9WFVhbvMj4J


texts: ArcGIS Resample Tool Crashes with "Not implemented"

texts:
I'm trying to run the Resample GP tool with a mosaic dataset as an input. The mosaic dataset contains a raster for Peru, the Democratic Republic of Congo, and Indonesia. Each raster is in WGS84, which is projected to Web Mercator within the Mosaic Dataset. Cell size for each is 0.00025 degrees.

Each raster has a raster function applied to transform its value from Julian day (i.e. 1 - 366) to year + Julian day (1/1/2015 becomes 15001).

I'm running resample using the MAJORITY option, with output pixels at 55m x 55m in WGS 84 Web Mercator. Eight minutes into the process, it crashes with the classic Error 999999: Error executing function, followed by Not implemented and Failed to execute (Resample).

I'm used to error 999999, but have never seen Not Implemented. Is there an issue with my inputs? Or is this kind of operation even possible with a raster of this size and pixel depth? What workarounds are available?


texts: PostGIS query of point feature ~ Find nearest point to query location?

texts:
I am querying a point feature using a PostGIS request fired in JavaScript.  The decoded URL looks like this and returns great results for polygon features.

http://test.tablelocation.com/api/v1/sql/?q=select '+tempCoi+' from '+tempTableName+' where ST_Intersects(the_geom, ST_GeomFromText('POINT('+clickLng+' '+clickLat+')',4326))&format=json&diagnostics=true&_maxage=86400&callback=?

As you can see, there are a few dynamic Javascript variables that are injected into each URL request.  Two of these variables (clickLng and clickLat) specify the lat/lng of my query.  They are generated based on where the user clicks on a google map.  I am returning one attribute from the point layer referenced via the tempCoi variable.  When I try to query a point feature, no data is returned.

I am guessing that this is because the location of my query is too specific (i.e. 6+ decimal points).

To make a long question shorter, is there a way to fuzz my data request.

For example, if I click at 45.123456N and -113.654987W how could I look for the closest point to that query?

Hope this is clear...


texts:
Is there any way to automatically restart the AutoIncrement () function based on the Field Value on ArcMap?

I am developing a model in  ModelBuilder to name parcels automatically. The polygons are already sorted on the order they should be. I need to execute the AutoIncrement function on each of the Field Values separately.
To simplify the problem, the result I need is shown below. I already have the !Name! field and based on that, the !Count! field result must be acquired automatically.


texts: Get feature at specific coordinate with PyQgis

texts:
I am using QGIS 3.2 and I am trying to get the features at a specific coordinate for a specific layer. when I use mouse position event the following code works fine but when I calculate the coordinate and use Identify tool. It doesn't work.
I do some calculation in my program to get the coordinates, then that coordinate goes to following code to get the features but the following codes returns empty.

for lyr in [layer for layer in QgsProject.instance().mapLayers().values()]:
        if lyr.name() == 'equipment':
            LayerID = lyr
            print(lyr)
canvas=iface.mapCanvas()
f=QgsMapToolIdentify(canvas).identify(151.111650,-33.793260,[LayerID],QgsMapToolIdentify.TopDownAll)
print(f)

The following picture shows the coordinate of the the equipment but this code returns empty list

I have tried the following:

Increasing the search Radius in QGIS option (even 50.0mm returns empty)
give the input as QgsPoint

To make it clear I should explain that the coordinate not necessarily intersect with the feature coordinate(in most cases it Does).
I managed to get around the issue by writing the following code. It does what I need but the issue is. It is too Slow with big data set. Does QGIS have a better way to handle this?

def tool_identifier(self,layer_name,point_x,point_y,buffer=3):
    """ The only reason I created this is because I couldnt do it with maptool Identifier
    this tool gets two coordinate x and y and checks a layer if there is a feature in the buffer it will be add to output list

    :param Layer_name: Qgis Layer name should be Already open.
    :type Layer_name: string

    :param point_x: x coordinate of the center of buffer.
    :type point_x: float (151.117596)

    :param point_x: x coordinate of the center of buffer.
    :type point_x: float (-33.794810)

    :param buffer: this is the buffer radious
    :type buffer: it should be a number and default value is 3m

    """

    identified_feature_list=[]
    distance = QgsDistanceArea()
    for lyr in [layer for layer in QgsProject.instance().mapLayers().values()]:
        if lyr.name() == layer_name:layer_obj = lyr

    for f in layer_obj.getFeatures():
        distance = QgsDistanceArea()
        distance.setEllipsoid('WGS84')
        geom = f.geometry()
        current_point=geom.asPoint()
        m = distance.measureLine(QgsPointXY(current_point.x(),current_point.y()), QgsPointXY(point_x,point_y))
        if m<buffer:
            identified_feature_list.append(f)

    return identified_feature_list


texts:
I have more than 3000 points that spread in in a very big area. I need to find all the points (blue) that are in distance of 50 minute drive time, from the stores (green stars). I want to calculate it in automate way on all the points, and not each point separately.

Not all the stores located on the road layer (brown colored lines are from the Open Street Map).

I read

Creating drive time polygons using Open Source tools?
Calculate driving distance and time from 2 lat/lon points

but didn't found any solution.

I also tried to use the Road graph plugin-with no success.


texts: Help needed using REGEXP for address string parsing

texts:
Using JavaScript, my aim is to create an "intelligent" address search function with a single textbox (rather than separate textboxes for house number, street name, street type, suburb, etc).

I'm attempting to parse the input string entered by the user, to determine which components of the address have been entered. A complication is that users may enter combinations such as:

1 Kent Street Sydney (number, name, type, suburb)
1 King Georges Road Sydney (number, 2x name, type, suburb)
1 New South Head Road Rose Bay (number, 3x name, type, 2x suburb)
Kent Street Sydney (name, type, suburb)
Sydney (suburb)
2000 (postcode)
275/2/1000 (lot/section/plan)

At the moment I'm breaking the input string into tokens using string.split(" ") and iterating through each token to determine which component it represents.

This is working most of the time (it fails on 1 St Georges Road Sydney), but a colleague suggested I try using a REGEXP to improve the process. I can't see any examples of where this has been done for an address search.

Does anyone have any examples where REGEXP has been used for address validation?

(Note that I'm not actually entering this information into a geocoder at this stage - I just need to break the input string into its component parts, and know what each part represents)

Thanks,
Steve


texts: ArcGIS desktop map package, with clipping

texts:
Package Map (Data Management) almost does what I want: take everything in my current map, wrap it up in a nice little parcel with all geometry, attributes, labels, symbols, etc. etc. configured just the way that's needed into a single file that can be handed off to someone else or archived as a snapshot in project time.

The one thing the out of the box tool does not do which I need is clip. Clip everything to dataframe, clip to selected polygon (or graphic), clip to visible area (plus a margin), and then apply Package Map magic.

How might this be done?


texts: Raster flipping latitude or not showing shapefile on basemap

texts:
I am trying to save a 2D numpy array as a raster, and then split the raster using shapefiles. When I create the raster, I either have the problem of the latitude being flipped upside down, or that the shapefiles won't plot on the basemap.

Code - Shapefiles plot correctly on basemap, though latitude is flipped upside down.

# Create raster
xmin, ymin, xmax, ymax = [lon.min(), lat.min(), lon.max(), lat.max()]
nrows, ncols = np.shape(array)
geotransform = ([xmin, 0.1, 0, ymax, 0, -0.1])
output_raster = gdal.GetDriverByName('GTiff').Create('raster.tif', ncols, nrows, 1, gdal.GDT_Float32)
output_raster.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
output_raster.SetProjection( srs.ExportToWkt() )
output_raster.GetRasterBand(1).WriteArray(array)
output_raster.FlushCache()
output_raster = None

# Plot
ds = gdal.Open('raster.tif')
zz = ds.ReadAsArray()
gt = ds.GetGeoTransform()
x = np.linspace(gt[0], gt[0] + gt[1] * zz.shape[1], zz.shape[1])
y = np.linspace(gt[3], gt[3] + gt[5] * zz.shape[0], zz.shape[0])
xx, yy = np.meshgrid(x, y)
plt.figure(1)
plt.subplot(111)
plt.tight_layout()
ax1 = plt.gca()
m = Basemap(llcrnrlon=xx.min(), llcrnrlat=yy.min(), urcrnrlon=xx.max(), urcrnrlat=yy.max(), projection = 'cyl')
m.readshapefile('shapefile', 'states')
cs = m.contourf(xx, yy, zz)

Code change to make latitude correct, though shapefiles do not plot at all on basemap.

geotransform = ([xmin, 0.1, 0, ymax, 0, 0.1])

I am after both the correct latitude and the shapefiles to plot.

Can anyone help?

EDIT:
The next step in the script is clip the raster data with a shapefile, creating new rasters for each polygon (Code below). When I use the first method above these new rasters contain the data (though upside down), and when using the second method the data in them is all 0.

# Creates tif files for each state polygon
inraster = 'raster.tif'
inshape = 'shapefile'

ds = ogr.Open(inshape)
lyr = ds.GetLayer(0)

lyr.ResetReading()
ft = lyr.GetNextFeature()

while ft:

    country_name = ft.GetFieldAsString('NAME_1')

    outraster = inraster.replace('.tif', '_%s.tif' % country_name.replace(' ', '_'))
    subprocess.call(['gdalwarp', inraster, outraster, '-cutline', inshape,
                 '-crop_to_cutline', '-cwhere', "'NAME_1'='%s'" % country_name])

    ft = lyr.GetNextFeature()

ds = None


texts: Converting map to use space oblique mercator projection using ArcGIS for Desktop?

texts:
I have a map document in ArcMap which I want to convert to using Space oblique mercator projection, but I am not sure if ArcMap supports it.

Any hints on how to achieve this?


texts: Delete Unnecessary Fields In Standalone Table Using Arcpy

texts:
My question is similar to that I asked here. The script below, however, is intended to delete unnecessary fields from a single standalone table rather than a database full of tables. The answer in the question referenced above has not performed consistently and I am thinking narrowing it down to a single table may provide some insight to why this is so.

import arcpy
from arcpy import env

table_name = (r'C:\MyDatabase.gdb\StandAloneTable')
field_names = [f.name for f in arcpy.ListFields(table_name) if not f.required]
keep = True

with arcpy.da.SearchCursor(table_name, field_names) as cursor:
    for row in cursor:
        v = row[0]
        if v is None or isinstance(v, basestring) and \
           (v == '<Null>' or \
            v == 0 or \
            v == '' or \
            v == ' ' or \
            v == '  '):
            keep = False
            break

    if keep is False:
        arcpy.DeleteField_management(table_name, field_names)

This script runs without error aside from failing to delete the fields populated only with the values identified in first if-statement as I had intended. No fields are being deleted. Any ideas on what I am doing wrong?

Also, while I am testing this script on standalone tables within a file GDB, the standalone tables I plan to eventually run this script for are housed within an ArcSDE geodatabase.


texts: Will QGIS 2.99 and QGIS 3.0 always be the same?

texts:
For a dev, QGIS 2.99.0-Master is the next QGIS version which will be released.

that why http://doc.qgis.org/api/index.html is so important to make the bridge.

But it is simple as that ?

What will become 2.99 when QGIS 3 will be released ?
Is the dev version of QGIS 3.x will be named QGIS 3.99.0-Master ?


texts: Mismatch between OpenLayers 5 and XYZ-layer on some zoom-levels

texts:
I have an image pyramid that I created with gdal2tiles. Then I displayed this layer on top of another map in OpenLayers 5 with the XYZ source.

On the map zoom levels 0-8 the XYZ layer is not aligning correctly with the underlying map, in the next zoom level it fits perfectly. (See the attached images below).

My guess is that I missed something in how to set up the OpenLayers layer but I can't figure out what.

Here is the OpenLayers part for the layer:

 private xyzlayer: TileLayer = new TileLayer({
    opacity: 1,
    extent: [613934.99711258523166, 6724227.18296865839511, 618936.50444324919954, 6729085.88076259661466],
    source: new olXYZ({
      maxZoom: 8,
      projection: 'EPSG:3006',
      url: 'http://localhost:10/georeferencerApp/georeferenced/gavle_test/{z}/{x}/{-y}.png',
      tileGrid: new xyzTileGrid({
        extent: [613934.997113, 6724227.18297, 618936.504443, 6729085.88076],
        resolutions: [25.38836208458858, 12.69418104229429, 6.34709052114714, 3.17354526057357, 1.58677263028679, 0.79338631514339, 0.39669315757170, 0.19834657878585, 0.09917328939292]
       })
    })
});

This is the XML file from gdal2tiles:

<TileMap version="1.0.0" tilemapservice="http://tms.osgeo.org/1.0.0">
  <Title>gavle_1954.tif</Title>
  <Abstract></Abstract>
  <SRS>EPSG:3006</SRS>
  <BoundingBox minx="613934.99711258523166" miny="6724227.18296865839511" maxx="618936.50444324919954" maxy="6729085.88076259661466"/>
  <Origin x="613934.99711258523166" y="6724227.18296865839511"/>
  <TileFormat width="256" height="256" mime-type="image/png" extension="png"/>
  <TileSets profile="raster">
    <TileSet href="0" units-per-pixel="25.38836208458858" order="0"/>
    <TileSet href="1" units-per-pixel="12.69418104229429" order="1"/>
    <TileSet href="2" units-per-pixel="6.34709052114714" order="2"/>
    <TileSet href="3" units-per-pixel="3.17354526057357" order="3"/>
    <TileSet href="4" units-per-pixel="1.58677263028679" order="4"/>
    <TileSet href="5" units-per-pixel="0.79338631514339" order="5"/>
    <TileSet href="6" units-per-pixel="0.39669315757170" order="6"/>
    <TileSet href="7" units-per-pixel="0.19834657878585" order="7"/>
    <TileSet href="8" units-per-pixel="0.09917328939292" order="8"/>
  </TileSets>
</TileMap>

And this is the gdal2tiles command:

gdal2tiles -p raster -z 0-8 -s EPSG:3006 gavle_1954.tif gavle_test

Could anyone give me a hint on what the problem might be?


texts:
I'm having difficulty setting up an S3 bucket as a blobstore through geowebcache.

I have installed geoserver 2.8.2 which has gwc 1.8.1

I am following the docs at http://geowebcache.org/docs/1.8.x/configuration/storage.html#amazon-simple-storage-service-s3-blob-store

My config for the blobstore (with names and keys removed):

  <backendTimeout>120</backendTimeout>
<blobStores>
    <S3BlobStore default="false">
        <id>S3Cache</id>
        <enabled>false</enabled>
        <bucket>bucketName</bucket>
        <prefix>gwc</prefix>
        <awsAccessKey>key</awsAccessKey>
        <awsSecretKey>secretkey</awsSecretKey>
        <maxConnections>50</maxConnections>
        <useHTTPS>true</useHTTPS>
        <useGzip>true</useGzip>
     </S3BlobStore>
</blobStores>

I have geoserver installed through tomcat7 so when I restart the service and watch my geoserver logs I get the following error:

2016-03-02 07:56:22,294 INFO [config.XMLConfiguration] - Configuration file validated fine.

2016-03-02 07:56:22,381 WARN [support.DisposableBeanAdapter] - Invocation of destroy method failed on bean with name 'geoServerLoader': java.lang.NullPointerException

2016-03-02 07:56:22,383 ERROR [context.ContextLoader] - Context initialization failed
org.springframework.beans.factory.BeanCreationException: Error creating bean with name 'GeoSeverTileLayerCatalog' defined in URL [jar:file:/var/lib/tomcat7/webapps/geoserver/WEB-INF/lib/gs-gwc-2.8.2.jar!/geowebcache-geoserver-context.xml]: Cannot resolve reference to bean 'gwcXmlConfig' while setting constructor argument; nested exception is org.springframework.beans.factory.BeanCreationException: Error creating bean with name 'gwcXmlConfig' defined in URL [jar:file:/var/lib/tomcat7/webapps/geoserver/WEB-INF/lib/gs-gwc-2.8.2.jar!/geowebcache-core-context.xml]: Invocation of init method failed; nested exception is com.thoughtworks.xstream.converters.ConversionException: S3BlobStore : S3BlobStore

---- Debugging information ----

message             : S3BlobStore

cause-exception     : com.thoughtworks.xstream.mapper.CannotResolveClassException

cause-message       : S3BlobStore

class               : java.util.ArrayList

required-type       : java.util.ArrayList

converter-type      : com.thoughtworks.xstream.converters.collections.CollectionConverter

class[1]            : org.geowebcache.config.GeoWebCacheConfiguration
converter-type[1]   : com.thoughtworks.xstream.converters.reflection.ReflectionConverter

version : f22447f59ca37f72bfe6bc746d834b0e73e7c3fc/f22447f59ca37f72bfe6bc746d834b0e73e7c3fc

The error to myself looks like a module is missing, which it might be.  The documentation for the latest versions of geoserver and gwc indicated to me that s3BlobStore functionality is now included.  There used to be an s3 blob store community plugin, but I can't find the download for this, and it looked like this was simply a mod for the geoserver gui to setup an s3 blob where configuration could already be used.

Please could someone point me in the right direction, as to what I might be missing.

EDIT: I've installed Geowebcache 1.8.1 as a standalone, and the config loads correctly and the s3 blob shows in the blobstores.  Now if the integrated geowebcache is different to the standalone, how can you get the s3 functionality included?


texts: Unable to compile using make test in Linux

texts:
While trying to compile my test plugin I encounter this error, for the moment I can't figure out why, last week I had some issues on windows but eventually solved them via OSGeo Shell but this time what should I do?

I typed in prompt: make test and got:

pyrcc4 -o resources_rc.py  resources.qrc
make: pyrcc4: Command not found
make: *** [resources_rc.py] Error 127

From the looks of it, I do not have the make command installed?


texts: Repair/Check for empty geometry in ArcSDE?

texts: I understand that the standard validation/repair tools are not for use with ArcSDE and that ESRI needs to fix this, but are there work-arounds for at least checking where there might be empty geometry in a feature using ArcObjects or a script? Anything really. I'm attempting to export features from layers in an ArcSDE instance but some of those layers in the database give the "Attempted an operation on empty geometry" error when attempting to extract. I'm running ArcMap 10.3


texts: ArcObjects Export Field with x y gives InvalidCastException?

texts:
I am try to export field and X Y but I get the following error. Is it possible for someone to help me with this?

  system.InvalidCastException: Unable to cast COM object of type 'System.__ComObject' to interface type 'ESRI.ArcGIS.Geometry.IPoint'.
  This operation failed because the QueryInterface call on the COM component for the interface with IID
  '{1B077915-38C5-11D0-92D2-00805F7C28B0}' failed due to the following error: No such interface supported (Exception from HRESULT: 0x80004002 (E_NOINTERFACE)).

The Code:

Dim delimiter As String = ","
Dim aFileName As String
aFileName = "c:\log.txt"
Dim pPt As IPoint

FileOpen(1, aFileName, OpenMode.Output, OpenAccess.Default, OpenShare.Shared)

' Get first layer in map which is the annotation layer
Dim pMXDocument As IMxDocument
'  pMXDocument = My.ArcMap.Document.ActiveView
pMXDocument = My.ArcMap.Document

Dim pMap As IMap
pMap = pMXDocument.FocusMap
Dim pLayer As ILayer
pLayer = pMap.Layer(0)
Dim pFeatureLayer As IFeatureLayer
pFeatureLayer = pLayer

' Create a cursor
Dim pFtrCsr As IFeatureCursor
pFtrCsr = pFeatureLayer.Search(Nothing, False)

' loop printing text string
Dim s As String
Dim pFtr As IFeature
pFtr = pFtrCsr.NextFeature
Do While Not pFtr Is Nothing
   ' On Error Resume Next
    pPt = pFtr.Shape

    s = CStr(pFtr.Value(pFtr.Fields.FindField("TEXTSTRING"))) + "," + CStr(pPt.X) + "," + CStr(pPt.Y)

    '  s = pFeature.Value(pFeature.Fields.FindField("TEXTSTRING"))
    '  Debug.Print(s)

    WriteLine(1, s)

    pFtr = pFtrCsr.NextFeature
Loop

FileClose(1)


texts:
I need to export data from an ArcGIS 10.2 File Geodatabase (GDB) to an older version, so that users running an older version can access the data. Similarly, I would like to save a .mxd file to an older version, ideally 10 or 10.1.

It seems that most people want to go the other direction. Has anyone found a way to do this with access to only a 10.2 license?


texts: Create a multiple ring buffer according to field value + fixed set of values (ArcPy, Python)

texts:
I'm trying to find a code that will get the values from a
specific columm (width) of a lines feature dataset and sum it with a fixed set of values (10, 20 and 30 meters) in order to apply a Multiple Ring Buffer.

Attribute table:

The code sample below uses the fixed set of values (10, 20 and 30 meters) but not the width values of my lines shapefile.

import arcpy
arcpy.env.workspace = "C:/data/Py.gdb"
outFeatureClass = "c:/output/multibuffer"
distances = [10,20,30]
bufferUnit = "meters"
Field_Name = ""
dissolveType = "ALL"
arcpy.MultipleRingBuffer_analysis("lines", outFeatureClass, distances, bufferUnit, Field_Name, dissolveType)


texts: OpenLayers: draw path with only two points

texts:
I want to draw a path with no more then 2 points. On the second point I want the path to be finished.

I'm using the OpenLayers DrawFeature, I set a callbacks for "point" where I count the created points and if the points are two I'm finishing the sketch with finishSketch() method. The problem is that when finishSketch() method calls the finishGeometry() method which removes the last added point from the path which cause the following error to show

"Uncaught TypeError: Cannot read property 'geometry' of null "

My question is the any legit way for finishing a path programmaticaly or this is a bug in the DrawFeature control.


texts: Running TWI over 84 DEM files (.asc) using Whitebox R-interface

texts:
I have 84 1-m resolution DEM files (~1.2GB each saved as .asc files) across a large area (~10,000 square km) and would like to calculate TWI. For each tool (i.e. FillDepressions, Slope, D8FlowAccumulation and TWI) I am looping over sequential files one by one but some of my outputs are not picking up rivers and streams where they should be.
How does Whitebox deal with the edges when looping over multiple different files? Must I mosaic the tiles first to ensure hydrological connectivity between tiles?
rL <- list_files_with_exts("input_filepath", "tif", full.names = TRUE)
outFolder <- "output_filepath"

for(i in 1:length(rL)){

  rLs <- rL[i]

  outn <- paste(outFolder, gsub(".*/", '', rLs), "_", "sca", ".tif", sep = "")

  wbt_d8_flow_accumulation(
    input = rLs,
    output = outn,
    out_type = 'specific contributing area',
    log = FALSE,
    clip = FALSE,
    pntr = FALSE,
    esri_pntr = FALSE,
    wd = NULL,
    verbose_mode = FALSE
  )

  print(paste(i, length(rL)))

}


texts: Merging adjacent polygons in shapefile that has been split at tile boundaries?

texts:
I'm currently using building data from the Ordnance Survey Vectormap District and I've noticed that a lot of building polygons are split in 2 because of the way the OS split up the shapefile.

This is what I'm seeing in QGIS – I've highlighted some of the buildings split at the tile boundary:

I can't manually merge the buildings as there are way too many of them that have been split at the boundaries (this is just a zoomed in example). I'm looking for an automated way to solve this.

Unfortunately, each side of the split buildings has a different id so I can't easily dissolve them.

How would you suggest I can automatically merge the building tiles?

Edit

I'm now using PostGIS to merge the polygons split at the tile boundaries. Here is the SQL statement that does this for me – it's many, many times faster than doing it with QGIS:

DROP TABLE merged;
CREATE TABLE merged AS
SELECT
  -- Merge polygons that are within 0.01 metres of each other
  (ST_Dump(ST_Union(ST_Buffer(the_geom, 0.01)))).geom AS the_geom
FROM unmerged;

-- Update the geometry_columns table
SELECT Populate_Geometry_Columns();


texts: please can you help me with writing script for viewshed that evaluate a count of visible cells from every pixel? Or is there any plugin that can do it?


texts: What is the best way to merge lots of small adjacents polygons ? (POSTGIS)

texts:
I am looking for a way to merge lots of small adjacents polygons.

I have this kind of features:

So if my polygons are adjacents AND have the same color attribute then I want to merge those polygons. Typically I want to merge the selected features of the picture. I try to add a "ST_TOUCHES" condition but all the polygons are not necessarily connected to each other.

Another question on this website is similar to this one but don't take into account that polygons need to be connected.

So the first step should be something like:

 SELECT ST_Union(ST_SnapToGrid(the_geom,0.0001))
 FROM my_poly
 GROUP BY color;

But how to add the "need to be connected" condition ?

Any ideas ?

EDIT

I found a solution:

Step 1) Merge all the polygons that have the same color.

Step 2) Split multiparts polygons to singleparts polygons.

SELECT (ST_DUMP(ST_UNION(ST_SNAPTOGRID(the_geom,0.0001)))).geom, color
FROM my_poly
GROUP BY color

But there is perhaps a better solution, because i'm loosing some informations during the step 1).


texts:
I don't know if it is a bug or my lack of python skills, but I keep getting 'ERROR 000816 The tool is not valid' in my python script when I call ExportMetadata_conversion

Initially, I was attempting this on an sde database, but then I tried to take a step back just match the example in the documentation and use a file geodb:

import arcpy
from arcpy import env
env.workspace = "C:/data"
#set local variables
dir = arcpy.GetInstallInfo("desktop")["InstallDir"]
translator = dir + "Metadata/Translator/ESRI_ISO2ISO19139.xml"
arcpy.ExportMetadata_conversion ("data.gdb/roads", translator,
    "roads_19139.xml")

But I received the same error.

If anyone is successfully calling this method, does your code match the documentation?


texts: How to convert PostGis Table (Grid data) to Continuous field using Qgis?

texts: I'd like to transform my Gridded Layer (generated form a PostGIS table of points) to a continuous and colored layer. How could I do this with Qgis or any other program?


texts: R-Scatterplot: How to match palette() to not incrementing color-column

texts:
I want to visualize the feature space in my training areas using scatterplots. I use a nested classification schema, so my class values are not incrementing but sorted in categories like 100 (mainclass), 110 (subclass)... I have color values already specified for each class by the GIS I'm using to digitalize the training areas, so I want to have them in the scatterplot.

I already tried to update the palette with my color vector but the results are not satisfying at all. After some experimentation I noticed that for incrementing class values the result would be satisfying.

Here is a representative example:

b1 <- c(1:50)
b2 <- runif(b1, b1-10, b1+10)
val1 <- c(rep(1,10), rep(2,10), rep(3,10), rep(4,10), rep(5,10))
val2 <- c(rep(110,10), rep(120,10), rep(200,10), rep(300,10), rep(800,10))
df <- data.frame(b1, b2, val1, val2)

par(mfrow=c(2,2))

palette("default")
plot(df$b1, df$b2, col= df$val1, pch = 20)
plot(df$b1, df$b2, col= df$val2, pch = 20)

palette(c("#FF0000", "#FB5DF6", "#AF3C96", "#857EC7", "#E8D1D1"))
plot(df$b1, df$b2, col= df$val1, pch = 20)
plot(df$b1, df$b2, col= df$val2, pch = 20)

which will result in:

The plot at bottomleft would be satisfying. How can I match my palette to the not incrementing categories on the bottomright image column?


texts: Strange issue with Tomcat and AGS 10 .NET

texts:
I have one of my applications being deploying through Tomcat of Java WebADF.  It's being done by someone else, on another site.

It, as it would, worls perfectly on my box, but on their application, it fails.  Initially, all that displays is a empty page, with a Javascript error highlighted:

Message: 'EsriPageElement' is undefined

    Line: 17
    Char: 2
    Code: 0
    URI: http://Myserver:8080/MyApp/js/MA_startup.js

However, on investigating the Tomcat logs, I find this:

SEVERE: Unable to initialize resource 'ags1'. This resource will be removed
com.esri.adf.web.ags.ADFAGSException: Cannot create MapServer instance for server object 'MyService' [hosts=[MyServer]]
    at com.esri.adf.web.ags.data.AGSLocalConnection.activateResource(Unknown Source)

...and I am a bit perplexed.  The user is logging on as himself.  The only thing I can think of, is that this is permissions related, in that when Tomcat is fired up, it hasn't been started by anyone with the permissions to 'see' AGS?

Any ideas peeps?  Permissions?  IT Mice under the hood?  Too much alcohol abuse in my youth?  Am I bad and this is Karma?  I genuinely don't have a clue.


texts: Error when executing TIN Interpolation with two sets in QGIS

texts:
I am trying to TIN two datasets using TIN Interpolation with XYZ points in QGIS.
I keep getting this error:
Traceback (most recent call last):
File "/usr/share/qgis/python/plugins/processing/algs/qgis/TinInterpolation.py", line 188, in
processAlgorithm
writer.writeFile(feedback)
Exception: unknown

Execution failed after 135.41 seconds

Loading resulting layers
Algorithm 'TIN interpolation' finished

If I do the sets seperatly, the error only occurs on one of the sets.
What am I doing wrong? I checked the sets and the data is all valid.
If I take the set from a different year, but with the same geographical extent, samen problem.


texts: in Arcmap can I save a graphic line pointer within annotation?

texts:
I have an annotation layer for street names but some of the streets are very short and the annotation has been moved aside to reveal the underlying road, then on maps I draw a graphic line arrow pointing from the annotation text to the street.

Is there a way to save this line within the annotation so that I don't have to remember to draw it in on every map where i show street names?


texts: Layers not displaying in the correct order

texts: I have a shapefile that is in my Layer tree above an aerial imagery layer. Yet, the map covers the shapefile with the aerial imagery. Any ideas how to fix this?


texts: How to get an existing feature's rotation angle in OpenLayers 3

texts:
I created two external controls for my OpenLayers 3 map: Rotate Left and Rotate Right.
I already stored the last dragged/clicked feature in a global variable, yet I'm struggling to retrieve the actual rotation of the feature to "add" or "substract" degrees to it.

Here is what the code of my "Rotate Right" button looks like:

var evt = lastEvent;
if (evt == "") {
    return false;
} else {
    var coordX = evt.coordinate[0];
    var coordY = evt.coordinate[1];

    var feature = map.forEachFeatureAtPixel(evt.pixel,
        function(feature, layer) {
            return feature;
        });

    var layer = map.forEachFeatureAtPixel(evt.pixel,
        function(feature, layer) {
            return layer;
        });

    // This is where I need to get the feature's rotation value.
    var rotation = Math.random() * 10 + 1;
    var origin = new ol.geom.Point(coordX, coordY);
    feature.setStyle(new ol.style.Style({
        image: new ol.style.Icon( /** @type {olx.style.IconOptions} */ ({
            anchor: [0.5, 0.5],
            rotation: rotation,
            anchorXUnits: 'fraction',
            anchorYUnits: 'fraction',
            src: 'https://cdn4.iconfinder.com/data/icons/pictype-free-vector-icons/16/home-128.png'
        }))
    }));
}

I've been browsing through the API documentation and the examples but I can't seem to find a way to access the feature's rotation value once it's already drawn on the layer.


texts: Exporting NDVI result data as a table using Earth Engine

texts:
How can I calculate and make NDVI result for each shapefile in one horizontal line in Excel?

//

var Akdarya = ee.FeatureCollection("users/xodjayevshovkat/Akdarya_GEE_35");

Map.centerObject(Akdarya)
//Define a geometry.
Map.addLayer(Akdarya);

// Make the geometry a feature collection.

// Import S2 image collection.
var S2 = ee.ImageCollection('COPERNICUS/S2')
  // Filter collection.
  .filterDate('2018-05-01','2018-09-30')
  .filterBounds(Akdarya)
  .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 20);

// Define a cloud masking function.
var maskcloud1 = function(image) {
  var QA60 = image.select(['QA60']);
  return image.updateMask(QA60.lt(1))
};

//

// Define an index function (return only NDVI).
var NDVI = function(image) {
  return image.expression(
    '(NIR - RED) / (NIR + RED)',
    {
      'NIR': image.select('B8'),
      'RED': image.select('B4'),
    }).rename('NDVI').copyProperties(image, image.propertyNames());
};

// Apply the cloud mask function.
var NDVIcol = S2.map(maskcloud1)
 // Apply the NDVI transformation.
 .map(NDVI);

// Calculate the per-pixel median/mean NDVI for all images in collection.
var NDVImed = NDVIcol.mean();

// Define vis params.
var NDVIvis = {min: -0.001, max: 0.9, palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
    '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'};

//

// create a collection of the NDVI

var NDVIclip = NDVImed.clipToCollection(Akdarya)

// get the mean value for the region from each image
var regions = Akdarya.toList(Akdarya.size())
var ts = NDVIcol.map(function(image){
  var date = image.date().format("Y-M-D");
  var index = ee.String(Akdarya.get('Farmers'))
  var mean = image.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: Akdarya,
    scale: 30,
   });

// and return a feature with 'null' geometry with properties (dictionary)
  return ee.Feature(null, {'mean': mean.get('NDVI'),
                            'date': date})
});

// Export a .csv table of date, mean NDVI for watershed
Export.table.toDrive({
  collection: ts,
  description: 'NDVI_test4',
  fileFormat: 'CSV'
});


texts: Which DOS method use to convert at-sensor radiance to at-surface values in GRASS (landsat data, i.landsat.toar modul)?

texts:
I am trying to calculate surface reflectance values from raw Landsat DN data downloaded from GLOVIS. I have imported Landsat data in GRASS 6 with header by i.landsat.import and reprojected the scene into different PCS. Then I can convert DN to at-sensor radiance/reflectance values using i.landsat.toar.

Subsequently, I would like to convert at-sensor values to at-surface values (needed for further i.topo.corr minnaert correction because of the mountainous area). The simple metod is implemented by DOS[1-4] methods in i.landsat.toar. However, there are 4 DOS methods (http://grass.osgeo.org/grass65/manuals/i.landsat.toar.html):

DOS1: TAUv = 1.0, TAUz = 1.0 and Esky = 0.0
DOS2: TAUv = 1.0, Esky =
0.0, and TAUz = sin(e) for all bands with maximum wave length less than 1. (i.e. bands 4-6 MSS, 1-4 TM, and 1-4 ETM+) other bands TAUz =
1.0
DOS3: TAUv = exp[-t/cos(sat_zenith)], TAUz = exp[-t/sin(e)], Esky = rayleigh
DOS4: TAUv = exp[-t/cos(sat_zenith)], TAUz = exp[-t/sin(e)], Esky = PI * radiance_dark

Do you have some experiences, which one of DOS methods is the simpliest one and provides the most reliable results in mountainous area?
My area of interest is varying from 800-2500 a.s.l. in Central Europe, surface about 1050 km2 (only a fraction of whole Landsat scene).


texts: Mapping building to only one grid which contains more than half of its area using QGIS?

texts:
I have a shapefile of a city map and I have created vector grid as polygons.

After intersecting the two shapefiles, a building will map to all the grids it overlaps.

However, if I just want it to map to the one grid which contains the majority of the building (>50%), how could I achieve it?

Is there a function or plugin in gqis for this task? Or do I have to write my own code to do this?

I am new to QGIS.


texts:
I'm trying to override the tooltip style of Leaflet 1.3.1.and want to completely remove the tooltip pointer as shown in figure.

This is my css:

.ToolTipProvince{
          font-size: 15px;
                font-weight:bold;
                padding: 0px;
                background: none;
                border: none;
                color: aqua;
                text-shadow: 0px 0px 10px black;
                box-shadow: none;
            outline: none;
        }

and applying style to layer tooltip:

layer.bindTooltip("Province Number: "+feature.properties.PROVINCE+"", {permanent: true, direction: 'centre', className: 'ToolTipProvince'});


texts: Make rasters to show the extent of certain data values

texts:
I'm trying make a raster of the area between 5 and 10 thousand feet in Colorado. I have DEMs, and I'd like to make a raster where all the area in this range is a value of 1, and the area outside the range (<5000 or >10,000) is 0. I tried (lower value < "raster-name" < upper value) in raster calculator, but this did not work (it made files with no data- I have a post about this too, but at this point I think it'd be more helpful to just start over and try a different way)

How should I approach this? And if there's a way I can avoid analyzing the ~60 DEMs one at a time, that would be even better

This isn't answered by a similar post because they wanted one 'less than' bound and a 'does not equal a value' (they wanted less than or equal to 15 and doesn't equal 0). I'm looking for something that works with an upper AND a lower bound


texts:
I am looking for a basemap of LA county showing aerial imagery similar to Google Earth. However, I also need it in a specific projection: NAD_1983_StatePlane_California_V_FIPS_0405_Feet.

Does anyone know where I can find such a basemap?

I can't seem to find it in California's GIS database


texts: Transferring data from WMS/REST/SOAP to my GeoServer

texts:
I spent many sleepless and desperate hours trying to solve my problem. All links in google are now purple, no blue unread links. You guys are my last hope in solving my seems to be a simple problem.

Okay, I have an online app which is a customised module of Drupal 7. When a user types any address in the text-field the module converts address into long/lat and then performs multiple json calls to various servers and databases to collect various data. Basically, the user gets a PDF-report with something like:

the nearest school is located about 1km away from this house
crime rate in this area is low
property prices here are high
bla bla bla
public transport is good
...and this is an aerial view of your house.

This is where the problem starts. I am trying to place a JPG file provided by our government ArcGIS server into my report. The government server is a bit slow. Sometimes it produces JPG view fast, but about 20% of requests are being ignored. Even if I open that URL in 4 browser windows and hit ENTER 3 of 4 images come nice, but 4th will fail. Or not fail today, but will fail tomorrow. It is just unpredictable and unreliable.

So, I decided not to depend on our government server and to get a copy of all imagery on my server. I am interested in about 50x50km block of land (1 big city).

After searching around I came across GeoServer as a suitable platform to host imagery and provide my app with map export capabilities. I was also considering MapServer and Mapnik. So, not yet sure that GeoServer is the best choice for me.

Okay, I am running my site on AWS EC2 server (Ubuntu).

Now, the question is how can I extract lowest layer of Imagery from government ArcGIS server and import it into my GeoServer? Can I request it all if geo tiff? Should I design some sort of script to simulate fly-over and saving its imagery on the fly?

The government ArcGIS server offers REST / SOAP / WMS interfaces. I can run query EXPORT to say capture 2000x2000px image of any BBOX and save it, but there are too many blocks. I could script a loop with about 500m step to walk through the map up and down, but I believe there should be a better way of doing that.

I have a qGIS which is able to connect to the server via WMS and even save some tiff of the current view. But I couldn't the way to save a whole layer into geo-tiff or something.

There are some vector layers also available (borders of properties, roads, parks, etc) so I'd love to save these to my GeoServer too.

As for licensing, the imagery is provided on Creative Commons base.


texts: PyQGIS - Densify geometries given an interval

texts:
I have question, is it possible to make python command for QGIS function "Densify geometries given an interval":

input layer: Polygon.shp
output layer: Densify_interval.shp
interval: 50 meters

for. example...

processing.runalg("qgis:"DENSIFY_GIVEN_AN_INTERVAL ???",input, output", interval)


texts:
I have a shapefile with size 200 Mb and when I convert it to KML it creates a 4.3 GB KML file, which I am unable to load in fusion table of Google. So I want to break it into a smaller KML file.

How to create small KML file from that shapefile or postgres table (that I create from shapefile using shp2pgsql command?


texts: Returning value from other field, using FieldCalculator

texts:
I am very new to Python and need what I believe is a very basic script in order for my field to represent another fields' input, if the field's value is 0. I tried to search for this, but gave up after two hours... I am currently trying the script on two fields that I made only for this purpose, using ArcGIS 10.0.
So what I want to do is returning the value of field "TESTFIELDBEREGNER" in "TESTFIELD" IF TESTFIELD=0, else return TESTFIELD.
I believe the error may lie in the definition, as I don't know how to define strings in this correctly? The geoprocessing results window tells me there is a syntax error in line 2.

def RemoveNULLS(x, v):
  if x == '0'
    return v
  else:
    return x

 RemoveNULLS( !TESTFIELD!, !TESTFIELDBEREGNER!)


texts: Arcpy Slow to Read Datasets from Connection (.sde) File

texts:
I'm looking to automate publishing map services to ArcGIS Server, and in doing so I'm creating map documents that have only one layer; each layer is from a SQL Server 2014 database that I am accessing through an sde connection file. To connect to the database I am using this code, which I found here.

mxd = arcpy.mapping.MapDocument("C:/fullpath/to/blank.mxd")
arcpy.env.workspace = "C:/full/path/to/sdefile.sde"
dat = arcpy.ListFeatureclasses("layerName")[0]
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
lyr = arcpy.mapping.Layer(dat)
arcpy.mapping.AddLayer(df, lyr)
mxd.save()

This all works. The issue is that the line dat = arcpy.ListFeatureclasses("layerName")[0] takes almost twenty minutes to run, sometimes forty minutes. Why is it taking so long to read that data?

It is not the connection itself as I can read the layer fine when I connect to it from ArcMap (though this is slow, it draws on the order of seconds and not minutes). Note that I have tried this without the search string and it does not change the time. Edit: the search is very slow (greater than 20 minutes) in a database with 228 tables and 146 views. However, it seems the search is tolerable (1 minute) in a database with 278 tables but 0 views.

My guesses:
The database is not an SDE database, just a table with a geometry column. When I try to bring these layers into ArcMap, it shrieks that it does not have a unique identifier field and that it must calculate the spatial extent itself.

The tables that take the longest are views, so that may have something to do with it. But these layers do still publish, they just take forever. Another Edit: I tested searching for a table and searching for a view in the same database and the search times were nearly identical. This must have something to do with needing to search through all the data.

I think its possible I have same issue as this OP with their database #2: Adding layer to MXD from different SDE slows down ArcPy script. I'm using ArcGIS Desktop 10.4.

Any confirmation on my guesses?


texts: When I duplicate a print composer, item properties is blank

texts:
I need to have several, identically placed and scaled, print composers.
One with a grid, and the others without. Each showing different layers.
I set up one, added layers, then used duplicate to replicate it on a new composer, but when I do this, the item properties is blank and I can't select the map. This means I can't remove the grid.
I also found that, although I could refresh the view on the first duplicate, I could not do so on subsequent duplicates.
Is this a bug?
How can I set up my print composers?
I'm exporting them to another program to edit and then overlay so they have to show the exact same area, at the same scale, just with different info within.


texts: Reproject coordinates in order to correct them in QGIS

texts:
I got two type of coordinates(correct values and incorrects(15m error)), and i want to specify transformation parameters between them in QGIS. I thought that finding Helmert transformation parameter values via least squares method, which is supported by QGIS, and add them to proj4 could resolve my problem.
I wondering if i can use similarity/afine transformation in order to reproject my coordinates? Any suggestion that enable to me correct my coordinates?

How can i guess transformation parameters for affine transformation plugin?


texts: How to calculate % of area for each slope class in a tiff file?

texts:
I have a slope tiff layer. The pixel type is floating points. I categorized them in 3 classes based on the values. Now I need to know the number of cells in each class but this layer does not have a attribute table.

How can I calculate the % of the study area in the 3 classes?

I am using ArcGIS 10.3 I was trying to write a con statement using raster calculator but it has not worked yet.


texts:
Existence of fiducials marks in a digital camera file

I have created a generic digital camera file using Erdas 2015 with the following data in it:

LPS Camera File Version: 2.0
Camera Name: "generic digital camera"
Camera Description: "generic digital camera"
Focal Length [mm]: 35
Principal Point:
    x0 [mm]: 0
    y0 [mm]: 0
Fiducial Coordinates:
    Number of Fiducials: 4
      ID       X[mm]         Y[mm]
       1           30            0
       2          -30            0
       3            0           30
       4            0          -30
Radial Distortion Observations:
    Type: "Radial Distance"
    Number of Observations: 2
      ID   Distance[mm] Distortion[microns]
       1            0                1
       2            1                2
Radial Lens Distortion Coefficients:
    k0: 0
    k1: 0
    k2: 0
Extended Parameters:
    Type: "NONE"

All  included data makes sense except the existence of the fiducials marks , As I now the fiducials are existed for farm cameras not for digital ones. And when I click on the fiducials tap in the edrads it says that no fiducials are required, but when is save the file I found the following Fiducial Coordinates

Fiducial Coordinates:
Number of Fiducials: 4
ID       X[mm]         Y[mm]
1           30            0
2          -30            0
3            0           30
4            0          -30

Any Explanations for their existence please.

Edit:

“Digital sensors are fixed, so you don't need any mark on the image to define the coordinate system” Yes all references says the same thing that the Fiducials are not needed in case of digital cameras, But the issue here is Why the Erdas includes the above fiducials in case we saved a generic camera file and sometimes I receive camera files with the same fiducials (they doesn’t change even if the camera and focal length was changed. I couldn’t assume that these are the coordinates of the image corners simply because they doesn’t change even if the image size (based on the Sensor size :6.8 μ, 5.20 μ, …) was changed

This is weird and I could not figure out the purpose of include them in the digital camera file.


texts: Exporting to GeoTIFF creates unreasonable large files

texts:
I'm trying to produce a high quality image of a geological structural map provided for at this exact WMS link (The link includes the dimensions I wanted and the max possible resolution they provide)

http://maps.greenmin.gl/geusmap/ows/32624.jsp?nocache=nocache&whoami=[email]&LAYERS=grl_geus_500k_geology_map&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=true&STYLES=&SRS=EPSG%3A32624&BBOX=499980.0,7190220.0,609780.0,7400040.0&WIDTH=1500&HEIGHT=3000

I tried asking the question here, but had no luck with anyone providing me with an answer. I've also tried working it out like this post, without luck. So instead I've tried a different approach of using the Layer -> Save as function in QGIS.
Which is partially successful. I am able to output a geo referenced TIFF of the region with the desired dimensions. But for some reason that is beyond me, the file size is 92 gb and every program I try to open it in, chokes hard. This map is an exact matching overlay of a Sentinel-2 satellite image, in high quality, geo referenced and with a file size of about 1 gb. It ought to be possible for me to do something like that.

It would also be fine if I could get it without a georeference, I can add that manually. I just cannot fathom why it is so hard to have QGIS output a large dimensional image without producing ridiculous file sizes.

Has anyone dealt with anything like this and found a suitable solution or work around?


texts: Can we avoid using SIP library when writing standalone QGIS python script

texts:
I wanted to write a small custom python script using QGIS functionality. So I started with the QGIS documentation on the custom application section.

After having some troubles with the import qgis.core statement, which I could fix, I ran into the error message

ImportError: No module named sip

This one made me somewhat of a headache, as you could see here.

As I could not fix this one, I wonder why it's so much of a problem. Also, the doc states

  "PyQGIS bindings depend on SIP and PyQt4."

So I wonder, why it actually complains. Shouldn't it then come with the QGIS distribution? (Actually, there is a sip folder containing PyQT4 folder and also a sip.exe. But this does not seem to be of any effect.)

After all, my question is, is it possible to work around it, getting QGIS python running without SIP?


texts:
I have a similar problem, but I can't get the syntax correct. This refers to the following question:

Merge polygons in CartoDB?

I'm trying to merge multiple polygons into one (in this case, just the 5 NYC boroughs into a single polygon, but for now I'm trying to merge all of them). This doesn't seem to work, as the result is still separate polygons.

SELECT

counties_ny
, cartodb_id
, max(county) as county
, max(fips) as fips
, max(county_fip) as county_fip

, st_union(the_geom_webmercator) as the_geom_webmercator

FROM counties_ny

group by county, counties_ny.cartodb_id


texts: Zoom level issue with openlayers 5 and custom tile generation using gdal2tiles

texts:
Now i have changed my code but the output is still not what is expected, because the resolution provided is of 0 - 8 zoom but my minZoom starts from 12 itself,
will i get my tiles rendered ?? how exactly openlayers will understand that at a zoom of 12 it needs to fetch 0(zoom) resolution tile image or at zoom 13 it needs to fetch 1(zoom) resolution???

var resolutions = [156367.7919628329, 78183.89598141646, 39091.94799070823, 19545.973995354114, 9772.986997677057,
 4886.4934988385285,2443.2467494192642,1221.6233747096321,610.8116873548161];

 var extent  = [454019.0871110868, 36291.71352911883, 454248.1414938449,
  36399.381018930886];

 const map = new ol.Map({
      target: 'map',
      loadTilesWhileAnimating: true,
      controls: [],
      view: new ol.View({
         projection: ol.proj.get('EPSG:21096'),
         center: ol.proj.transform([32.6000, 0.3290], 'EPSG:4326', 'EPSG:21096'),
         resolutions: resolutions,
         zoom: 12
    })
});

const tileGrid = new ol.tilegrid.TileGrid({
     extent : extent,
     resolutions : resolutions})

   var xyzLayer = new ol.layer.Tile({
    source: new ol.source.XYZ({
            url: "http://localhost:8089/gis-web/map/tiles/{z}/{x}/{y}.png",
            tileGrid: tileGrid
    })
});

map.addLayer(xyzLayer);

 const setOlNonHierarchical = (layer) => {
 var vectorSource = new ol.source.Vector({
   format: new ol.format.GeoJSON(),
   loader: function(extent, resolution, projection) {
    showLoader(loader);
    let startIndex, maxFeatures, workspace;
    workspace = $('#prop_workspace').val();
    maxFeatures = $('#prop_max').val();
    startIndex = $('#prop_start_index').val();
    var proj = projection.getCode();
      var url = `${common_layers_url}${workspace}/${layer.geom}/${layer.title}/${maxFeatures}/${startIndex}`;
      var xhr = new XMLHttpRequest();
      xhr.open('GET', url);
      var onError = function() {
        vectorSource.removeLoadedExtent(extent);
      }
      xhr.onerror = onError;
      xhr.onload = function() {

        if (xhr.status == 200) {
          fadeLoader(loader);
          vectorSource.addFeatures(
              vectorSource.getFormat().readFeatures(xhr.responseText));
               //Set House No for Building
               if(layer.title=='Buildings'){
                      if(!isHouseNoSet)
                      setTimeout(setHouseNo,100);
               }
               //Open Draw Event for New Building Survey
               if(layer.title==$("#newBuildingSurveyId").attr('layer-text') && openSurveyNew==1){
                setTimeout(()=>{
                    selectedOpenEditFeatures();
                },100);
               }
        } else {
          onError();
        }
      }
      xhr.send();
    }
  });

 map.addLayer(new ol.layer.Vector({
       zIndex: layer.layerZindex,
       name_: layer.title,
       shortcode_:layer.layerShortCode,
       geomType: layer.geom,
       renderMode: 'image',
       customStyle: false,
       source:vectorSource
    }));

}

Okk so this the code from where my vector sources are getting added which show up at zoom 12


texts: QGIS : Weighted join attributes by location

texts:
In the example below, I have Urban Areas (UA1) and NUTS (NUTS1, NUTS2, NUTS3) regions in Europe.

Using the population of NUTS1, NUTS2 and NUTS3, I want to estimate the population of UA1.

I can join the population attribute of the NUTS region, but how does QGIS calculate the mean ? Is it a "basic" mean (suming the 3 populations and dividing by 3), or is it a weighted mean according to the surface of the urban area on each region ( in the example, it would be roughly 80% NUTS3 + 10% NUTS1 + 10% NUTS2)?

I could find this information on the QGIS documentation.


texts: Using advanced field calculator in QGIS Python console

texts:
I am trying to calculate area using field calculator in the Python console of QGIS.
Here is the expression which I wrote.
output8 = processing.runandload('qgis:advancedpythonfieldcalculator', "E:\GIS-Temp\Result\shapefile_2017\Clip1modified_2017.shp", "Area", 1, 5, 5, ' ', $geom.area(), "E:\GIS-Temp\Result\shapefile_2017\Clip2modified_2017.shp")

The error comes as syntax error.


texts: Clip image collection to imported shapefile (region of interest)?

texts:
I am working in Google Earth Engine. Is it possible to clip the Landsat scenes to the boundary of my area of interest rather than just filtering down the number of scenes to those that touch the area of interest?

I have figured out the filtering, but not how to clip. This is what my script looks like:

// Define a region of interest  (i.e. the shapefile you imported)
var roi = K2C_aoi;

// Load the Landsat 8 scaled radiance image collection.
var landsatCollection = ee.ImageCollection('LANDSAT/LC08/C01/T1')
    .filterDate('2019-06-10', '2019-07-12')
    .filterBounds(K2C_aoi)

// Make a cloud-free composite.
var composite = ee.Algorithms.Landsat.simpleComposite({
  collection: landsatCollection,
  asFloat: true
});

// Visualize the Composite using different spectral band combinations to help you identify object
Map.addLayer(composite, {bands: ['B4', 'B3', 'B2'], max: 0.5, gamma: 2}, 'Natural color', true);


texts: GRASS GIS v.surf.idw not working correctly

texts:
I have recently been struggling with an issue in GRASS. I am a new GRASS user and recently downloaded version 7.8.5 for my MacBook (macOS Catalina).  I converted this sample dataset to a shapefile. I am able to load and visualize the file and set the region to be the set of points. When I run v.surf.idw, I get the following error:
Input vector map <coal_test@PERMANENT> is 2D - using categories to interpolate
WARNING: No record for point (cat = 1)
WARNING: No record for point (cat = 2)
WARNING: No record for point (cat = 3)
.
.
.
WARNING: No record for point (cat = 96)
0 points loaded
ERROR: No points found

I have tried to use v.surf.idw with different datasets as well and am still getting the same error. I've also tried re-installing GRASS GIS to no avail.


texts: Using function to query text field in ArcGIS Select By Attributes?

texts:
I can do this:

And get result:

On the other hand, when I do this:

I get an error message:

Why query using function is working on numeric field and is not working on text in both shapefile and file geodatabase?


texts: Clustering adjacent polygons with values over a given threshold with Python

texts:
I have a GeoPandas dataframe that looks like the following
    geometry                                            raster_val
0   POLYGON ((1.85626 49.12003, 1.85626 49.11171, ...   28.530001
1   POLYGON ((1.87290 49.12003, 1.87290 49.11171, ...   28.930000
2   POLYGON ((1.88123 49.12003, 1.88123 49.11171, ...   28.670000
3   POLYGON ((1.89788 49.12003, 1.89788 49.11171, ...   28.510000
4   POLYGON ((1.90620 49.12003, 1.90620 49.11171, ...   29.969999

that I can show
f,ax=plt.subplots(figsize=(8,8))
gdf.plot(column='raster_val', scheme='quantiles', legend=True, ax=ax)

I would like to cluster all the adjacent polygons with a given values and create macropolygons based on such aggregration. Let say I want to create macropolygons that have raster_val > 33


texts: ArcGIS JS API 4.x - Polyline size(extrusion) in 3D View

texts: I've a PolyLine , I need to adjust its height and display it in SceneView (As a paper with 0 thickness) . I was able to display the PolyLine as a Line and I could be able to adjust its elevation (relative-to-ground) , but I couldn't control the visualVariables.size , Is there any special symbology - renderer I can use ?


texts: Splitting Polygons at midpoint using ArcPy?

texts:
I'm attempting to split around 4000 polygons at their midpoint, perpendicular to their longest axis (ie across the width at the midpoint), as in the diagram below.

Ideally, I'd like to do this automatically and avoid manually splitting each polygon. I have extracted the midpoint of polygons by converting the longest lines that can be drawn in each, I just need to determine a method to draw a width line across this point automatically.

Polygons vary in their width, and hence tools which split polygons by defining width lines of a certain length isn't really what I'm looking for.

Any ideas?


texts: DistanceTo Geometry printing generic error

texts:
Tried reading through everything else related to the distanceTo tool on SE, but im still running in to trouble.

Basically, my input Feature Class is a collection of polygons with begin/end lat/long.  There is one BEGIN_LAT,BEGIN_LON,END_LAT,END_LON for each record.  These lat/long pairs do not relate to the actual shape of the polygon (they might, but most will be different).  Trying to find the minimum distance value for each polygon to the beginning LAT/LON point and end LAT/LON point.  Eventually I am going to use the distance to compute time values, but for now I would just like to see the value print out properly :)

Code block as follows:

cursor = arcpy.da.SearchCursor(inFC,
['SHAPE@',"BEGIN_LAT","BEGIN_LON","END_LAT","END_LON"])

for row in cursor:
    polygonshape = row[0]
    beginpoint = arcpy.Point(row[2],row[1])
    begingeometry = arcpy.PointGeometry(beginpoint)
    distance = polygonshape.distanceTo(begingeometry)
    print distance

I'm getting the generic error of:

ValueError: <geoprocessing describe geometry object object at 0x05F17D40>

I also print the type of polygonshape and begingeometry....  Polygon and Point respectively.


texts: Is there a standard to construct a paginated API response for GeoJSON?

texts:
To return a paginated GeoJSON payload from an API, what is the best practice or most common convention for including page information?  The pagination applies to the number of features in a feature collection.  If the standards specify pagination, links to the details would be appreciated; AFAIK the pagination for JSON/GeoJSON is not specified in specs e.g.

https://tools.ietf.org/html/rfc7946 (geojson)
https://tools.ietf.org/html/rfc8142 (geojsons)

For example, consider these two options.  The first option is to add a "pages" data-structure to the GeoJSON, e.g.
{
    "type": "FeatureCollection",
    "features": [{
        "geometry": {
            "coordinates": [
                [
                    [-123.3948022, 41.254734],
                    [-123.3948022, 41.2537144],
                    [-123.3907336, 41.2537144],
                    [-123.3907336, 41.254734],
                    [-123.3948022, 41.254734]
                ]
            ],
            "type": "Polygon"
        },
        "properties": {
            "level": 0
        },
        "type": "Feature"
    }, {
        "geometry": {
            "coordinates": [
                [
                    [-123.3907336, 41.254734],
                    [-123.3907336, 41.2537144],
                    [-123.3893774, 41.2537144],
                    [-123.3893774, 41.254734],
                    [-123.3907336, 41.254734]
                ]
            ],
            "type": "Polygon"
        },
        "properties": {
            "level": 2
        },
        "type": "Feature"
    }],
    "pages": {
        "next": {
            "limit": 2,
            "offset": 2
        },
        "page": {
            "limit": 2,
            "offset": 0
        },
        "prev": null
    }
}

The second option creates a nested JSON payload with separate geojson and pages fields, i.e. the pseudo-JSON is like:
{
  "geojson": { ... },
  "pages": {
    "next": {
        "limit": 2,
        "offset": 2
    },
    "page": {
        "limit": 2,
        "offset": 0
    },
    "prev": null
  }
}

Or similarly, but without the extra "pages" nesting:
{
  "geojson": { ... },
  "next": {
        "limit": 2,
        "offset": 2
  },
  "page": {
        "limit": 2,
        "offset": 0
  },
  "prev": null
}


texts: What are the requirements to set "Supports Advanced Queries" true on a Feature Layer hosted on ArcGis Server 10.1

texts:
I'm trying to use the ArcGIS query feature, OrderByFields, but I've read that for that feature to work the feature layer must show "Supports Advanced Queries: true". The layer we are hosting says that "Supports Advanced Queries: false", but I can't find any documentation that says what is required for this property to be true.

We are currently using ArcMap 10.1 with shape files to publish our service definition files, and we are hosting the map services on ArcGIS server 10.1


texts: Can we use the strategies OpenLayers.Strategy.Fixed and OpenLayers.Strategy.BBOX simulataneously on same vector layer?


texts: Generate a list of all theme keywords in ArcCatalog?

texts:
Is it possible to generate a list of all themes keywords in ArcCatalog v9.3.1? similarly, I would like to do the same for place keywords as well?

My office only have about 10 standalone ArcMap licenses, without any databases or fancy servers, other than a NAS server.


texts: Getting ImportError no module named arcpy when importing into PythonWin?

texts:
Here is the error i get when i use the import arcpy command in pythonwin:

  Traceback (most recent call last):   File "", line
  1, in  ImportError: No module named arcpy

Can any one tell me how to fix this?


texts: How to render a static map from CartoDB?

texts:
CartoDB is working great for dynamic maps, but I have a use case to render a totally static map, just an image online. Google offers its static map api. Afaik there is no similar API for CartoDB. But does anyone have a work around that would let me display a single, non-zooming non-panning image?

It'd be ideal if I could just put it in an html img tag, but something with javascript could potentially work. It'd be ideal if it's lightweight, as many of these users have very slow connections.


texts: Merging Only Selected Lines in ArcGIS for Desktop?

texts:
Working in ArcGIS 10.1 for Desktop, is there a way to only merge selected lines from a layer?

For example, I want to merge these three red lines (all from one layer).


texts: ArcGIS Online app launcher display non-licensed apps

texts:
I invited a new member to my ArcGIS Online organization with an "Advanced" user type and "Publisher" role. The new member is not licensed with any field apps or the essential apps bundle.
But the member can see the experience builder, field maps(beta), and solutions in the ArcGIS Online app launcher.
So, my concern is if the member is not licensed with these apps, he/she should not be seeing these apps in the app launcher.


texts: Copying entire geodatabase with relationship classes and domains and reproject all feature classes from wgs84 to state plane?

texts:
I have an SDE geodatabase with several feature classes, tables, relationship classes, and domains.  All the feature classes are currently in WGS84 coordinate system.

I need an exact copy of the entire geodatabase, but with all feature classes projected in the appropriate state plane coordinate system.

I used the feature class to feature class tool to copy the feature classes from one gdb to another and set the environment output coordinate system to state plane, but I lost all of the relationships and tables.

Is there a way to do this in one step without having to recreate all the relationship classes?


texts:
I have set up a project in QGIS and can export maps using the Print Composer manually, but I want to have a single OSX Mac terminal command line that specifies the coordinates, the scale, and the dots per inch of a required map and it generates a PNG from my project?

Is this possible?


texts:
I want to import he5 irradiance files from giovanni. See a list here.
This post shows how to import he5 files into QGIS, but I'd like to import them into R. Is it possible? I tried package rhdf5 but couldn't find how to do it.


texts:
I have a geotiff file that I made in Python. After I created it, I realized it was for all of North America, and I only need the United States. I then got what I assumed was a shapefile of the US border, and tried to clip with

gdalwarp -cutline shapefile.shp input.tif output.tif

When I do this, it creates an output file, but the whole file is filled with zeros. I have tried a few other methods, but they all end up being filled with zeros. I think it could have to do with the spatial reference or the extent of either the shp or tif. Anyone have any ideas? Files included below.

Here is the geotiff file.

Here is the shapefile.

I have also tried a different shapefile, but that hasn't worked either.


texts:
I have installed PyCharm Edu and am trying to run a python script calling on ArcPy.

I get an import arcpy module error:

However, when I use the python scripter in ArcGIS Desktop, I see that it is able to find arcpy

I tried finding the actual arcpy folder in C:\Python27\ArcGIS10.3\Lib, but am unable to find it.

Why do I get the arcpy module error on Pycharm?

Is there an installation error in ArcGIS Desktop or PyCharm?


texts: Mapping a collection to extract monthly surface calculations to a table

texts:
I'm totally new using GEE and Javascript.
I'm trying to calculate a monthly based time series of a lagoon surface area. Besides, I want to export the results to a .csv or similar to compute some statistics outside GEE.

I think I'm on the right way, but I need some help to keep going.

So far, I can calculate the surface area for indivuald months whith this code:

Map.centerObject(area, 9);
Map.addLayer(area)

var thr = 0.35 // trheshold for water and non-water pixels

var h20mask = function(image) { // Water mask
  var water = image.gt(thr)  // convert to 1 those pixels above 'trheshold'
  return water.updateMask(water.gt(0));// mask pixels of non-water
};

var addMNDWI_L8 = function(image) {
  var mndwi = image.normalizedDifference(['B3', 'B6']).rename('MNDWI');
  return image.addBands(mndwi);
};

var addMNDWI_L5 = function(image) {
  var mndwi = image.normalizedDifference(['B2', 'B5']).rename('MNDWI');
  return image.addBands(mndwi);
};

var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
  .filterBounds(area)
  .filter(ee.Filter.lessThan('CLOUD_COVER',10));

var l5 = ee.ImageCollection('LANDSAT/LT05/C01/T2_TOA')
  .filterBounds(area)
  .filter(ee.Filter.lessThan('CLOUD_COVER',10));

var l8withMNDWI = l8.map(addMNDWI_L8);
var l5withMNDWI = l5.map(addMNDWI_L5);

// Surface area January 2018
var l8jan2018 = l8withMNDWI
  .filterDate('2018-01-01','2018-02-01')
  .mean()  // Just in case there is more than one image per month
  .select('MNDWI')
  .clip(area)

var water_l8jan2018 = h20mask(l8jan2018)

var waterArea = water_l8jan2018.reduceRegion({
   reducer: ee.Reducer.sum(),
   geometry: area,
   scale: 30,
   maxPixels: 1e13
});

print('January 2018')
print(waterArea);

In the case of Landsat 8, which is available from 2013-04-11 to 2020-01-11. I know I should map a months series using somethink like this:

var start = ee.Date.fromYMD(2013,04,01); // from 2013-04-01
var months = ee.List.sequence(0, 69); // to 2020-01-01
var startDates = months.map(function(d) {
  return start.advance(d, 'month');
});

But I don't know how to do it and neither how to build and export a table or matrix containing two columns, date and area.


texts:
i have a problem in ArcGIS with colormap. It is not active. I follow the steps written here http://desktop.arcgis.com/en/arcmap/10.3/manage-data/raster-and-images/creating-a-color-map-clr-file.htm, but in my case colormap is not active.

How can I make it active?


texts: Difference between Nominatim and Overpass-Api?

texts:
What is the difference between Nominatim and Overpass-Api?

What can't I do with Nominatim that I can do with the overpass-api ?


texts:
After the DEM creation, I wanted to create slope layer.

However, the result has not satisfied me. I don't think the slope map is correct. Have you ever encountered false calculated terrain map?

This is the DEM:

And this is the slope layer.


texts: ArcGIS Pro Unable to Pairwise Intersect Point with Polygon?

texts:
I am trying to create a point from XY coordinates then use Pairwise Intersect to determine the corresponding parcel. The intersect is not returning any values. Any thoughts on the error in this code?

parcel_tax_data is a simple polygon shapefile with GCS WGS 1984 coordinate system. It appears to successfully create the PointGeometry using the XY coordinates, convert to feature class and creates the intersect output file, but there is no contents. It does not recognize the point is intersected with any of of the parcels.

spatial_ref = arcpy.Describe(parcel_tax_data).spatialReference

point = arcpy.Point(lng_lat[0], lng_lat[1])
geo_point = arcpy.PointGeometry(point, spatial_ref)

arcpy.CreateFeatureclass_management(r"C:\path\to\project", "point_feature_class", "POINT", "", "DISABLED", "DISABLED", spatial_ref)

cursor = arcpy.da.InsertCursor("point_feature_class", ['SHAPE@XY'])
print("Cursor Created")
cursor.insertRow(geo_point)

inFeatures = ["point_feature_class", parcel_tax_data]
intersectOutput = "point"
print(inFeatures)

arcpy.PairwiseIntersect_analysis(inFeatures, intersectOutput, "ALL")


texts:
I'm trying to export a raster image from Earth Engine to Google Drive, however it cuts off the top and bottom part of the raster I'd like to have (the portion of the raster marked red in the image).

Here is the code I am using:

//Add MODIS Land Coverage dataset
var dataset = ee.ImageCollection('MODIS/006/MCD12Q1')
  .filterDate('2001-01-01', '2001-12-31');

var LandCover = dataset.select('LC_Prop2');

//Use bands 20, 25, 35 and 36: .eq(36).max(img.eq(35)).max(img.eq(25)).max(img.eq(30))
//Band 36: pure cropland, when cropland is lost, a change to 30/35/25
var LandCoverAgricultureMasked = LandCover.map(function(img) {
                   var mask = img.select(['LC_Prop2']).eq(36).max(img.eq(35)).max(img.eq(25));
                   return img.updateMask(mask)
                        .set('system:time_start', img.get('system:time_start'));
                   });

//Nice colors
var LandCoverVis = {
  min: 1.0,
  max: 40.0,
  palette: [
    '05450a', '086a10', '54a708', '78d203', '009900', 'c6b044', 'dcd159',
    'dade48', 'fbff13', 'b6ff05', '27ff87', 'c24f44', 'a5a5a5', 'ff6d4c',
    '69fff8', 'f9ffa4', '1c0dff', '05450a', '086a10', '54a708', '78d203',
    'dade48', 'fbff13', 'b6ff05', '27ff87', 'c24f44', 'a5a5a5', 'ff6d4c',
    '69fff8', 'f9ffa4', '1c0dff', '05450a', '086a10', '54a708', '78d203',
    'dade48', 'fbff13', 'b6ff05', '27ff87', '27ff87'
  ],
};

//Add fusion table to map
//var thesis_fusion = ee.FeatureCollection('ft:1Mn5wAIvtkzHzUMrG2yQryc9U4qSOa0ycJNEbywj6', 'geometry');
Map.addLayer(roi, {}, "Region of Interest");

//Clip Crop Coverage raster with region of interest (here: Object_ID XX)
var LandCoverAgricultureMaskedMosaic = LandCoverAgricultureMasked.mosaic()
var LandCoverAgricultureMaskedClip = LandCoverAgricultureMaskedMosaic.clip(roi);

Map.addLayer(LandCoverAgricultureMaskedClip, LandCoverVis, 'Agriculture');

//Time series of cropland change over all regions
//print(ui.Chart.image.seriesByRegion(LandCoverAgricultureMasked, thesis_fusion.filter(ee.Filter.eq("system:index", '130')), ee.Reducer.count()));

Export.image.toDrive({
  image: LandCoverAgricultureMaskedClip,
  description: 'Landcover Agriculture in ROI',
  scale: 500
});

Does anyone see the issue?


texts:
i am trying to pass polyline and polygon feature classes
i am not sure how to pass  the list of features to the Intersect object:

i tried this method:

ArrayList<Object> featuresList = new ArrayList<>();

        FeatureClass polygonfeatureClass = openFeatureClassFromShapeFile(srcShapefilePath, "polygon.shp");
        FeatureClass polylinefeatureClass = openFeatureClassFromShapeFile(srcShapefilePath, "lines.shp");

        Intersect newintersect = new Intersect();
        featuresList.add(polylinefeatureClass);
        featuresList.add(polygonfeatureClass);

        gp = new GeoProcessor();
        gp.setOverwriteOutput(true);
        gp.setTemporaryMapLayers(false);
        gp.setAddOutputsToMap(false);
        gp.setEnvironmentValue("workspace", srcShapefilePath);

newintersect.setInFeatures(featuresList)

i get this error :

Item not found in this collection. in 'DAO.Fields'

com.esri.arcgis.geoprocessing.GeoProcessor.execute(Unknown Source)

i tried to pass the shape files pathes as string too>
i got the same error
then i tried this method:
after i created 2 feature classes in gdb

String in1=srcShapefilePath + File.separator + "layers.gdb" + File.separator +            "lines";

       String in2=srcShapefilePath + File.separator +  "layers.gdb" + File.separator + "polygon";
       newintersect.setInFeatures(in1 + ";" + in2);
       gp.execute(newintersect, null);

and i got the same error

AutomationException: 0x80004005 - Unspecified error
    at com.esri.arcgis.geoprocessing.GeoProcessor.execute(Unknown Source)
    at com.esri.arcgis.geoprocessing.GeoProcessor.execute(Unknown Source)

Any Idea?


texts: Is there a way to quickly reduce precision in a number of DOUBLE fields?

texts:
I was wondering if there is a way to reduce the precision for display of a DOUBLE field of a number of fields, in a shapefile. I have 25 fields and would like to do it as automatically as possible but I'm not certain where to start. I don't particularly want to create new fields and transfer over the values manually. I am displaying the text using dynamic text on data driven pages . I just want to display it with two rather than 5-6 decimal places.

This is what the attribute table looks like, along with the relevant portion of layout:

Here is the field properties box:


texts:
I am using ArcMap 10.4.1 for desktop.
I have two layers with points in them and I am trying to get the points of layer A to "snap" to the closest points of layer B. I am using the Integrate tool with a tolerance value that should match the greatest distance between the points of the two layers. The tool will run fine but the end result is that the points in both layer A and B are changed rather than just the points in layer A snapping to layer B. Any idea what is going wrong here?


texts: overlap two GeoTiffs to fill missing information (transparency)

texts:
I have several hundered geotiff files. Some of them (almost all of them) have black areas (irregular) on the edges. However, these black areas are covered by some other file in the file set.

I read this post but it looks that the black in my file is not marked as NODATA (only 3 bands in the VRT file).

Is there away to join/overlap two geotiff files (into one) and make the black color as transparent so missing areas are filled from the other file?

Thanks.


texts: Representing projected coordinates from several UTM zones into the same one

texts:
I have a database with coordinates in WGS84 (4326) that all of them belong to Italy. Italy is split in two more or less equal UTM zones, 32 and 33. I need to transform those geographic coordinates to a projected system, but I don't want them to be in two different systems (ETRS89 UTM 32N AND ETRS89 UTM 33N), I want to keep them in a single representative SRID. My guess is that I should stick to something based on ETRS89.

I have the same problem with Mexico, which has several official UTM zones and I need to represent coordinates in just one. After a lot of searching, I found ITRF2008 LCC, but I don't know if that should be the one to go with.

So what should be the projected system to use for point representation and distance calculation for the whole country of Italy and Mexico?
Is there a method to make this kind of decisions when you run into an area with several equal UTM zones? Is there a consensus about this?

It is clear that when a country has most of its area in a single UTM zone (like Spain), you just pick the one that covers the biggest area like asked here. What about these two countries?

Also, this is for professional use, so the chosen system should be something somewhat "official" or "standard" if possible.


texts: How do I draw polygons from aerial imagery?

texts: I will be working on a project to automatically identify ponds from satellite imagery. I have uploaded the basemap into ArcMap 10.7, but I can't create a feature layer to start editing. How do I add a blank base layer? I will be tracing the ponds from the basemap (as this is the most detailed imagery I can find). Any other suggestions would be helpful. Could I use another layer instead of the basemap? This is for imagery in Zambia, so any US specific function will not work.


texts: Spatial Data Cluster in sql-server ? ( like ST_ClusterWithin in Postgis )

texts:
Postgres 9.5 providing ST_ClusterWithin for Cluster the spatial data in DB level.

Is there any available solution for data clustering in MS SQL (NOT SERVER CLUSTER)

OR

Is there Any Function ( like ST_ClusterWithin ) in MS SQL Server


texts:
I'd like to create a spatialite database and load it with shapefiles using the Mississippi Transverse Mercator (MSTM) projection, but it doesn't have an EPSG definition in the spatialite library of projections.
When loaded into a QGIS session, it correctly projects the shapefiles, automatically building a custom proj4 code for it from the ESRI prj files that come along with the data.  I checked that QGIS accurately reprojects this data, on the fly, to online WGS84 Orthophoto (Bing Maps) and it does.
I'm using qspatialite to load these shapefiles into a single newly created db.sqlite; I try to populate the mandatory SRID field with the custom projection code definition number that QGIS generates, for example SRID= "100003" etc, which is on the dropdown picklist to choose from in the qspatialite plugin.  Qspatialite can't load the shapefiles with any custom codes stored in QGIS though; it issues an error and fails to load the data--I figure because that isn't in the Spatialite EPSG library of projections.   To double check that the plugin installation works, I did load the shapefiles successfully using a standard EPSG number, but my shapefiles are really in a custom projection CRS.
Is there a way I can copy/paste QGIS's custom proj4 definition to the spatialite EPSG library and create a brand new custom SRID definition, then load my data using that new SRID number in the mandatory data loader field?  Or, would spatialite-gui be better than qspatialite?
These are the MSTM projection parameters (not Proj4):
Scale Factor: 0.9998335
Central Meridian: -89 45' 00''
Central Latitude: 32 30' 00''
False Northing: 1,300,000 meters
False Easting: 500,000 meters
Projection Units: meters
Spheroid: GRS80
Datum: NAD83

--EDIT OF POST 1/24/2013 BELOW--
I just wanted to add a workaround method for bulk-batch importing shapefiles into a spatialite database, using Qspatialite's bulk load OGR-file button. This would be useful if the IMPORT OGR FILE command won't accept a custom EPSG number that you need to use.   For myself, I just entered a similar but commonly used EPSG number using that also uses the same DATUM.  For example, my custom EPSG: 3814 apparently gets little use and hasn't yet been included in QGIS spatial reference database (Lisboa 1.8.0).  I found that the "custom" coordinate system uses NAD83.  Instead of entering the SRID of 3814 that causes an import failure, I just used a similar NAD83, in this case EPSG:26995, that's preloaded into QGIS.  They import just fine like that.  Of course the imported shapefiles will not be projected correctly in QGIS relative to a basemap.
To fix that I created a custom coordinate system in QGIS using the +proj code that QGIS generates to correctly project the individual shapefiles that are now importing to a spatialite database, that also included a 3814.prj spatial reference file.  THEN after creating that CUSTOM CRS, I was able to load the files from my spatialite database and individually set the projection under the  of each layer to my custom spatial reference code that I saved in QGIS.  Presto, everything aligns perfectly with my basemap while using spatialite layers.  And, I was able to use Qspatialite to bulk-batch over a hundred shapefiles and then a couple of hundred dbf tables of additional information.  One by one would have been time consuming.


texts: Assign buffer an attribute based on the largest polygon inside

texts:
I would like to assign a buffer polygon a value from another layer. I have cliped the two layers together. I then tried to dissolve the new layer by buffer ID, but I could not manage to get the largest area inside to be the value for the whole buffer.

I tried to do it by using calculate field, but im not very good at python.


texts: Carto Builder is enabled but I'm still seeing the old editor

texts:
After taking a break from mapping, I've just logged back into Carto. My Carto account has builder enabled, but I'm still seeing the old editor when I look at my datasets and maps? The profile account dashboard has been updated but not the editor?

How do I fix this? I've tried creating and publishing a new map to see if that would "activate" the new builder but it hasn't.

Note: I created the account years ago when it was CartoDB

Suggestions?